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 routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
27: /// symmetric tridiagonal matrix using the implicit QL or QR method.
28: /// The eigenvectors of a full or band symmetric matrix can also be found
29: /// if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
30: /// tridiagonal form.
31: ///
32: ///</summary>
33: public class DSTEQR
34: {
35:
36:
37: #region Dependencies
38:
39: LSAME _lsame; DLAMCH _dlamch; DLANST _dlanst; DLAPY2 _dlapy2; DLAE2 _dlae2; DLAEV2 _dlaev2; DLARTG _dlartg;
40: DLASCL _dlascl;DLASET _dlaset; DLASR _dlasr; DLASRT _dlasrt; DSWAP _dswap; XERBLA _xerbla;
41:
42: #endregion
43:
44:
45: #region Fields
46:
47: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double THREE = 3.0E0;
48: const int MAXIT = 30;int I = 0; int ICOMPZ = 0; int II = 0; int ISCALE = 0; int J = 0; int JTOT = 0; int K = 0; int L = 0;
49: int L1 = 0;int LEND = 0; int LENDM1 = 0; int LENDP1 = 0; int LENDSV = 0; int LM1 = 0; int LSV = 0; int M = 0; int MM = 0;
50: int MM1 = 0;int NM1 = 0; int NMAXIT = 0; double ANORM = 0; double B = 0; double C = 0; double EPS = 0; double EPS2 = 0;
51: double F = 0;double G = 0; double P = 0; double R = 0; double RT1 = 0; double RT2 = 0; double S = 0; double SAFMAX = 0;
52: double SAFMIN = 0;double SSFMAX = 0; double SSFMIN = 0; double TST = 0;
53:
54: #endregion
55:
56: public DSTEQR(LSAME lsame, DLAMCH dlamch, DLANST dlanst, DLAPY2 dlapy2, DLAE2 dlae2, DLAEV2 dlaev2, DLARTG dlartg, DLASCL dlascl, DLASET dlaset, DLASR dlasr
57: , DLASRT dlasrt, DSWAP dswap, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._lsame = lsame; this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2;
64: this._dlaev2 = dlaev2;this._dlartg = dlartg; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasr = dlasr;
65: this._dlasrt = dlasrt;this._dswap = dswap; this._xerbla = xerbla;
66:
67: #endregion
68:
69: }
70:
71: public DSTEQR()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: LSAME lsame = new LSAME();
78: DLAMC3 dlamc3 = new DLAMC3();
79: DLASSQ dlassq = new DLASSQ();
80: DLAPY2 dlapy2 = new DLAPY2();
81: DLAE2 dlae2 = new DLAE2();
82: DLAEV2 dlaev2 = new DLAEV2();
83: XERBLA xerbla = new XERBLA();
84: DSWAP dswap = new DSWAP();
85: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
86: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
87: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
88: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
89: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
90: DLANST dlanst = new DLANST(lsame, dlassq);
91: DLARTG dlartg = new DLARTG(dlamch);
92: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
93: DLASET dlaset = new DLASET(lsame);
94: DLASR dlasr = new DLASR(lsame, xerbla);
95: DLASRT dlasrt = new DLASRT(lsame, xerbla);
96:
97: #endregion
98:
99:
100: #region Set Dependencies
101:
102: this._lsame = lsame; this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2;
103: this._dlaev2 = dlaev2;this._dlartg = dlartg; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasr = dlasr;
104: this._dlasrt = dlasrt;this._dswap = dswap; this._xerbla = xerbla;
105:
106: #endregion
107:
108: }
109: /// <summary>
110: /// Purpose
111: /// =======
112: ///
113: /// DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
114: /// symmetric tridiagonal matrix using the implicit QL or QR method.
115: /// The eigenvectors of a full or band symmetric matrix can also be found
116: /// if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
117: /// tridiagonal form.
118: ///
119: ///</summary>
120: /// <param name="COMPZ">
121: /// (input) CHARACTER*1
122: /// = 'N': Compute eigenvalues only.
123: /// = 'V': Compute eigenvalues and eigenvectors of the original
124: /// symmetric matrix. On entry, Z must contain the
125: /// orthogonal matrix used to reduce the original matrix
126: /// to tridiagonal form.
127: /// = 'I': Compute eigenvalues and eigenvectors of the
128: /// tridiagonal matrix. Z is initialized to the identity
129: /// matrix.
130: ///</param>
131: /// <param name="N">
132: /// (input) INTEGER
133: /// The order of the matrix. N .GE. 0.
134: ///</param>
135: /// <param name="D">
136: /// (input/output) DOUBLE PRECISION array, dimension (N)
137: /// On entry, the diagonal elements of the tridiagonal matrix.
138: /// On exit, if INFO = 0, the eigenvalues in ascending order.
139: ///</param>
140: /// <param name="E">
141: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
142: /// On entry, the (n-1) subdiagonal elements of the tridiagonal
143: /// matrix.
144: /// On exit, E has been destroyed.
145: ///</param>
146: /// <param name="Z">
147: /// (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
148: /// On entry, if COMPZ = 'V', then Z contains the orthogonal
149: /// matrix used in the reduction to tridiagonal form.
150: /// On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
151: /// orthonormal eigenvectors of the original symmetric matrix,
152: /// and if COMPZ = 'I', Z contains the orthonormal eigenvectors
153: /// of the symmetric tridiagonal matrix.
154: /// If COMPZ = 'N', then Z is not referenced.
155: ///</param>
156: /// <param name="LDZ">
157: /// (input) INTEGER
158: /// The leading dimension of the array Z. LDZ .GE. 1, and if
159: /// eigenvectors are desired, then LDZ .GE. max(1,N).
160: ///</param>
161: /// <param name="WORK">
162: /// (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
163: /// If COMPZ = 'N', then WORK is not referenced.
164: ///</param>
165: /// <param name="INFO">
166: /// (output) INTEGER
167: /// = 0: successful exit
168: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
169: /// .GT. 0: the algorithm has failed to find all the eigenvalues in
170: /// a total of 30*N iterations; if INFO = i, then i
171: /// elements of E have not converged to zero; on exit, D
172: /// and E contain the elements of a symmetric tridiagonal
173: /// matrix which is orthogonally similar to the original
174: /// matrix.
175: ///</param>
176: 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
177: , ref double[] WORK, int offset_work, ref int INFO)
178: {
179:
180: #region Array Index Correction
181:
182: 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;
183:
184: #endregion
185:
186:
187: #region Strings
188:
189: COMPZ = COMPZ.Substring(0, 1);
190:
191: #endregion
192:
193:
194: #region Prolog
195:
196: // *
197: // * -- LAPACK routine (version 3.1) --
198: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
199: // * November 2006
200: // *
201: // * .. Scalar Arguments ..
202: // * ..
203: // * .. Array Arguments ..
204: // * ..
205: // *
206: // * Purpose
207: // * =======
208: // *
209: // * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
210: // * symmetric tridiagonal matrix using the implicit QL or QR method.
211: // * The eigenvectors of a full or band symmetric matrix can also be found
212: // * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
213: // * tridiagonal form.
214: // *
215: // * Arguments
216: // * =========
217: // *
218: // * COMPZ (input) CHARACTER*1
219: // * = 'N': Compute eigenvalues only.
220: // * = 'V': Compute eigenvalues and eigenvectors of the original
221: // * symmetric matrix. On entry, Z must contain the
222: // * orthogonal matrix used to reduce the original matrix
223: // * to tridiagonal form.
224: // * = 'I': Compute eigenvalues and eigenvectors of the
225: // * tridiagonal matrix. Z is initialized to the identity
226: // * matrix.
227: // *
228: // * N (input) INTEGER
229: // * The order of the matrix. N >= 0.
230: // *
231: // * D (input/output) DOUBLE PRECISION array, dimension (N)
232: // * On entry, the diagonal elements of the tridiagonal matrix.
233: // * On exit, if INFO = 0, the eigenvalues in ascending order.
234: // *
235: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
236: // * On entry, the (n-1) subdiagonal elements of the tridiagonal
237: // * matrix.
238: // * On exit, E has been destroyed.
239: // *
240: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
241: // * On entry, if COMPZ = 'V', then Z contains the orthogonal
242: // * matrix used in the reduction to tridiagonal form.
243: // * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
244: // * orthonormal eigenvectors of the original symmetric matrix,
245: // * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
246: // * of the symmetric tridiagonal matrix.
247: // * If COMPZ = 'N', then Z is not referenced.
248: // *
249: // * LDZ (input) INTEGER
250: // * The leading dimension of the array Z. LDZ >= 1, and if
251: // * eigenvectors are desired, then LDZ >= max(1,N).
252: // *
253: // * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
254: // * If COMPZ = 'N', then WORK is not referenced.
255: // *
256: // * INFO (output) INTEGER
257: // * = 0: successful exit
258: // * < 0: if INFO = -i, the i-th argument had an illegal value
259: // * > 0: the algorithm has failed to find all the eigenvalues in
260: // * a total of 30*N iterations; if INFO = i, then i
261: // * elements of E have not converged to zero; on exit, D
262: // * and E contain the elements of a symmetric tridiagonal
263: // * matrix which is orthogonally similar to the original
264: // * matrix.
265: // *
266: // * =====================================================================
267: // *
268: // * .. Parameters ..
269: // * ..
270: // * .. Local Scalars ..
271: // * ..
272: // * .. External Functions ..
273: // * ..
274: // * .. External Subroutines ..
275: // * ..
276: // * .. Intrinsic Functions ..
277: // INTRINSIC ABS, MAX, SIGN, SQRT;
278: // * ..
279: // * .. Executable Statements ..
280: // *
281: // * Test the input parameters.
282: // *
283:
284: #endregion
285:
286:
287: #region Body
288:
289: INFO = 0;
290: // *
291: if (this._lsame.Run(COMPZ, "N"))
292: {
293: ICOMPZ = 0;
294: }
295: else
296: {
297: if (this._lsame.Run(COMPZ, "V"))
298: {
299: ICOMPZ = 1;
300: }
301: else
302: {
303: if (this._lsame.Run(COMPZ, "I"))
304: {
305: ICOMPZ = 2;
306: }
307: else
308: {
309: ICOMPZ = - 1;
310: }
311: }
312: }
313: if (ICOMPZ < 0)
314: {
315: INFO = - 1;
316: }
317: else
318: {
319: if (N < 0)
320: {
321: INFO = - 2;
322: }
323: else
324: {
325: if ((LDZ < 1) || (ICOMPZ > 0 && LDZ < Math.Max(1, N)))
326: {
327: INFO = - 6;
328: }
329: }
330: }
331: if (INFO != 0)
332: {
333: this._xerbla.Run("DSTEQR", - INFO);
334: return;
335: }
336: // *
337: // * Quick return if possible
338: // *
339: if (N == 0) return;
340: // *
341: if (N == 1)
342: {
343: if (ICOMPZ == 2) Z[1+1 * LDZ + o_z] = ONE;
344: return;
345: }
346: // *
347: // * Determine the unit roundoff and over/underflow thresholds.
348: // *
349: EPS = this._dlamch.Run("E");
350: EPS2 = Math.Pow(EPS,2);
351: SAFMIN = this._dlamch.Run("S");
352: SAFMAX = ONE / SAFMIN;
353: SSFMAX = Math.Sqrt(SAFMAX) / THREE;
354: SSFMIN = Math.Sqrt(SAFMIN) / EPS2;
355: // *
356: // * Compute the eigenvalues and eigenvectors of the tridiagonal
357: // * matrix.
358: // *
359: if (ICOMPZ == 2)
360: {
361: this._dlaset.Run("Full", N, N, ZERO, ONE, ref Z, offset_z
362: , LDZ);
363: }
364: // *
365: NMAXIT = N * MAXIT;
366: JTOT = 0;
367: // *
368: // * Determine where the matrix splits and choose QL or QR iteration
369: // * for each block, according to whether top or bottom diagonal
370: // * element is smaller.
371: // *
372: L1 = 1;
373: NM1 = N - 1;
374: // *
375: LABEL10:;
376: if (L1 > N) goto LABEL160;
377: if (L1 > 1) E[L1 - 1 + o_e] = ZERO;
378: if (L1 <= NM1)
379: {
380: for (M = L1; M <= NM1; M++)
381: {
382: TST = Math.Abs(E[M + o_e]);
383: if (TST == ZERO) goto LABEL30;
384: if (TST <= (Math.Sqrt(Math.Abs(D[M + o_d])) * Math.Sqrt(Math.Abs(D[M + 1 + o_d]))) * EPS)
385: {
386: E[M + o_e] = ZERO;
387: goto LABEL30;
388: }
389: }
390: }
391: M = N;
392: // *
393: LABEL30:;
394: L = L1;
395: LSV = L;
396: LEND = M;
397: LENDSV = LEND;
398: L1 = M + 1;
399: if (LEND == L) goto LABEL10;
400: // *
401: // * Scale submatrix in rows and columns L to LEND
402: // *
403: ANORM = this._dlanst.Run("I", LEND - L + 1, D, L + o_d, E, L + o_e);
404: ISCALE = 0;
405: if (ANORM == ZERO) goto LABEL10;
406: if (ANORM > SSFMAX)
407: {
408: ISCALE = 1;
409: this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L + 1
410: , 1, ref D, L + o_d, N, ref INFO);
411: this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L
412: , 1, ref E, L + o_e, N, ref INFO);
413: }
414: else
415: {
416: if (ANORM < SSFMIN)
417: {
418: ISCALE = 2;
419: this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L + 1
420: , 1, ref D, L + o_d, N, ref INFO);
421: this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L
422: , 1, ref E, L + o_e, N, ref INFO);
423: }
424: }
425: // *
426: // * Choose between QL and QR iteration
427: // *
428: if (Math.Abs(D[LEND + o_d]) < Math.Abs(D[L + o_d]))
429: {
430: LEND = LSV;
431: L = LENDSV;
432: }
433: // *
434: if (LEND > L)
435: {
436: // *
437: // * QL Iteration
438: // *
439: // * Look for small subdiagonal element.
440: // *
441: LABEL40:;
442: if (L != LEND)
443: {
444: LENDM1 = LEND - 1;
445: for (M = L; M <= LENDM1; M++)
446: {
447: TST = Math.Pow(Math.Abs(E[M + o_e]),2);
448: if (TST <= (EPS2 * Math.Abs(D[M + o_d])) * Math.Abs(D[M + 1 + o_d]) + SAFMIN) goto LABEL60;
449: }
450: }
451: // *
452: M = LEND;
453: // *
454: LABEL60:;
455: if (M < LEND) E[M + o_e] = ZERO;
456: P = D[L + o_d];
457: if (M == L) goto LABEL80;
458: // *
459: // * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
460: // * to compute its eigensystem.
461: // *
462: if (M == L + 1)
463: {
464: if (ICOMPZ > 0)
465: {
466: this._dlaev2.Run(D[L + o_d], E[L + o_e], D[L + 1 + o_d], ref RT1, ref RT2, ref C
467: , ref S);
468: WORK[L + o_work] = C;
469: WORK[N - 1 + L + o_work] = S;
470: this._dlasr.Run("R", "V", "B", N, 2, WORK, L + o_work
471: , WORK, N - 1 + L + o_work, ref Z, 1+L * LDZ + o_z, LDZ);
472: }
473: else
474: {
475: this._dlae2.Run(D[L + o_d], E[L + o_e], D[L + 1 + o_d], ref RT1, ref RT2);
476: }
477: D[L + o_d] = RT1;
478: D[L + 1 + o_d] = RT2;
479: E[L + o_e] = ZERO;
480: L = L + 2;
481: if (L <= LEND) goto LABEL40;
482: goto LABEL140;
483: }
484: // *
485: if (JTOT == NMAXIT) goto LABEL140;
486: JTOT = JTOT + 1;
487: // *
488: // * Form shift.
489: // *
490: G = (D[L + 1 + o_d] - P) / (TWO * E[L + o_e]);
491: R = this._dlapy2.Run(G, ONE);
492: G = D[M + o_d] - P + (E[L + o_e] / (G + FortranLib.Sign(R,G)));
493: // *
494: S = ONE;
495: C = ONE;
496: P = ZERO;
497: // *
498: // * Inner loop
499: // *
500: MM1 = M - 1;
501: for (I = MM1; I >= L; I += - 1)
502: {
503: F = S * E[I + o_e];
504: B = C * E[I + o_e];
505: this._dlartg.Run(G, F, ref C, ref S, ref R);
506: if (I != M - 1) E[I + 1 + o_e] = R;
507: G = D[I + 1 + o_d] - P;
508: R = (D[I + o_d] - G) * S + TWO * C * B;
509: P = S * R;
510: D[I + 1 + o_d] = G + P;
511: G = C * R - B;
512: // *
513: // * If eigenvectors are desired, then save rotations.
514: // *
515: if (ICOMPZ > 0)
516: {
517: WORK[I + o_work] = C;
518: WORK[N - 1 + I + o_work] = - S;
519: }
520: // *
521: }
522: // *
523: // * If eigenvectors are desired, then apply saved rotations.
524: // *
525: if (ICOMPZ > 0)
526: {
527: MM = M - L + 1;
528: this._dlasr.Run("R", "V", "B", N, MM, WORK, L + o_work
529: , WORK, N - 1 + L + o_work, ref Z, 1+L * LDZ + o_z, LDZ);
530: }
531: // *
532: D[L + o_d] = D[L + o_d] - P;
533: E[L + o_e] = G;
534: goto LABEL40;
535: // *
536: // * Eigenvalue found.
537: // *
538: LABEL80:;
539: D[L + o_d] = P;
540: // *
541: L = L + 1;
542: if (L <= LEND) goto LABEL40;
543: goto LABEL140;
544: // *
545: }
546: else
547: {
548: // *
549: // * QR Iteration
550: // *
551: // * Look for small superdiagonal element.
552: // *
553: LABEL90:;
554: if (L != LEND)
555: {
556: LENDP1 = LEND + 1;
557: for (M = L; M >= LENDP1; M += - 1)
558: {
559: TST = Math.Pow(Math.Abs(E[M - 1 + o_e]),2);
560: if (TST <= (EPS2 * Math.Abs(D[M + o_d])) * Math.Abs(D[M - 1 + o_d]) + SAFMIN) goto LABEL110;
561: }
562: }
563: // *
564: M = LEND;
565: // *
566: LABEL110:;
567: if (M > LEND) E[M - 1 + o_e] = ZERO;
568: P = D[L + o_d];
569: if (M == L) goto LABEL130;
570: // *
571: // * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
572: // * to compute its eigensystem.
573: // *
574: if (M == L - 1)
575: {
576: if (ICOMPZ > 0)
577: {
578: this._dlaev2.Run(D[L - 1 + o_d], E[L - 1 + o_e], D[L + o_d], ref RT1, ref RT2, ref C
579: , ref S);
580: WORK[M + o_work] = C;
581: WORK[N - 1 + M + o_work] = S;
582: this._dlasr.Run("R", "V", "F", N, 2, WORK, M + o_work
583: , WORK, N - 1 + M + o_work, ref Z, 1+(L - 1) * LDZ + o_z, LDZ);
584: }
585: else
586: {
587: this._dlae2.Run(D[L - 1 + o_d], E[L - 1 + o_e], D[L + o_d], ref RT1, ref RT2);
588: }
589: D[L - 1 + o_d] = RT1;
590: D[L + o_d] = RT2;
591: E[L - 1 + o_e] = ZERO;
592: L = L - 2;
593: if (L >= LEND) goto LABEL90;
594: goto LABEL140;
595: }
596: // *
597: if (JTOT == NMAXIT) goto LABEL140;
598: JTOT = JTOT + 1;
599: // *
600: // * Form shift.
601: // *
602: G = (D[L - 1 + o_d] - P) / (TWO * E[L - 1 + o_e]);
603: R = this._dlapy2.Run(G, ONE);
604: G = D[M + o_d] - P + (E[L - 1 + o_e] / (G + FortranLib.Sign(R,G)));
605: // *
606: S = ONE;
607: C = ONE;
608: P = ZERO;
609: // *
610: // * Inner loop
611: // *
612: LM1 = L - 1;
613: for (I = M; I <= LM1; I++)
614: {
615: F = S * E[I + o_e];
616: B = C * E[I + o_e];
617: this._dlartg.Run(G, F, ref C, ref S, ref R);
618: if (I != M) E[I - 1 + o_e] = R;
619: G = D[I + o_d] - P;
620: R = (D[I + 1 + o_d] - G) * S + TWO * C * B;
621: P = S * R;
622: D[I + o_d] = G + P;
623: G = C * R - B;
624: // *
625: // * If eigenvectors are desired, then save rotations.
626: // *
627: if (ICOMPZ > 0)
628: {
629: WORK[I + o_work] = C;
630: WORK[N - 1 + I + o_work] = S;
631: }
632: // *
633: }
634: // *
635: // * If eigenvectors are desired, then apply saved rotations.
636: // *
637: if (ICOMPZ > 0)
638: {
639: MM = L - M + 1;
640: this._dlasr.Run("R", "V", "F", N, MM, WORK, M + o_work
641: , WORK, N - 1 + M + o_work, ref Z, 1+M * LDZ + o_z, LDZ);
642: }
643: // *
644: D[L + o_d] = D[L + o_d] - P;
645: E[LM1 + o_e] = G;
646: goto LABEL90;
647: // *
648: // * Eigenvalue found.
649: // *
650: LABEL130:;
651: D[L + o_d] = P;
652: // *
653: L = L - 1;
654: if (L >= LEND) goto LABEL90;
655: goto LABEL140;
656: // *
657: }
658: // *
659: // * Undo scaling if necessary
660: // *
661: LABEL140:;
662: if (ISCALE == 1)
663: {
664: this._dlascl.Run("G", 0, 0, SSFMAX, ANORM, LENDSV - LSV + 1
665: , 1, ref D, LSV + o_d, N, ref INFO);
666: this._dlascl.Run("G", 0, 0, SSFMAX, ANORM, LENDSV - LSV
667: , 1, ref E, LSV + o_e, N, ref INFO);
668: }
669: else
670: {
671: if (ISCALE == 2)
672: {
673: this._dlascl.Run("G", 0, 0, SSFMIN, ANORM, LENDSV - LSV + 1
674: , 1, ref D, LSV + o_d, N, ref INFO);
675: this._dlascl.Run("G", 0, 0, SSFMIN, ANORM, LENDSV - LSV
676: , 1, ref E, LSV + o_e, N, ref INFO);
677: }
678: }
679: // *
680: // * Check for no convergence to an eigenvalue after a total
681: // * of N*MAXIT iterations.
682: // *
683: if (JTOT < NMAXIT) goto LABEL10;
684: for (I = 1; I <= N - 1; I++)
685: {
686: if (E[I + o_e] != ZERO) INFO = INFO + 1;
687: }
688: goto LABEL190;
689: // *
690: // * Order eigenvalues and eigenvectors.
691: // *
692: LABEL160:;
693: if (ICOMPZ == 0)
694: {
695: // *
696: // * Use Quick Sort
697: // *
698: this._dlasrt.Run("I", N, ref D, offset_d, ref INFO);
699: // *
700: }
701: else
702: {
703: // *
704: // * Use Selection Sort to minimize swaps of eigenvectors
705: // *
706: for (II = 2; II <= N; II++)
707: {
708: I = II - 1;
709: K = I;
710: P = D[I + o_d];
711: for (J = II; J <= N; J++)
712: {
713: if (D[J + o_d] < P)
714: {
715: K = J;
716: P = D[J + o_d];
717: }
718: }
719: if (K != I)
720: {
721: D[K + o_d] = D[I + o_d];
722: D[I + o_d] = P;
723: this._dswap.Run(N, ref Z, 1+I * LDZ + o_z, 1, ref Z, 1+K * LDZ + o_z, 1);
724: }
725: }
726: }
727: // *
728: LABEL190:;
729: return;
730: // *
731: // * End of DSTEQR
732: // *
733:
734: #endregion
735:
736: }
737: }
738: }