1: #region Translated by Jose Antonio De Santiago-Castillo.
2:
3: //Translated by Jose Antonio De Santiago-Castillo.
4: //E-mail:JAntonioDeSantiago@gmail.com
5: //Web: www.DotNumerics.com
6: //
7: //Fortran to C# Translation.
8: //Translated by:
9: //F2CSharp Version 0.71 (November 10, 2009)
10: //Code Optimizations: None
11: //
12: #endregion
13:
14: using System;
15: using DotNumerics.FortranLibrary;
16:
17: namespace DotNumerics.CSLapack
18: {
19: /// <summary>
20: /// -- LAPACK driver routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
27: /// symmetric tridiagonal matrix using the divide and conquer method.
28: /// The eigenvectors of a full or band real symmetric matrix can also be
29: /// found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
30: /// matrix to tridiagonal form.
31: ///
32: /// This code makes very mild assumptions about floating point
33: /// arithmetic. It will work on machines with a guard digit in
34: /// add/subtract, or on those binary machines without guard digits
35: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
36: /// It could conceivably fail on hexadecimal or decimal machines
37: /// without guard digits, but we know of none. See DLAED3 for details.
38: ///
39: ///</summary>
40: public class DSTEDC
41: {
42:
43:
44: #region Dependencies
45:
46: LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANST _dlanst; DGEMM _dgemm; DLACPY _dlacpy; DLAED0 _dlaed0;
47: DLASCL _dlascl;DLASET _dlaset; DLASRT _dlasrt; DSTEQR _dsteqr; DSTERF _dsterf; DSWAP _dswap; XERBLA _xerbla;
48:
49: #endregion
50:
51:
52: #region Fields
53:
54: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; bool LQUERY = false; int FINISH = 0;
55: int I = 0;int ICOMPZ = 0; int II = 0; int J = 0; int K = 0; int LGN = 0; int LIWMIN = 0; int LWMIN = 0; int M = 0;
56: int SMLSIZ = 0;int START = 0; int STOREZ = 0; int STRTRW = 0; double EPS = 0; double ORGNRM = 0; double P = 0;
57: double TINY = 0;
58:
59: #endregion
60:
61: public DSTEDC(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANST dlanst, DGEMM dgemm, DLACPY dlacpy, DLAED0 dlaed0, DLASCL dlascl, DLASET dlaset, DLASRT dlasrt
62: , DSTEQR dsteqr, DSTERF dsterf, DSWAP dswap, XERBLA xerbla)
63: {
64:
65:
66: #region Set Dependencies
67:
68: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dgemm = dgemm;
69: this._dlacpy = dlacpy;this._dlaed0 = dlaed0; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasrt = dlasrt;
70: this._dsteqr = dsteqr;this._dsterf = dsterf; this._dswap = dswap; this._xerbla = xerbla;
71:
72: #endregion
73:
74: }
75:
76: public DSTEDC()
77: {
78:
79:
80: #region Dependencies (Initialization)
81:
82: LSAME lsame = new LSAME();
83: IEEECK ieeeck = new IEEECK();
84: IPARMQ iparmq = new IPARMQ();
85: DLAMC3 dlamc3 = new DLAMC3();
86: DLASSQ dlassq = new DLASSQ();
87: XERBLA xerbla = new XERBLA();
88: DCOPY dcopy = new DCOPY();
89: IDAMAX idamax = new IDAMAX();
90: DLAPY2 dlapy2 = new DLAPY2();
91: DLAMRG dlamrg = new DLAMRG();
92: DROT drot = new DROT();
93: DSCAL dscal = new DSCAL();
94: DNRM2 dnrm2 = new DNRM2();
95: DLAED5 dlaed5 = new DLAED5();
96: DLAE2 dlae2 = new DLAE2();
97: DLAEV2 dlaev2 = new DLAEV2();
98: DSWAP dswap = new DSWAP();
99: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
100: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
101: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
102: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
103: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
104: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
105: DLANST dlanst = new DLANST(lsame, dlassq);
106: DGEMM dgemm = new DGEMM(lsame, xerbla);
107: DLACPY dlacpy = new DLACPY(lsame);
108: DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
109: DLAED6 dlaed6 = new DLAED6(dlamch);
110: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
111: DLASET dlaset = new DLASET(lsame);
112: DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
113: DLAED1 dlaed1 = new DLAED1(dcopy, dlaed2, dlaed3, dlamrg, xerbla);
114: DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
115: DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
116: DGEMV dgemv = new DGEMV(lsame, xerbla);
117: DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
118: DLAED7 dlaed7 = new DLAED7(dgemm, dlaed8, dlaed9, dlaeda, dlamrg, xerbla);
119: DLARTG dlartg = new DLARTG(dlamch);
120: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
121: DLASR dlasr = new DLASR(lsame, xerbla);
122: DLASRT dlasrt = new DLASRT(lsame, xerbla);
123: DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
124: , dlasrt, dswap, xerbla);
125: DLAED0 dlaed0 = new DLAED0(dcopy, dgemm, dlacpy, dlaed1, dlaed7, dsteqr, xerbla, ilaenv);
126: DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
127:
128: #endregion
129:
130:
131: #region Set Dependencies
132:
133: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dgemm = dgemm;
134: this._dlacpy = dlacpy;this._dlaed0 = dlaed0; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasrt = dlasrt;
135: this._dsteqr = dsteqr;this._dsterf = dsterf; this._dswap = dswap; this._xerbla = xerbla;
136:
137: #endregion
138:
139: }
140: /// <summary>
141: /// Purpose
142: /// =======
143: ///
144: /// DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
145: /// symmetric tridiagonal matrix using the divide and conquer method.
146: /// The eigenvectors of a full or band real symmetric matrix can also be
147: /// found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
148: /// matrix to tridiagonal form.
149: ///
150: /// This code makes very mild assumptions about floating point
151: /// arithmetic. It will work on machines with a guard digit in
152: /// add/subtract, or on those binary machines without guard digits
153: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
154: /// It could conceivably fail on hexadecimal or decimal machines
155: /// without guard digits, but we know of none. See DLAED3 for details.
156: ///
157: ///</summary>
158: /// <param name="COMPZ">
159: /// (input) CHARACTER*1
160: /// = 'N': Compute eigenvalues only.
161: /// = 'I': Compute eigenvectors of tridiagonal matrix also.
162: /// = 'V': Compute eigenvectors of original dense symmetric
163: /// matrix also. On entry, Z contains the orthogonal
164: /// matrix used to reduce the original matrix to
165: /// tridiagonal form.
166: ///</param>
167: /// <param name="N">
168: /// (input) INTEGER
169: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
170: ///</param>
171: /// <param name="D">
172: /// (input/output) DOUBLE PRECISION array, dimension (N)
173: /// On entry, the diagonal elements of the tridiagonal matrix.
174: /// On exit, if INFO = 0, the eigenvalues in ascending order.
175: ///</param>
176: /// <param name="E">
177: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
178: /// On entry, the subdiagonal elements of the tridiagonal matrix.
179: /// On exit, E has been destroyed.
180: ///</param>
181: /// <param name="Z">
182: /// (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
183: /// On entry, if COMPZ = 'V', then Z contains the orthogonal
184: /// matrix used in the reduction to tridiagonal form.
185: /// On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
186: /// orthonormal eigenvectors of the original symmetric matrix,
187: /// and if COMPZ = 'I', Z contains the orthonormal eigenvectors
188: /// of the symmetric tridiagonal matrix.
189: /// If COMPZ = 'N', then Z is not referenced.
190: ///</param>
191: /// <param name="LDZ">
192: /// (input) INTEGER
193: /// The leading dimension of the array Z. LDZ .GE. 1.
194: /// If eigenvectors are desired, then LDZ .GE. max(1,N).
195: ///</param>
196: /// <param name="WORK">
197: /// (workspace/output) DOUBLE PRECISION array,
198: /// dimension (LWORK)
199: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
200: ///</param>
201: /// <param name="LWORK">
202: /// (input) INTEGER
203: /// The dimension of the array WORK.
204: /// If COMPZ = 'N' or N .LE. 1 then LWORK must be at least 1.
205: /// If COMPZ = 'V' and N .GT. 1 then LWORK must be at least
206: /// ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
207: /// where lg( N ) = smallest integer k such
208: /// that 2**k .GE. N.
209: /// If COMPZ = 'I' and N .GT. 1 then LWORK must be at least
210: /// ( 1 + 4*N + N**2 ).
211: /// Note that for COMPZ = 'I' or 'V', then if N is less than or
212: /// equal to the minimum divide size, usually 25, then LWORK need
213: /// only be max(1,2*(N-1)).
214: ///
215: /// If LWORK = -1, then a workspace query is assumed; the routine
216: /// only calculates the optimal size of the WORK array, returns
217: /// this value as the first entry of the WORK array, and no error
218: /// message related to LWORK is issued by XERBLA.
219: ///</param>
220: /// <param name="IWORK">
221: /// (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
222: /// On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
223: ///</param>
224: /// <param name="LIWORK">
225: /// (input) INTEGER
226: /// The dimension of the array IWORK.
227: /// If COMPZ = 'N' or N .LE. 1 then LIWORK must be at least 1.
228: /// If COMPZ = 'V' and N .GT. 1 then LIWORK must be at least
229: /// ( 6 + 6*N + 5*N*lg N ).
230: /// If COMPZ = 'I' and N .GT. 1 then LIWORK must be at least
231: /// ( 3 + 5*N ).
232: /// Note that for COMPZ = 'I' or 'V', then if N is less than or
233: /// equal to the minimum divide size, usually 25, then LIWORK
234: /// need only be 1.
235: ///
236: /// If LIWORK = -1, then a workspace query is assumed; the
237: /// routine only calculates the optimal size of the IWORK array,
238: /// returns this value as the first entry of the IWORK array, and
239: /// no error message related to LIWORK is issued by XERBLA.
240: ///</param>
241: /// <param name="INFO">
242: /// (output) INTEGER
243: /// = 0: successful exit.
244: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
245: /// .GT. 0: The algorithm failed to compute an eigenvalue while
246: /// working on the submatrix lying in rows and columns
247: /// INFO/(N+1) through mod(INFO,N+1).
248: ///</param>
249: public void Run(string COMPZ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Z, int offset_z, int LDZ
250: , ref double[] WORK, int offset_work, int LWORK, ref int[] IWORK, int offset_iwork, int LIWORK, ref int INFO)
251: {
252:
253: #region Array Index Correction
254:
255: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_z = -1 - LDZ + offset_z; int o_work = -1 + offset_work;
256: int o_iwork = -1 + offset_iwork;
257:
258: #endregion
259:
260:
261: #region Strings
262:
263: COMPZ = COMPZ.Substring(0, 1);
264:
265: #endregion
266:
267:
268: #region Prolog
269:
270: // *
271: // * -- LAPACK driver routine (version 3.1) --
272: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
273: // * November 2006
274: // *
275: // * .. Scalar Arguments ..
276: // * ..
277: // * .. Array Arguments ..
278: // * ..
279: // *
280: // * Purpose
281: // * =======
282: // *
283: // * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
284: // * symmetric tridiagonal matrix using the divide and conquer method.
285: // * The eigenvectors of a full or band real symmetric matrix can also be
286: // * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
287: // * matrix to tridiagonal form.
288: // *
289: // * This code makes very mild assumptions about floating point
290: // * arithmetic. It will work on machines with a guard digit in
291: // * add/subtract, or on those binary machines without guard digits
292: // * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
293: // * It could conceivably fail on hexadecimal or decimal machines
294: // * without guard digits, but we know of none. See DLAED3 for details.
295: // *
296: // * Arguments
297: // * =========
298: // *
299: // * COMPZ (input) CHARACTER*1
300: // * = 'N': Compute eigenvalues only.
301: // * = 'I': Compute eigenvectors of tridiagonal matrix also.
302: // * = 'V': Compute eigenvectors of original dense symmetric
303: // * matrix also. On entry, Z contains the orthogonal
304: // * matrix used to reduce the original matrix to
305: // * tridiagonal form.
306: // *
307: // * N (input) INTEGER
308: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
309: // *
310: // * D (input/output) DOUBLE PRECISION array, dimension (N)
311: // * On entry, the diagonal elements of the tridiagonal matrix.
312: // * On exit, if INFO = 0, the eigenvalues in ascending order.
313: // *
314: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
315: // * On entry, the subdiagonal elements of the tridiagonal matrix.
316: // * On exit, E has been destroyed.
317: // *
318: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
319: // * On entry, if COMPZ = 'V', then Z contains the orthogonal
320: // * matrix used in the reduction to tridiagonal form.
321: // * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
322: // * orthonormal eigenvectors of the original symmetric matrix,
323: // * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
324: // * of the symmetric tridiagonal matrix.
325: // * If COMPZ = 'N', then Z is not referenced.
326: // *
327: // * LDZ (input) INTEGER
328: // * The leading dimension of the array Z. LDZ >= 1.
329: // * If eigenvectors are desired, then LDZ >= max(1,N).
330: // *
331: // * WORK (workspace/output) DOUBLE PRECISION array,
332: // * dimension (LWORK)
333: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
334: // *
335: // * LWORK (input) INTEGER
336: // * The dimension of the array WORK.
337: // * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
338: // * If COMPZ = 'V' and N > 1 then LWORK must be at least
339: // * ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
340: // * where lg( N ) = smallest integer k such
341: // * that 2**k >= N.
342: // * If COMPZ = 'I' and N > 1 then LWORK must be at least
343: // * ( 1 + 4*N + N**2 ).
344: // * Note that for COMPZ = 'I' or 'V', then if N is less than or
345: // * equal to the minimum divide size, usually 25, then LWORK need
346: // * only be max(1,2*(N-1)).
347: // *
348: // * If LWORK = -1, then a workspace query is assumed; the routine
349: // * only calculates the optimal size of the WORK array, returns
350: // * this value as the first entry of the WORK array, and no error
351: // * message related to LWORK is issued by XERBLA.
352: // *
353: // * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
354: // * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
355: // *
356: // * LIWORK (input) INTEGER
357: // * The dimension of the array IWORK.
358: // * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
359: // * If COMPZ = 'V' and N > 1 then LIWORK must be at least
360: // * ( 6 + 6*N + 5*N*lg N ).
361: // * If COMPZ = 'I' and N > 1 then LIWORK must be at least
362: // * ( 3 + 5*N ).
363: // * Note that for COMPZ = 'I' or 'V', then if N is less than or
364: // * equal to the minimum divide size, usually 25, then LIWORK
365: // * need only be 1.
366: // *
367: // * If LIWORK = -1, then a workspace query is assumed; the
368: // * routine only calculates the optimal size of the IWORK array,
369: // * returns this value as the first entry of the IWORK array, and
370: // * no error message related to LIWORK is issued by XERBLA.
371: // *
372: // * INFO (output) INTEGER
373: // * = 0: successful exit.
374: // * < 0: if INFO = -i, the i-th argument had an illegal value.
375: // * > 0: The algorithm failed to compute an eigenvalue while
376: // * working on the submatrix lying in rows and columns
377: // * INFO/(N+1) through mod(INFO,N+1).
378: // *
379: // * Further Details
380: // * ===============
381: // *
382: // * Based on contributions by
383: // * Jeff Rutter, Computer Science Division, University of California
384: // * at Berkeley, USA
385: // * Modified by Francoise Tisseur, University of Tennessee.
386: // *
387: // * =====================================================================
388: // *
389: // * .. Parameters ..
390: // * ..
391: // * .. Local Scalars ..
392: // * ..
393: // * .. External Functions ..
394: // * ..
395: // * .. External Subroutines ..
396: // * ..
397: // * .. Intrinsic Functions ..
398: // INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT;
399: // * ..
400: // * .. Executable Statements ..
401: // *
402: // * Test the input parameters.
403: // *
404:
405: #endregion
406:
407:
408: #region Body
409:
410: INFO = 0;
411: LQUERY = (LWORK == - 1 || LIWORK == - 1);
412: // *
413: if (this._lsame.Run(COMPZ, "N"))
414: {
415: ICOMPZ = 0;
416: }
417: else
418: {
419: if (this._lsame.Run(COMPZ, "V"))
420: {
421: ICOMPZ = 1;
422: }
423: else
424: {
425: if (this._lsame.Run(COMPZ, "I"))
426: {
427: ICOMPZ = 2;
428: }
429: else
430: {
431: ICOMPZ = - 1;
432: }
433: }
434: }
435: if (ICOMPZ < 0)
436: {
437: INFO = - 1;
438: }
439: else
440: {
441: if (N < 0)
442: {
443: INFO = - 2;
444: }
445: else
446: {
447: if ((LDZ < 1) || (ICOMPZ > 0 && LDZ < Math.Max(1, N)))
448: {
449: INFO = - 6;
450: }
451: }
452: }
453: // *
454: if (INFO == 0)
455: {
456: // *
457: // * Compute the workspace requirements
458: // *
459: SMLSIZ = this._ilaenv.Run(9, "DSTEDC", " ", 0, 0, 0, 0);
460: if (N <= 1 || ICOMPZ == 0)
461: {
462: LIWMIN = 1;
463: LWMIN = 1;
464: }
465: else
466: {
467: if (N <= SMLSIZ)
468: {
469: LIWMIN = 1;
470: LWMIN = 2 * (N - 1);
471: }
472: else
473: {
474: LGN = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N)) / Math.Log(TWO)));
475: if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
476: if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
477: if (ICOMPZ == 1)
478: {
479: LWMIN = 1 + 3 * N + 2 * N * LGN + 3 * (int)Math.Pow(N, 2);
480: LIWMIN = 6 + 6 * N + 5 * N * LGN;
481: }
482: else
483: {
484: if (ICOMPZ == 2)
485: {
486: LWMIN = 1 + 4 * N + (int)Math.Pow(N, 2);
487: LIWMIN = 3 + 5 * N;
488: }
489: }
490: }
491: }
492: WORK[1 + o_work] = LWMIN;
493: IWORK[1 + o_iwork] = LIWMIN;
494: // *
495: if (LWORK < LWMIN && !LQUERY)
496: {
497: INFO = - 8;
498: }
499: else
500: {
501: if (LIWORK < LIWMIN && !LQUERY)
502: {
503: INFO = - 10;
504: }
505: }
506: }
507: // *
508: if (INFO != 0)
509: {
510: this._xerbla.Run("DSTEDC", - INFO);
511: return;
512: }
513: else
514: {
515: if (LQUERY)
516: {
517: return;
518: }
519: }
520: // *
521: // * Quick return if possible
522: // *
523: if (N == 0) return;
524: if (N == 1)
525: {
526: if (ICOMPZ != 0) Z[1+1 * LDZ + o_z] = ONE;
527: return;
528: }
529: // *
530: // * If the following conditional clause is removed, then the routine
531: // * will use the Divide and Conquer routine to compute only the
532: // * eigenvalues, which requires (3N + 3N**2) real workspace and
533: // * (2 + 5N + 2N lg(N)) integer workspace.
534: // * Since on many architectures DSTERF is much faster than any other
535: // * algorithm for finding eigenvalues only, it is used here
536: // * as the default. If the conditional clause is removed, then
537: // * information on the size of workspace needs to be changed.
538: // *
539: // * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
540: // *
541: if (ICOMPZ == 0)
542: {
543: this._dsterf.Run(N, ref D, offset_d, ref E, offset_e, ref INFO);
544: goto LABEL50;
545: }
546: // *
547: // * If N is smaller than the minimum divide size (SMLSIZ+1), then
548: // * solve the problem with another solver.
549: // *
550: if (N <= SMLSIZ)
551: {
552: // *
553: this._dsteqr.Run(COMPZ, N, ref D, offset_d, ref E, offset_e, ref Z, offset_z, LDZ
554: , ref WORK, offset_work, ref INFO);
555: // *
556: }
557: else
558: {
559: // *
560: // * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
561: // * use.
562: // *
563: if (ICOMPZ == 1)
564: {
565: STOREZ = 1 + N * N;
566: }
567: else
568: {
569: STOREZ = 1;
570: }
571: // *
572: if (ICOMPZ == 2)
573: {
574: this._dlaset.Run("Full", N, N, ZERO, ONE, ref Z, offset_z
575: , LDZ);
576: }
577: // *
578: // * Scale.
579: // *
580: ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);
581: if (ORGNRM == ZERO) goto LABEL50;
582: // *
583: EPS = this._dlamch.Run("Epsilon");
584: // *
585: START = 1;
586: // *
587: // * while ( START <= N )
588: // *
589: LABEL10:;
590: if (START <= N)
591: {
592: // *
593: // * Let FINISH be the position of the next subdiagonal entry
594: // * such that E( FINISH ) <= TINY or FINISH = N if no such
595: // * subdiagonal exists. The matrix identified by the elements
596: // * between START and FINISH constitutes an independent
597: // * sub-problem.
598: // *
599: FINISH = START;
600: LABEL20:;
601: if (FINISH < N)
602: {
603: TINY = EPS * Math.Sqrt(Math.Abs(D[FINISH + o_d])) * Math.Sqrt(Math.Abs(D[FINISH + 1 + o_d]));
604: if (Math.Abs(E[FINISH + o_e]) > TINY)
605: {
606: FINISH = FINISH + 1;
607: goto LABEL20;
608: }
609: }
610: // *
611: // * (Sub) Problem determined. Compute its size and solve it.
612: // *
613: M = FINISH - START + 1;
614: if (M == 1)
615: {
616: START = FINISH + 1;
617: goto LABEL10;
618: }
619: if (M > SMLSIZ)
620: {
621: // *
622: // * Scale.
623: // *
624: ORGNRM = this._dlanst.Run("M", M, D, START + o_d, E, START + o_e);
625: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, M
626: , 1, ref D, START + o_d, M, ref INFO);
627: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, M - 1
628: , 1, ref E, START + o_e, M - 1, ref INFO);
629: // *
630: if (ICOMPZ == 1)
631: {
632: STRTRW = 1;
633: }
634: else
635: {
636: STRTRW = START;
637: }
638: this._dlaed0.Run(ICOMPZ, N, M, ref D, START + o_d, ref E, START + o_e, ref Z, STRTRW+START * LDZ + o_z
639: , LDZ, ref WORK, 1 + o_work, N, ref WORK, STOREZ + o_work, ref IWORK, offset_iwork, ref INFO);
640: if (INFO != 0)
641: {
642: INFO = (INFO / (M + 1) + START - 1) * (N + 1) + FortranLib.Mod(INFO,(M + 1)) + START - 1;
643: goto LABEL50;
644: }
645: // *
646: // * Scale back.
647: // *
648: this._dlascl.Run("G", 0, 0, ONE, ORGNRM, M
649: , 1, ref D, START + o_d, M, ref INFO);
650: // *
651: }
652: else
653: {
654: if (ICOMPZ == 1)
655: {
656: // *
657: // * Since QR won't update a Z matrix which is larger than
658: // * the length of D, we must solve the sub-problem in a
659: // * workspace and then multiply back into Z.
660: // *
661: this._dsteqr.Run("I", M, ref D, START + o_d, ref E, START + o_e, ref WORK, offset_work, M
662: , ref WORK, M * M + 1 + o_work, ref INFO);
663: this._dlacpy.Run("A", N, M, Z, 1+START * LDZ + o_z, LDZ, ref WORK, STOREZ + o_work
664: , N);
665: this._dgemm.Run("N", "N", N, M, M, ONE
666: , WORK, STOREZ + o_work, N, WORK, offset_work, M, ZERO, ref Z, 1+START * LDZ + o_z
667: , LDZ);
668: }
669: else
670: {
671: if (ICOMPZ == 2)
672: {
673: this._dsteqr.Run("I", M, ref D, START + o_d, ref E, START + o_e, ref Z, START+START * LDZ + o_z, LDZ
674: , ref WORK, offset_work, ref INFO);
675: }
676: else
677: {
678: this._dsterf.Run(M, ref D, START + o_d, ref E, START + o_e, ref INFO);
679: }
680: }
681: if (INFO != 0)
682: {
683: INFO = START * (N + 1) + FINISH;
684: goto LABEL50;
685: }
686: }
687: // *
688: START = FINISH + 1;
689: goto LABEL10;
690: }
691: // *
692: // * endwhile
693: // *
694: // * If the problem split any number of times, then the eigenvalues
695: // * will not be properly ordered. Here we permute the eigenvalues
696: // * (and the associated eigenvectors) into ascending order.
697: // *
698: if (M != N)
699: {
700: if (ICOMPZ == 0)
701: {
702: // *
703: // * Use Quick Sort
704: // *
705: this._dlasrt.Run("I", N, ref D, offset_d, ref INFO);
706: // *
707: }
708: else
709: {
710: // *
711: // * Use Selection Sort to minimize swaps of eigenvectors
712: // *
713: for (II = 2; II <= N; II++)
714: {
715: I = II - 1;
716: K = I;
717: P = D[I + o_d];
718: for (J = II; J <= N; J++)
719: {
720: if (D[J + o_d] < P)
721: {
722: K = J;
723: P = D[J + o_d];
724: }
725: }
726: if (K != I)
727: {
728: D[K + o_d] = D[I + o_d];
729: D[I + o_d] = P;
730: this._dswap.Run(N, ref Z, 1+I * LDZ + o_z, 1, ref Z, 1+K * LDZ + o_z, 1);
731: }
732: }
733: }
734: }
735: }
736: // *
737: LABEL50:;
738: WORK[1 + o_work] = LWMIN;
739: IWORK[1 + o_iwork] = LIWMIN;
740: // *
741: return;
742: // *
743: // * End of DSTEDC
744: // *
745:
746: #endregion
747:
748: }
749: }
750: }