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: /// DLASQ2 computes all the eigenvalues of the symmetric positive
27: /// definite tridiagonal matrix associated with the qd array Z to high
28: /// relative accuracy are computed to high relative accuracy, in the
29: /// absence of denormalization, underflow and overflow.
30: ///
31: /// To see the relation of Z to the tridiagonal matrix, let L be a
32: /// unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
33: /// let U be an upper bidiagonal matrix with 1's above and diagonal
34: /// Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
35: /// symmetric tridiagonal to which it is similar.
36: ///
37: /// Note : DLASQ2 defines a logical variable, IEEE, which is true
38: /// on machines which follow ieee-754 floating-point standard in their
39: /// handling of infinities and NaNs, and false otherwise. This variable
40: /// is passed to DLAZQ3.
41: ///
42: ///</summary>
43: public class DLASQ2
44: {
45:
46:
47: #region Dependencies
48:
49: DLAZQ3 _dlazq3; DLASRT _dlasrt; XERBLA _xerbla; DLAMCH _dlamch; ILAENV _ilaenv;
50:
51: #endregion
52:
53:
54: #region Fields
55:
56: const double CBIAS = 1.50E0; const double ZERO = 0.0E0; const double HALF = 0.5E0; const double ONE = 1.0E0;
57: const double TWO = 2.0E0;const double FOUR = 4.0E0; const double HUNDRD = 100.0E0; bool IEEE = false; int I0 = 0;
58: int I4 = 0;int IINFO = 0; int IPN4 = 0; int ITER = 0; int IWHILA = 0; int IWHILB = 0; int K = 0; int N0 = 0; int NBIG = 0;
59: int NDIV = 0;int NFAIL = 0; int PP = 0; int SPLT = 0; int TTYPE = 0; double D = 0; double DESIG = 0; double DMIN = 0;
60: double DMIN1 = 0;double DMIN2 = 0; double DN = 0; double DN1 = 0; double DN2 = 0; double E = 0; double EMAX = 0;
61: double EMIN = 0;double EPS = 0; double OLDEMN = 0; double QMAX = 0; double QMIN = 0; double S = 0; double SAFMIN = 0;
62: double SIGMA = 0;double T = 0; double TAU = 0; double TEMP = 0; double TOL = 0; double TOL2 = 0; double TRACE = 0;
63: double ZMAX = 0;
64:
65: #endregion
66:
67: public DLASQ2(DLAZQ3 dlazq3, DLASRT dlasrt, XERBLA xerbla, DLAMCH dlamch, ILAENV ilaenv)
68: {
69:
70:
71: #region Set Dependencies
72:
73: this._dlazq3 = dlazq3; this._dlasrt = dlasrt; this._xerbla = xerbla; this._dlamch = dlamch; this._ilaenv = ilaenv;
74:
75: #endregion
76:
77: }
78:
79: public DLASQ2()
80: {
81:
82:
83: #region Dependencies (Initialization)
84:
85: DLASQ5 dlasq5 = new DLASQ5();
86: LSAME lsame = new LSAME();
87: DLAMC3 dlamc3 = new DLAMC3();
88: DLAZQ4 dlazq4 = new DLAZQ4();
89: XERBLA xerbla = new XERBLA();
90: IEEECK ieeeck = new IEEECK();
91: IPARMQ iparmq = new IPARMQ();
92: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
93: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
94: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
95: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
96: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
97: DLASQ6 dlasq6 = new DLASQ6(dlamch);
98: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
99: DLASRT dlasrt = new DLASRT(lsame, xerbla);
100: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
101:
102: #endregion
103:
104:
105: #region Set Dependencies
106:
107: this._dlazq3 = dlazq3; this._dlasrt = dlasrt; this._xerbla = xerbla; this._dlamch = dlamch; this._ilaenv = ilaenv;
108:
109: #endregion
110:
111: }
112: /// <summary>
113: /// Purpose
114: /// =======
115: ///
116: /// DLASQ2 computes all the eigenvalues of the symmetric positive
117: /// definite tridiagonal matrix associated with the qd array Z to high
118: /// relative accuracy are computed to high relative accuracy, in the
119: /// absence of denormalization, underflow and overflow.
120: ///
121: /// To see the relation of Z to the tridiagonal matrix, let L be a
122: /// unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
123: /// let U be an upper bidiagonal matrix with 1's above and diagonal
124: /// Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
125: /// symmetric tridiagonal to which it is similar.
126: ///
127: /// Note : DLASQ2 defines a logical variable, IEEE, which is true
128: /// on machines which follow ieee-754 floating-point standard in their
129: /// handling of infinities and NaNs, and false otherwise. This variable
130: /// is passed to DLAZQ3.
131: ///
132: ///</summary>
133: /// <param name="N">
134: /// (input) INTEGER
135: /// The number of rows and columns in the matrix. N .GE. 0.
136: ///</param>
137: /// <param name="Z">
138: /// (workspace) DOUBLE PRECISION array, dimension ( 4*N )
139: /// On entry Z holds the qd array. On exit, entries 1 to N hold
140: /// the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
141: /// trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
142: /// N .GT. 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
143: /// holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
144: /// shifts that failed.
145: ///</param>
146: /// <param name="INFO">
147: /// (output) INTEGER
148: /// = 0: successful exit
149: /// .LT. 0: if the i-th argument is a scalar and had an illegal
150: /// value, then INFO = -i, if the i-th argument is an
151: /// array and the j-entry had an illegal value, then
152: /// INFO = -(i*100+j)
153: /// .GT. 0: the algorithm failed
154: /// = 1, a split was marked by a positive value in E
155: /// = 2, current block of Z not diagonalized after 30*N
156: /// iterations (in inner while loop)
157: /// = 3, termination criterion of outer while loop not met
158: /// (program created more than N unreduced blocks)
159: ///</param>
160: public void Run(int N, ref double[] Z, int offset_z, ref int INFO)
161: {
162:
163: #region Array Index Correction
164:
165: int o_z = -1 + offset_z;
166:
167: #endregion
168:
169:
170: #region Prolog
171:
172: // *
173: // * -- LAPACK routine (version 3.1) --
174: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
175: // * November 2006
176: // *
177: // * Modified to call DLAZQ3 in place of DLASQ3, 13 Feb 03, SJH.
178: // *
179: // * .. Scalar Arguments ..
180: // * ..
181: // * .. Array Arguments ..
182: // * ..
183: // *
184: // * Purpose
185: // * =======
186: // *
187: // * DLASQ2 computes all the eigenvalues of the symmetric positive
188: // * definite tridiagonal matrix associated with the qd array Z to high
189: // * relative accuracy are computed to high relative accuracy, in the
190: // * absence of denormalization, underflow and overflow.
191: // *
192: // * To see the relation of Z to the tridiagonal matrix, let L be a
193: // * unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
194: // * let U be an upper bidiagonal matrix with 1's above and diagonal
195: // * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
196: // * symmetric tridiagonal to which it is similar.
197: // *
198: // * Note : DLASQ2 defines a logical variable, IEEE, which is true
199: // * on machines which follow ieee-754 floating-point standard in their
200: // * handling of infinities and NaNs, and false otherwise. This variable
201: // * is passed to DLAZQ3.
202: // *
203: // * Arguments
204: // * =========
205: // *
206: // * N (input) INTEGER
207: // * The number of rows and columns in the matrix. N >= 0.
208: // *
209: // * Z (workspace) DOUBLE PRECISION array, dimension ( 4*N )
210: // * On entry Z holds the qd array. On exit, entries 1 to N hold
211: // * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
212: // * trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
213: // * N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
214: // * holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
215: // * shifts that failed.
216: // *
217: // * INFO (output) INTEGER
218: // * = 0: successful exit
219: // * < 0: if the i-th argument is a scalar and had an illegal
220: // * value, then INFO = -i, if the i-th argument is an
221: // * array and the j-entry had an illegal value, then
222: // * INFO = -(i*100+j)
223: // * > 0: the algorithm failed
224: // * = 1, a split was marked by a positive value in E
225: // * = 2, current block of Z not diagonalized after 30*N
226: // * iterations (in inner while loop)
227: // * = 3, termination criterion of outer while loop not met
228: // * (program created more than N unreduced blocks)
229: // *
230: // * Further Details
231: // * ===============
232: // * Local Variables: I0:N0 defines a current unreduced segment of Z.
233: // * The shifts are accumulated in SIGMA. Iteration count is in ITER.
234: // * Ping-pong is controlled by PP (alternates between 0 and 1).
235: // *
236: // * =====================================================================
237: // *
238: // * .. Parameters ..
239: // * ..
240: // * .. Local Scalars ..
241: // * ..
242: // * .. External Subroutines ..
243: // * ..
244: // * .. External Functions ..
245: // * ..
246: // * .. Intrinsic Functions ..
247: // INTRINSIC ABS, DBLE, MAX, MIN, SQRT;
248: // * ..
249: // * .. Executable Statements ..
250: // *
251: // * Test the input arguments.
252: // * (in case DLASQ2 is not called by DLASQ1)
253: // *
254:
255: #endregion
256:
257:
258: #region Body
259:
260: INFO = 0;
261: EPS = this._dlamch.Run("Precision");
262: SAFMIN = this._dlamch.Run("Safe minimum");
263: TOL = EPS * HUNDRD;
264: TOL2 = Math.Pow(TOL,2);
265: // *
266: if (N < 0)
267: {
268: INFO = - 1;
269: this._xerbla.Run("DLASQ2", 1);
270: return;
271: }
272: else
273: {
274: if (N == 0)
275: {
276: return;
277: }
278: else
279: {
280: if (N == 1)
281: {
282: // *
283: // * 1-by-1 case.
284: // *
285: if (Z[1 + o_z] < ZERO)
286: {
287: INFO = - 201;
288: this._xerbla.Run("DLASQ2", 2);
289: }
290: return;
291: }
292: else
293: {
294: if (N == 2)
295: {
296: // *
297: // * 2-by-2 case.
298: // *
299: if (Z[2 + o_z] < ZERO || Z[3 + o_z] < ZERO)
300: {
301: INFO = - 2;
302: this._xerbla.Run("DLASQ2", 2);
303: return;
304: }
305: else
306: {
307: if (Z[3 + o_z] > Z[1 + o_z])
308: {
309: D = Z[3 + o_z];
310: Z[3 + o_z] = Z[1 + o_z];
311: Z[1 + o_z] = D;
312: }
313: }
314: Z[5 + o_z] = Z[1 + o_z] + Z[2 + o_z] + Z[3 + o_z];
315: if (Z[2 + o_z] > Z[3 + o_z] * TOL2)
316: {
317: T = HALF * ((Z[1 + o_z] - Z[3 + o_z]) + Z[2 + o_z]);
318: S = Z[3 + o_z] * (Z[2 + o_z] / T);
319: if (S <= T)
320: {
321: S = Z[3 + o_z] * (Z[2 + o_z] / (T * (ONE + Math.Sqrt(ONE + S / T))));
322: }
323: else
324: {
325: S = Z[3 + o_z] * (Z[2 + o_z] / (T + Math.Sqrt(T) * Math.Sqrt(T + S)));
326: }
327: T = Z[1 + o_z] + (S + Z[2 + o_z]);
328: Z[3 + o_z] = Z[3 + o_z] * (Z[1 + o_z] / T);
329: Z[1 + o_z] = T;
330: }
331: Z[2 + o_z] = Z[3 + o_z];
332: Z[6 + o_z] = Z[2 + o_z] + Z[1 + o_z];
333: return;
334: }
335: }
336: }
337: }
338: // *
339: // * Check for negative data and compute sums of q's and e's.
340: // *
341: Z[2 * N + o_z] = ZERO;
342: EMIN = Z[2 + o_z];
343: QMAX = ZERO;
344: ZMAX = ZERO;
345: D = ZERO;
346: E = ZERO;
347: // *
348: for (K = 1; K <= 2 * (N - 1); K += 2)
349: {
350: if (Z[K + o_z] < ZERO)
351: {
352: INFO = - (200 + K);
353: this._xerbla.Run("DLASQ2", 2);
354: return;
355: }
356: else
357: {
358: if (Z[K + 1 + o_z] < ZERO)
359: {
360: INFO = - (200 + K + 1);
361: this._xerbla.Run("DLASQ2", 2);
362: return;
363: }
364: }
365: D = D + Z[K + o_z];
366: E = E + Z[K + 1 + o_z];
367: QMAX = Math.Max(QMAX, Z[K + o_z]);
368: EMIN = Math.Min(EMIN, Z[K + 1 + o_z]);
369: ZMAX = Math.Max(QMAX, Math.Max(ZMAX, Z[K + 1 + o_z]));
370: }
371: if (Z[2 * N - 1 + o_z] < ZERO)
372: {
373: INFO = - (200 + 2 * N - 1);
374: this._xerbla.Run("DLASQ2", 2);
375: return;
376: }
377: D = D + Z[2 * N - 1 + o_z];
378: QMAX = Math.Max(QMAX, Z[2 * N - 1 + o_z]);
379: ZMAX = Math.Max(QMAX, ZMAX);
380: // *
381: // * Check for diagonality.
382: // *
383: if (E == ZERO)
384: {
385: for (K = 2; K <= N; K++)
386: {
387: Z[K + o_z] = Z[2 * K - 1 + o_z];
388: }
389: this._dlasrt.Run("D", N, ref Z, offset_z, ref IINFO);
390: Z[2 * N - 1 + o_z] = D;
391: return;
392: }
393: // *
394: TRACE = D + E;
395: // *
396: // * Check for zero data.
397: // *
398: if (TRACE == ZERO)
399: {
400: Z[2 * N - 1 + o_z] = ZERO;
401: return;
402: }
403: // *
404: // * Check whether the machine is IEEE conformable.
405: // *
406: IEEE = this._ilaenv.Run(10, "DLASQ2", "N", 1, 2, 3, 4) == 1 && this._ilaenv.Run(11, "DLASQ2", "N", 1, 2, 3, 4) == 1;
407: // *
408: // * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
409: // *
410: for (K = 2 * N; K >= 2; K += - 2)
411: {
412: Z[2 * K + o_z] = ZERO;
413: Z[2 * K - 1 + o_z] = Z[K + o_z];
414: Z[2 * K - 2 + o_z] = ZERO;
415: Z[2 * K - 3 + o_z] = Z[K - 1 + o_z];
416: }
417: // *
418: I0 = 1;
419: N0 = N;
420: // *
421: // * Reverse the qd-array, if warranted.
422: // *
423: if (CBIAS * Z[4 * I0 - 3 + o_z] < Z[4 * N0 - 3 + o_z])
424: {
425: IPN4 = 4 * (I0 + N0);
426: for (I4 = 4 * I0; I4 <= 2 * (I0 + N0 - 1); I4 += 4)
427: {
428: TEMP = Z[I4 - 3 + o_z];
429: Z[I4 - 3 + o_z] = Z[IPN4 - I4 - 3 + o_z];
430: Z[IPN4 - I4 - 3 + o_z] = TEMP;
431: TEMP = Z[I4 - 1 + o_z];
432: Z[I4 - 1 + o_z] = Z[IPN4 - I4 - 5 + o_z];
433: Z[IPN4 - I4 - 5 + o_z] = TEMP;
434: }
435: }
436: // *
437: // * Initial split checking via dqd and Li's test.
438: // *
439: PP = 0;
440: // *
441: for (K = 1; K <= 2; K++)
442: {
443: // *
444: D = Z[4 * N0 + PP - 3 + o_z];
445: for (I4 = 4 * (N0 - 1) + PP; I4 >= 4 * I0 + PP; I4 += - 4)
446: {
447: if (Z[I4 - 1 + o_z] <= TOL2 * D)
448: {
449: Z[I4 - 1 + o_z] = - ZERO;
450: D = Z[I4 - 3 + o_z];
451: }
452: else
453: {
454: D = Z[I4 - 3 + o_z] * (D / (D + Z[I4 - 1 + o_z]));
455: }
456: }
457: // *
458: // * dqd maps Z to ZZ plus Li's test.
459: // *
460: EMIN = Z[4 * I0 + PP + 1 + o_z];
461: D = Z[4 * I0 + PP - 3 + o_z];
462: for (I4 = 4 * I0 + PP; I4 <= 4 * (N0 - 1) + PP; I4 += 4)
463: {
464: Z[I4 - 2 * PP - 2 + o_z] = D + Z[I4 - 1 + o_z];
465: if (Z[I4 - 1 + o_z] <= TOL2 * D)
466: {
467: Z[I4 - 1 + o_z] = - ZERO;
468: Z[I4 - 2 * PP - 2 + o_z] = D;
469: Z[I4 - 2 * PP + o_z] = ZERO;
470: D = Z[I4 + 1 + o_z];
471: }
472: else
473: {
474: if (SAFMIN * Z[I4 + 1 + o_z] < Z[I4 - 2 * PP - 2 + o_z] && SAFMIN * Z[I4 - 2 * PP - 2 + o_z] < Z[I4 + 1 + o_z])
475: {
476: TEMP = Z[I4 + 1 + o_z] / Z[I4 - 2 * PP - 2 + o_z];
477: Z[I4 - 2 * PP + o_z] = Z[I4 - 1 + o_z] * TEMP;
478: D = D * TEMP;
479: }
480: else
481: {
482: Z[I4 - 2 * PP + o_z] = Z[I4 + 1 + o_z] * (Z[I4 - 1 + o_z] / Z[I4 - 2 * PP - 2 + o_z]);
483: D = Z[I4 + 1 + o_z] * (D / Z[I4 - 2 * PP - 2 + o_z]);
484: }
485: }
486: EMIN = Math.Min(EMIN, Z[I4 - 2 * PP + o_z]);
487: }
488: Z[4 * N0 - PP - 2 + o_z] = D;
489: // *
490: // * Now find qmax.
491: // *
492: QMAX = Z[4 * I0 - PP - 2 + o_z];
493: for (I4 = 4 * I0 - PP + 2; I4 <= 4 * N0 - PP - 2; I4 += 4)
494: {
495: QMAX = Math.Max(QMAX, Z[I4 + o_z]);
496: }
497: // *
498: // * Prepare for the next iteration on K.
499: // *
500: PP = 1 - PP;
501: }
502: // *
503: // * Initialise variables to pass to DLAZQ3
504: // *
505: TTYPE = 0;
506: DMIN1 = ZERO;
507: DMIN2 = ZERO;
508: DN = ZERO;
509: DN1 = ZERO;
510: DN2 = ZERO;
511: TAU = ZERO;
512: // *
513: ITER = 2;
514: NFAIL = 0;
515: NDIV = 2 * (N0 - I0);
516: // *
517: for (IWHILA = 1; IWHILA <= N + 1; IWHILA++)
518: {
519: if (N0 < 1) goto LABEL150;
520: // *
521: // * While array unfinished do
522: // *
523: // * E(N0) holds the value of SIGMA when submatrix in I0:N0
524: // * splits from the rest of the array, but is negated.
525: // *
526: DESIG = ZERO;
527: if (N0 == N)
528: {
529: SIGMA = ZERO;
530: }
531: else
532: {
533: SIGMA = - Z[4 * N0 - 1 + o_z];
534: }
535: if (SIGMA < ZERO)
536: {
537: INFO = 1;
538: return;
539: }
540: // *
541: // * Find last unreduced submatrix's top index I0, find QMAX and
542: // * EMIN. Find Gershgorin-type bound if Q's much greater than E's.
543: // *
544: EMAX = ZERO;
545: if (N0 > I0)
546: {
547: EMIN = Math.Abs(Z[4 * N0 - 5 + o_z]);
548: }
549: else
550: {
551: EMIN = ZERO;
552: }
553: QMIN = Z[4 * N0 - 3 + o_z];
554: QMAX = QMIN;
555: for (I4 = 4 * N0; I4 >= 8; I4 += - 4)
556: {
557: if (Z[I4 - 5 + o_z] <= ZERO) goto LABEL100;
558: if (QMIN >= FOUR * EMAX)
559: {
560: QMIN = Math.Min(QMIN, Z[I4 - 3 + o_z]);
561: EMAX = Math.Max(EMAX, Z[I4 - 5 + o_z]);
562: }
563: QMAX = Math.Max(QMAX, Z[I4 - 7 + o_z] + Z[I4 - 5 + o_z]);
564: EMIN = Math.Min(EMIN, Z[I4 - 5 + o_z]);
565: }
566: I4 = 4;
567: // *
568: LABEL100:;
569: I0 = I4 / 4;
570: // *
571: // * Store EMIN for passing to DLAZQ3.
572: // *
573: Z[4 * N0 - 1 + o_z] = EMIN;
574: // *
575: // * Put -(initial shift) into DMIN.
576: // *
577: DMIN = - Math.Max(ZERO, QMIN - TWO * Math.Sqrt(QMIN) * Math.Sqrt(EMAX));
578: // *
579: // * Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.
580: // *
581: PP = 0;
582: // *
583: NBIG = 30 * (N0 - I0 + 1);
584: for (IWHILB = 1; IWHILB <= NBIG; IWHILB++)
585: {
586: if (I0 > N0) goto LABEL130;
587: // *
588: // * While submatrix unfinished take a good dqds step.
589: // *
590: this._dlazq3.Run(I0, ref N0, ref Z, offset_z, PP, ref DMIN, ref SIGMA
591: , ref DESIG, ref QMAX, ref NFAIL, ref ITER, ref NDIV, IEEE
592: , ref TTYPE, ref DMIN1, ref DMIN2, ref DN, ref DN1, ref DN2
593: , ref TAU);
594: // *
595: PP = 1 - PP;
596: // *
597: // * When EMIN is very small check for splits.
598: // *
599: if (PP == 0 && N0 - I0 >= 3)
600: {
601: if (Z[4 * N0 + o_z] <= TOL2 * QMAX || Z[4 * N0 - 1 + o_z] <= TOL2 * SIGMA)
602: {
603: SPLT = I0 - 1;
604: QMAX = Z[4 * I0 - 3 + o_z];
605: EMIN = Z[4 * I0 - 1 + o_z];
606: OLDEMN = Z[4 * I0 + o_z];
607: for (I4 = 4 * I0; I4 <= 4 * (N0 - 3); I4 += 4)
608: {
609: if (Z[I4 + o_z] <= TOL2 * Z[I4 - 3 + o_z] || Z[I4 - 1 + o_z] <= TOL2 * SIGMA)
610: {
611: Z[I4 - 1 + o_z] = - SIGMA;
612: SPLT = I4 / 4;
613: QMAX = ZERO;
614: EMIN = Z[I4 + 3 + o_z];
615: OLDEMN = Z[I4 + 4 + o_z];
616: }
617: else
618: {
619: QMAX = Math.Max(QMAX, Z[I4 + 1 + o_z]);
620: EMIN = Math.Min(EMIN, Z[I4 - 1 + o_z]);
621: OLDEMN = Math.Min(OLDEMN, Z[I4 + o_z]);
622: }
623: }
624: Z[4 * N0 - 1 + o_z] = EMIN;
625: Z[4 * N0 + o_z] = OLDEMN;
626: I0 = SPLT + 1;
627: }
628: }
629: // *
630: }
631: // *
632: INFO = 2;
633: return;
634: // *
635: // * end IWHILB
636: // *
637: LABEL130:;
638: // *
639: }
640: // *
641: INFO = 3;
642: return;
643: // *
644: // * end IWHILA
645: // *
646: LABEL150:;
647: // *
648: // * Move q's to the front.
649: // *
650: for (K = 2; K <= N; K++)
651: {
652: Z[K + o_z] = Z[4 * K - 3 + o_z];
653: }
654: // *
655: // * Sort and compute sum of eigenvalues.
656: // *
657: this._dlasrt.Run("D", N, ref Z, offset_z, ref IINFO);
658: // *
659: E = ZERO;
660: for (K = N; K >= 1; K += - 1)
661: {
662: E = E + Z[K + o_z];
663: }
664: // *
665: // * Store trace, sum(eigenvalues) and information on performance.
666: // *
667: Z[2 * N + 1 + o_z] = TRACE;
668: Z[2 * N + 2 + o_z] = E;
669: Z[2 * N + 3 + o_z] = Convert.ToDouble(ITER);
670: Z[2 * N + 4 + o_z] = Convert.ToDouble(NDIV) / Convert.ToDouble(Math.Pow(N,2));
671: Z[2 * N + 5 + o_z] = HUNDRD * NFAIL / Convert.ToDouble(ITER);
672: return;
673: // *
674: // * End of DLASQ2
675: // *
676:
677: #endregion
678:
679: }
680: }
681: }