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: /// DTREVC computes some or all of the right and/or left eigenvectors of
27: /// a real upper quasi-triangular matrix T.
28: /// Matrices of this type are produced by the Schur factorization of
29: /// a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
30: ///
31: /// The right eigenvector x and the left eigenvector y of T corresponding
32: /// to an eigenvalue w are defined by:
33: ///
34: /// T*x = w*x, (y**H)*T = w*(y**H)
35: ///
36: /// where y**H denotes the conjugate transpose of y.
37: /// The eigenvalues are not input to this routine, but are read directly
38: /// from the diagonal blocks of T.
39: ///
40: /// This routine returns the matrices X and/or Y of right and left
41: /// eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
42: /// input matrix. If Q is the orthogonal factor that reduces a matrix
43: /// A to Schur form T, then Q*X and Q*Y are the matrices of right and
44: /// left eigenvectors of A.
45: ///
46: ///</summary>
47: public class DTREVC
48: {
49:
50:
51: #region Dependencies
52:
53: LSAME _lsame; IDAMAX _idamax; DDOT _ddot; DLAMCH _dlamch; DAXPY _daxpy; DCOPY _dcopy; DGEMV _dgemv; DLALN2 _dlaln2;
54: DSCAL _dscal;XERBLA _xerbla; DLABAD _dlabad;
55:
56: #endregion
57:
58:
59: #region Fields
60:
61: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool ALLV = false; bool BOTHV = false; bool LEFTV = false;
62: bool OVER = false;bool PAIR = false; bool RIGHTV = false; bool SOMEV = false; int I = 0; int IERR = 0; int II = 0;
63: int IP = 0;int IS = 0; int J = 0; int J1 = 0; int J2 = 0; int JNXT = 0; int K = 0; int KI = 0; int N2 = 0;
64: double BETA = 0;double BIGNUM = 0; double EMAX = 0; double OVFL = 0; double REC = 0; double REMAX = 0; double SCALE = 0;
65: double SMIN = 0;double SMLNUM = 0; double ULP = 0; double UNFL = 0; double VCRIT = 0; double VMAX = 0; double WI = 0;
66: double WR = 0;double XNORM = 0; double[] X = new double[2 * 2]; int offset_x = 0; int o_x = -3;
67:
68: #endregion
69:
70: public DTREVC(LSAME lsame, IDAMAX idamax, DDOT ddot, DLAMCH dlamch, DAXPY daxpy, DCOPY dcopy, DGEMV dgemv, DLALN2 dlaln2, DSCAL dscal, XERBLA xerbla
71: , DLABAD dlabad)
72: {
73:
74:
75: #region Set Dependencies
76:
77: this._lsame = lsame; this._idamax = idamax; this._ddot = ddot; this._dlamch = dlamch; this._daxpy = daxpy;
78: this._dcopy = dcopy;this._dgemv = dgemv; this._dlaln2 = dlaln2; this._dscal = dscal; this._xerbla = xerbla;
79: this._dlabad = dlabad;
80:
81: #endregion
82:
83: }
84:
85: public DTREVC()
86: {
87:
88:
89: #region Dependencies (Initialization)
90:
91: LSAME lsame = new LSAME();
92: IDAMAX idamax = new IDAMAX();
93: DDOT ddot = new DDOT();
94: DLAMC3 dlamc3 = new DLAMC3();
95: DAXPY daxpy = new DAXPY();
96: DCOPY dcopy = new DCOPY();
97: XERBLA xerbla = new XERBLA();
98: DLADIV dladiv = new DLADIV();
99: DSCAL dscal = new DSCAL();
100: DLABAD dlabad = new DLABAD();
101: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
102: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
103: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
104: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
105: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
106: DGEMV dgemv = new DGEMV(lsame, xerbla);
107: DLALN2 dlaln2 = new DLALN2(dlamch, dladiv);
108:
109: #endregion
110:
111:
112: #region Set Dependencies
113:
114: this._lsame = lsame; this._idamax = idamax; this._ddot = ddot; this._dlamch = dlamch; this._daxpy = daxpy;
115: this._dcopy = dcopy;this._dgemv = dgemv; this._dlaln2 = dlaln2; this._dscal = dscal; this._xerbla = xerbla;
116: this._dlabad = dlabad;
117:
118: #endregion
119:
120: }
121: /// <summary>
122: /// Purpose
123: /// =======
124: ///
125: /// DTREVC computes some or all of the right and/or left eigenvectors of
126: /// a real upper quasi-triangular matrix T.
127: /// Matrices of this type are produced by the Schur factorization of
128: /// a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
129: ///
130: /// The right eigenvector x and the left eigenvector y of T corresponding
131: /// to an eigenvalue w are defined by:
132: ///
133: /// T*x = w*x, (y**H)*T = w*(y**H)
134: ///
135: /// where y**H denotes the conjugate transpose of y.
136: /// The eigenvalues are not input to this routine, but are read directly
137: /// from the diagonal blocks of T.
138: ///
139: /// This routine returns the matrices X and/or Y of right and left
140: /// eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
141: /// input matrix. If Q is the orthogonal factor that reduces a matrix
142: /// A to Schur form T, then Q*X and Q*Y are the matrices of right and
143: /// left eigenvectors of A.
144: ///
145: ///</summary>
146: /// <param name="SIDE">
147: /// (input) CHARACTER*1
148: /// = 'R': compute right eigenvectors only;
149: /// = 'L': compute left eigenvectors only;
150: /// = 'B': compute both right and left eigenvectors.
151: ///</param>
152: /// <param name="HOWMNY">
153: /// (input) CHARACTER*1
154: /// = 'A': compute all right and/or left eigenvectors;
155: /// = 'B': compute all right and/or left eigenvectors,
156: /// backtransformed by the matrices in VR and/or VL;
157: /// = 'S': compute selected right and/or left eigenvectors,
158: /// as indicated by the logical array SELECT.
159: ///</param>
160: /// <param name="SELECT">
161: /// (input/output) LOGICAL array, dimension (N)
162: /// If HOWMNY = 'S', SELECT specifies the eigenvectors to be
163: /// computed.
164: /// If w(j) is a real eigenvalue, the corresponding real
165: /// eigenvector is computed if SELECT(j) is .TRUE..
166: /// If w(j) and w(j+1) are the real and imaginary parts of a
167: /// complex eigenvalue, the corresponding complex eigenvector is
168: /// computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
169: /// on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
170: /// .FALSE..
171: /// Not referenced if HOWMNY = 'A' or 'B'.
172: ///</param>
173: /// <param name="N">
174: /// (input) INTEGER
175: /// The order of the matrix T. N .GE. 0.
176: ///</param>
177: /// <param name="T">
178: /// (input) DOUBLE PRECISION array, dimension (LDT,N)
179: /// The upper quasi-triangular matrix T in Schur canonical form.
180: ///</param>
181: /// <param name="LDT">
182: /// (input) INTEGER
183: /// The leading dimension of the array T. LDT .GE. max(1,N).
184: ///</param>
185: /// <param name="VL">
186: /// (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
187: /// On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
188: /// contain an N-by-N matrix Q (usually the orthogonal matrix Q
189: /// of Schur vectors returned by DHSEQR).
190: /// On exit, if SIDE = 'L' or 'B', VL contains:
191: /// if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
192: /// if HOWMNY = 'B', the matrix Q*Y;
193: /// if HOWMNY = 'S', the left eigenvectors of T specified by
194: /// SELECT, stored consecutively in the columns
195: /// of VL, in the same order as their
196: /// eigenvalues.
197: /// A complex eigenvector corresponding to a complex eigenvalue
198: /// is stored in two consecutive columns, the first holding the
199: /// real part, and the second the imaginary part.
200: /// Not referenced if SIDE = 'R'.
201: ///</param>
202: /// <param name="LDVL">
203: /// (input) INTEGER
204: /// The leading dimension of the array VL. LDVL .GE. 1, and if
205: /// SIDE = 'L' or 'B', LDVL .GE. N.
206: ///</param>
207: /// <param name="VR">
208: /// (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
209: /// On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
210: /// contain an N-by-N matrix Q (usually the orthogonal matrix Q
211: /// of Schur vectors returned by DHSEQR).
212: /// On exit, if SIDE = 'R' or 'B', VR contains:
213: /// if HOWMNY = 'A', the matrix X of right eigenvectors of T;
214: /// if HOWMNY = 'B', the matrix Q*X;
215: /// if HOWMNY = 'S', the right eigenvectors of T specified by
216: /// SELECT, stored consecutively in the columns
217: /// of VR, in the same order as their
218: /// eigenvalues.
219: /// A complex eigenvector corresponding to a complex eigenvalue
220: /// is stored in two consecutive columns, the first holding the
221: /// real part and the second the imaginary part.
222: /// Not referenced if SIDE = 'L'.
223: ///</param>
224: /// <param name="LDVR">
225: /// (input) INTEGER
226: /// The leading dimension of the array VR. LDVR .GE. 1, and if
227: /// SIDE = 'R' or 'B', LDVR .GE. N.
228: ///</param>
229: /// <param name="MM">
230: /// (input) INTEGER
231: /// The number of columns in the arrays VL and/or VR. MM .GE. M.
232: ///</param>
233: /// <param name="M">
234: /// (output) INTEGER
235: /// The number of columns in the arrays VL and/or VR actually
236: /// used to store the eigenvectors.
237: /// If HOWMNY = 'A' or 'B', M is set to N.
238: /// Each selected real eigenvector occupies one column and each
239: /// selected complex eigenvector occupies two columns.
240: ///</param>
241: /// <param name="WORK">
242: /// (workspace) DOUBLE PRECISION array, dimension (3*N)
243: ///</param>
244: /// <param name="INFO">
245: /// (output) INTEGER
246: /// = 0: successful exit
247: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
248: ///</param>
249: public void Run(string SIDE, string HOWMNY, ref bool[] SELECT, int offset_select, int N, double[] T, int offset_t, int LDT
250: , ref double[] VL, int offset_vl, int LDVL, ref double[] VR, int offset_vr, int LDVR, int MM, ref int M
251: , ref double[] WORK, int offset_work, ref int INFO)
252: {
253:
254: #region Array Index Correction
255:
256: int o_select = -1 + offset_select; int o_t = -1 - LDT + offset_t; int o_vl = -1 - LDVL + offset_vl;
257: int o_vr = -1 - LDVR + offset_vr; int o_work = -1 + offset_work;
258:
259: #endregion
260:
261:
262: #region Strings
263:
264: SIDE = SIDE.Substring(0, 1); HOWMNY = HOWMNY.Substring(0, 1);
265:
266: #endregion
267:
268:
269: #region Prolog
270:
271: // *
272: // * -- LAPACK routine (version 3.1) --
273: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
274: // * November 2006
275: // *
276: // * .. Scalar Arguments ..
277: // * ..
278: // * .. Array Arguments ..
279: // * ..
280: // *
281: // * Purpose
282: // * =======
283: // *
284: // * DTREVC computes some or all of the right and/or left eigenvectors of
285: // * a real upper quasi-triangular matrix T.
286: // * Matrices of this type are produced by the Schur factorization of
287: // * a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
288: // *
289: // * The right eigenvector x and the left eigenvector y of T corresponding
290: // * to an eigenvalue w are defined by:
291: // *
292: // * T*x = w*x, (y**H)*T = w*(y**H)
293: // *
294: // * where y**H denotes the conjugate transpose of y.
295: // * The eigenvalues are not input to this routine, but are read directly
296: // * from the diagonal blocks of T.
297: // *
298: // * This routine returns the matrices X and/or Y of right and left
299: // * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
300: // * input matrix. If Q is the orthogonal factor that reduces a matrix
301: // * A to Schur form T, then Q*X and Q*Y are the matrices of right and
302: // * left eigenvectors of A.
303: // *
304: // * Arguments
305: // * =========
306: // *
307: // * SIDE (input) CHARACTER*1
308: // * = 'R': compute right eigenvectors only;
309: // * = 'L': compute left eigenvectors only;
310: // * = 'B': compute both right and left eigenvectors.
311: // *
312: // * HOWMNY (input) CHARACTER*1
313: // * = 'A': compute all right and/or left eigenvectors;
314: // * = 'B': compute all right and/or left eigenvectors,
315: // * backtransformed by the matrices in VR and/or VL;
316: // * = 'S': compute selected right and/or left eigenvectors,
317: // * as indicated by the logical array SELECT.
318: // *
319: // * SELECT (input/output) LOGICAL array, dimension (N)
320: // * If HOWMNY = 'S', SELECT specifies the eigenvectors to be
321: // * computed.
322: // * If w(j) is a real eigenvalue, the corresponding real
323: // * eigenvector is computed if SELECT(j) is .TRUE..
324: // * If w(j) and w(j+1) are the real and imaginary parts of a
325: // * complex eigenvalue, the corresponding complex eigenvector is
326: // * computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
327: // * on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
328: // * .FALSE..
329: // * Not referenced if HOWMNY = 'A' or 'B'.
330: // *
331: // * N (input) INTEGER
332: // * The order of the matrix T. N >= 0.
333: // *
334: // * T (input) DOUBLE PRECISION array, dimension (LDT,N)
335: // * The upper quasi-triangular matrix T in Schur canonical form.
336: // *
337: // * LDT (input) INTEGER
338: // * The leading dimension of the array T. LDT >= max(1,N).
339: // *
340: // * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
341: // * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
342: // * contain an N-by-N matrix Q (usually the orthogonal matrix Q
343: // * of Schur vectors returned by DHSEQR).
344: // * On exit, if SIDE = 'L' or 'B', VL contains:
345: // * if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
346: // * if HOWMNY = 'B', the matrix Q*Y;
347: // * if HOWMNY = 'S', the left eigenvectors of T specified by
348: // * SELECT, stored consecutively in the columns
349: // * of VL, in the same order as their
350: // * eigenvalues.
351: // * A complex eigenvector corresponding to a complex eigenvalue
352: // * is stored in two consecutive columns, the first holding the
353: // * real part, and the second the imaginary part.
354: // * Not referenced if SIDE = 'R'.
355: // *
356: // * LDVL (input) INTEGER
357: // * The leading dimension of the array VL. LDVL >= 1, and if
358: // * SIDE = 'L' or 'B', LDVL >= N.
359: // *
360: // * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
361: // * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
362: // * contain an N-by-N matrix Q (usually the orthogonal matrix Q
363: // * of Schur vectors returned by DHSEQR).
364: // * On exit, if SIDE = 'R' or 'B', VR contains:
365: // * if HOWMNY = 'A', the matrix X of right eigenvectors of T;
366: // * if HOWMNY = 'B', the matrix Q*X;
367: // * if HOWMNY = 'S', the right eigenvectors of T specified by
368: // * SELECT, stored consecutively in the columns
369: // * of VR, in the same order as their
370: // * eigenvalues.
371: // * A complex eigenvector corresponding to a complex eigenvalue
372: // * is stored in two consecutive columns, the first holding the
373: // * real part and the second the imaginary part.
374: // * Not referenced if SIDE = 'L'.
375: // *
376: // * LDVR (input) INTEGER
377: // * The leading dimension of the array VR. LDVR >= 1, and if
378: // * SIDE = 'R' or 'B', LDVR >= N.
379: // *
380: // * MM (input) INTEGER
381: // * The number of columns in the arrays VL and/or VR. MM >= M.
382: // *
383: // * M (output) INTEGER
384: // * The number of columns in the arrays VL and/or VR actually
385: // * used to store the eigenvectors.
386: // * If HOWMNY = 'A' or 'B', M is set to N.
387: // * Each selected real eigenvector occupies one column and each
388: // * selected complex eigenvector occupies two columns.
389: // *
390: // * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
391: // *
392: // * INFO (output) INTEGER
393: // * = 0: successful exit
394: // * < 0: if INFO = -i, the i-th argument had an illegal value
395: // *
396: // * Further Details
397: // * ===============
398: // *
399: // * The algorithm used in this program is basically backward (forward)
400: // * substitution, with scaling to make the the code robust against
401: // * possible overflow.
402: // *
403: // * Each eigenvector is normalized so that the element of largest
404: // * magnitude has magnitude 1; here the magnitude of a complex number
405: // * (x,y) is taken to be |x| + |y|.
406: // *
407: // * =====================================================================
408: // *
409: // * .. Parameters ..
410: // * ..
411: // * .. Local Scalars ..
412: // * ..
413: // * .. External Functions ..
414: // * ..
415: // * .. External Subroutines ..
416: // * ..
417: // * .. Intrinsic Functions ..
418: // INTRINSIC ABS, MAX, SQRT;
419: // * ..
420: // * .. Local Arrays ..
421: // * ..
422: // * .. Executable Statements ..
423: // *
424: // * Decode and test the input parameters
425: // *
426:
427: #endregion
428:
429:
430: #region Body
431:
432: BOTHV = this._lsame.Run(SIDE, "B");
433: RIGHTV = this._lsame.Run(SIDE, "R") || BOTHV;
434: LEFTV = this._lsame.Run(SIDE, "L") || BOTHV;
435: // *
436: ALLV = this._lsame.Run(HOWMNY, "A");
437: OVER = this._lsame.Run(HOWMNY, "B");
438: SOMEV = this._lsame.Run(HOWMNY, "S");
439: // *
440: INFO = 0;
441: if (!RIGHTV && !LEFTV)
442: {
443: INFO = - 1;
444: }
445: else
446: {
447: if (!ALLV && !OVER && !SOMEV)
448: {
449: INFO = - 2;
450: }
451: else
452: {
453: if (N < 0)
454: {
455: INFO = - 4;
456: }
457: else
458: {
459: if (LDT < Math.Max(1, N))
460: {
461: INFO = - 6;
462: }
463: else
464: {
465: if (LDVL < 1 || (LEFTV && LDVL < N))
466: {
467: INFO = - 8;
468: }
469: else
470: {
471: if (LDVR < 1 || (RIGHTV && LDVR < N))
472: {
473: INFO = - 10;
474: }
475: else
476: {
477: // *
478: // * Set M to the number of columns required to store the selected
479: // * eigenvectors, standardize the array SELECT if necessary, and
480: // * test MM.
481: // *
482: if (SOMEV)
483: {
484: M = 0;
485: PAIR = false;
486: for (J = 1; J <= N; J++)
487: {
488: if (PAIR)
489: {
490: PAIR = false;
491: SELECT[J + o_select] = false;
492: }
493: else
494: {
495: if (J < N)
496: {
497: if (T[J + 1+J * LDT + o_t] == ZERO)
498: {
499: if (SELECT[J + o_select]) M = M + 1;
500: }
501: else
502: {
503: PAIR = true;
504: if (SELECT[J + o_select] || SELECT[J + 1 + o_select])
505: {
506: SELECT[J + o_select] = true;
507: M = M + 2;
508: }
509: }
510: }
511: else
512: {
513: if (SELECT[N + o_select]) M = M + 1;
514: }
515: }
516: }
517: }
518: else
519: {
520: M = N;
521: }
522: // *
523: if (MM < M)
524: {
525: INFO = - 11;
526: }
527: }
528: }
529: }
530: }
531: }
532: }
533: if (INFO != 0)
534: {
535: this._xerbla.Run("DTREVC", - INFO);
536: return;
537: }
538: // *
539: // * Quick return if possible.
540: // *
541: if (N == 0) return;
542: // *
543: // * Set the constants to control overflow.
544: // *
545: UNFL = this._dlamch.Run("Safe minimum");
546: OVFL = ONE / UNFL;
547: this._dlabad.Run(ref UNFL, ref OVFL);
548: ULP = this._dlamch.Run("Precision");
549: SMLNUM = UNFL * (N / ULP);
550: BIGNUM = (ONE - ULP) / SMLNUM;
551: // *
552: // * Compute 1-norm of each column of strictly upper triangular
553: // * part of T to control overflow in triangular solver.
554: // *
555: WORK[1 + o_work] = ZERO;
556: for (J = 2; J <= N; J++)
557: {
558: WORK[J + o_work] = ZERO;
559: for (I = 1; I <= J - 1; I++)
560: {
561: WORK[J + o_work] = WORK[J + o_work] + Math.Abs(T[I+J * LDT + o_t]);
562: }
563: }
564: // *
565: // * Index IP is used to specify the real or complex eigenvalue:
566: // * IP = 0, real eigenvalue,
567: // * 1, first of conjugate complex pair: (wr,wi)
568: // * -1, second of conjugate complex pair: (wr,wi)
569: // *
570: N2 = 2 * N;
571: // *
572: if (RIGHTV)
573: {
574: // *
575: // * Compute right eigenvectors.
576: // *
577: IP = 0;
578: IS = M;
579: for (KI = N; KI >= 1; KI += - 1)
580: {
581: // *
582: if (IP == 1) goto LABEL130;
583: if (KI == 1) goto LABEL40;
584: if (T[KI+(KI - 1) * LDT + o_t] == ZERO) goto LABEL40;
585: IP = - 1;
586: // *
587: LABEL40:;
588: if (SOMEV)
589: {
590: if (IP == 0)
591: {
592: if (!SELECT[KI + o_select]) goto LABEL130;
593: }
594: else
595: {
596: if (!SELECT[KI - 1 + o_select]) goto LABEL130;
597: }
598: }
599: // *
600: // * Compute the KI-th eigenvalue (WR,WI).
601: // *
602: WR = T[KI+KI * LDT + o_t];
603: WI = ZERO;
604: if (IP != 0) WI = Math.Sqrt(Math.Abs(T[KI+(KI - 1) * LDT + o_t])) * Math.Sqrt(Math.Abs(T[KI - 1+KI * LDT + o_t]));
605: SMIN = Math.Max(ULP * (Math.Abs(WR) + Math.Abs(WI)), SMLNUM);
606: // *
607: if (IP == 0)
608: {
609: // *
610: // * Real right eigenvector
611: // *
612: WORK[KI + N + o_work] = ONE;
613: // *
614: // * Form right-hand side
615: // *
616: for (K = 1; K <= KI - 1; K++)
617: {
618: WORK[K + N + o_work] = - T[K+KI * LDT + o_t];
619: }
620: // *
621: // * Solve the upper quasi-triangular system:
622: // * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
623: // *
624: JNXT = KI - 1;
625: for (J = KI - 1; J >= 1; J += - 1)
626: {
627: if (J > JNXT) goto LABEL60;
628: J1 = J;
629: J2 = J;
630: JNXT = J - 1;
631: if (J > 1)
632: {
633: if (T[J+(J - 1) * LDT + o_t] != ZERO)
634: {
635: J1 = J - 1;
636: JNXT = J - 2;
637: }
638: }
639: // *
640: if (J1 == J2)
641: {
642: // *
643: // * 1-by-1 diagonal block
644: // *
645: this._dlaln2.Run(false, 1, 1, SMIN, ONE, T, J+J * LDT + o_t
646: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
647: , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
648: // *
649: // * Scale X(1,1) to avoid overflow when updating
650: // * the right-hand side.
651: // *
652: if (XNORM > ONE)
653: {
654: if (WORK[J + o_work] > BIGNUM / XNORM)
655: {
656: X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;
657: SCALE = SCALE / XNORM;
658: }
659: }
660: // *
661: // * Scale if necessary
662: // *
663: if (SCALE != ONE) this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);
664: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
665: // *
666: // * Update right-hand side
667: // *
668: this._daxpy.Run(J - 1, - X[1+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
669: // *
670: }
671: else
672: {
673: // *
674: // * 2-by-2 diagonal block
675: // *
676: this._dlaln2.Run(false, 2, 1, SMIN, ONE, T, J - 1+(J - 1) * LDT + o_t
677: , LDT, ONE, ONE, WORK, J - 1 + N + o_work, N, WR
678: , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
679: // *
680: // * Scale X(1,1) and X(2,1) to avoid overflow when
681: // * updating the right-hand side.
682: // *
683: if (XNORM > ONE)
684: {
685: BETA = Math.Max(WORK[J - 1 + o_work], WORK[J + o_work]);
686: if (BETA > BIGNUM / XNORM)
687: {
688: X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;
689: X[2+1 * 2 + o_x] = X[2+1 * 2 + o_x] / XNORM;
690: SCALE = SCALE / XNORM;
691: }
692: }
693: // *
694: // * Scale if necessary
695: // *
696: if (SCALE != ONE) this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);
697: WORK[J - 1 + N + o_work] = X[1+1 * 2 + o_x];
698: WORK[J + N + o_work] = X[2+1 * 2 + o_x];
699: // *
700: // * Update right-hand side
701: // *
702: this._daxpy.Run(J - 2, - X[1+1 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
703: this._daxpy.Run(J - 2, - X[2+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
704: }
705: LABEL60:;
706: }
707: // *
708: // * Copy the vector x or Q*x to VR and normalize.
709: // *
710: if (!OVER)
711: {
712: this._dcopy.Run(KI, WORK, 1 + N + o_work, 1, ref VR, 1+IS * LDVR + o_vr, 1);
713: // *
714: II = this._idamax.Run(KI, VR, 1+IS * LDVR + o_vr, 1);
715: REMAX = ONE / Math.Abs(VR[II+IS * LDVR + o_vr]);
716: this._dscal.Run(KI, REMAX, ref VR, 1+IS * LDVR + o_vr, 1);
717: // *
718: for (K = KI + 1; K <= N; K++)
719: {
720: VR[K+IS * LDVR + o_vr] = ZERO;
721: }
722: }
723: else
724: {
725: if (KI > 1)
726: {
727: this._dgemv.Run("N", N, KI - 1, ONE, VR, offset_vr, LDVR
728: , WORK, 1 + N + o_work, 1, WORK[KI + N + o_work], ref VR, 1+KI * LDVR + o_vr, 1);
729: }
730: // *
731: II = this._idamax.Run(N, VR, 1+KI * LDVR + o_vr, 1);
732: REMAX = ONE / Math.Abs(VR[II+KI * LDVR + o_vr]);
733: this._dscal.Run(N, REMAX, ref VR, 1+KI * LDVR + o_vr, 1);
734: }
735: // *
736: }
737: else
738: {
739: // *
740: // * Complex right eigenvector.
741: // *
742: // * Initial solve
743: // * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
744: // * [ (T(KI,KI-1) T(KI,KI) ) ]
745: // *
746: if (Math.Abs(T[KI - 1+KI * LDT + o_t]) >= Math.Abs(T[KI+(KI - 1) * LDT + o_t]))
747: {
748: WORK[KI - 1 + N + o_work] = ONE;
749: WORK[KI + N2 + o_work] = WI / T[KI - 1+KI * LDT + o_t];
750: }
751: else
752: {
753: WORK[KI - 1 + N + o_work] = - WI / T[KI+(KI - 1) * LDT + o_t];
754: WORK[KI + N2 + o_work] = ONE;
755: }
756: WORK[KI + N + o_work] = ZERO;
757: WORK[KI - 1 + N2 + o_work] = ZERO;
758: // *
759: // * Form right-hand side
760: // *
761: for (K = 1; K <= KI - 2; K++)
762: {
763: WORK[K + N + o_work] = - WORK[KI - 1 + N + o_work] * T[K+(KI - 1) * LDT + o_t];
764: WORK[K + N2 + o_work] = - WORK[KI + N2 + o_work] * T[K+KI * LDT + o_t];
765: }
766: // *
767: // * Solve upper quasi-triangular system:
768: // * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
769: // *
770: JNXT = KI - 2;
771: for (J = KI - 2; J >= 1; J += - 1)
772: {
773: if (J > JNXT) goto LABEL90;
774: J1 = J;
775: J2 = J;
776: JNXT = J - 1;
777: if (J > 1)
778: {
779: if (T[J+(J - 1) * LDT + o_t] != ZERO)
780: {
781: J1 = J - 1;
782: JNXT = J - 2;
783: }
784: }
785: // *
786: if (J1 == J2)
787: {
788: // *
789: // * 1-by-1 diagonal block
790: // *
791: this._dlaln2.Run(false, 1, 2, SMIN, ONE, T, J+J * LDT + o_t
792: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
793: , WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
794: // *
795: // * Scale X(1,1) and X(1,2) to avoid overflow when
796: // * updating the right-hand side.
797: // *
798: if (XNORM > ONE)
799: {
800: if (WORK[J + o_work] > BIGNUM / XNORM)
801: {
802: X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;
803: X[1+2 * 2 + o_x] = X[1+2 * 2 + o_x] / XNORM;
804: SCALE = SCALE / XNORM;
805: }
806: }
807: // *
808: // * Scale if necessary
809: // *
810: if (SCALE != ONE)
811: {
812: this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);
813: this._dscal.Run(KI, SCALE, ref WORK, 1 + N2 + o_work, 1);
814: }
815: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
816: WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];
817: // *
818: // * Update the right-hand side
819: // *
820: this._daxpy.Run(J - 1, - X[1+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
821: this._daxpy.Run(J - 1, - X[1+2 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);
822: // *
823: }
824: else
825: {
826: // *
827: // * 2-by-2 diagonal block
828: // *
829: this._dlaln2.Run(false, 2, 2, SMIN, ONE, T, J - 1+(J - 1) * LDT + o_t
830: , LDT, ONE, ONE, WORK, J - 1 + N + o_work, N, WR
831: , WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
832: // *
833: // * Scale X to avoid overflow when updating
834: // * the right-hand side.
835: // *
836: if (XNORM > ONE)
837: {
838: BETA = Math.Max(WORK[J - 1 + o_work], WORK[J + o_work]);
839: if (BETA > BIGNUM / XNORM)
840: {
841: REC = ONE / XNORM;
842: X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] * REC;
843: X[1+2 * 2 + o_x] = X[1+2 * 2 + o_x] * REC;
844: X[2+1 * 2 + o_x] = X[2+1 * 2 + o_x] * REC;
845: X[2+2 * 2 + o_x] = X[2+2 * 2 + o_x] * REC;
846: SCALE = SCALE * REC;
847: }
848: }
849: // *
850: // * Scale if necessary
851: // *
852: if (SCALE != ONE)
853: {
854: this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);
855: this._dscal.Run(KI, SCALE, ref WORK, 1 + N2 + o_work, 1);
856: }
857: WORK[J - 1 + N + o_work] = X[1+1 * 2 + o_x];
858: WORK[J + N + o_work] = X[2+1 * 2 + o_x];
859: WORK[J - 1 + N2 + o_work] = X[1+2 * 2 + o_x];
860: WORK[J + N2 + o_work] = X[2+2 * 2 + o_x];
861: // *
862: // * Update the right-hand side
863: // *
864: this._daxpy.Run(J - 2, - X[1+1 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
865: this._daxpy.Run(J - 2, - X[2+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);
866: this._daxpy.Run(J - 2, - X[1+2 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);
867: this._daxpy.Run(J - 2, - X[2+2 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);
868: }
869: LABEL90:;
870: }
871: // *
872: // * Copy the vector x or Q*x to VR and normalize.
873: // *
874: if (!OVER)
875: {
876: this._dcopy.Run(KI, WORK, 1 + N + o_work, 1, ref VR, 1+(IS - 1) * LDVR + o_vr, 1);
877: this._dcopy.Run(KI, WORK, 1 + N2 + o_work, 1, ref VR, 1+IS * LDVR + o_vr, 1);
878: // *
879: EMAX = ZERO;
880: for (K = 1; K <= KI; K++)
881: {
882: EMAX = Math.Max(EMAX, Math.Abs(VR[K+(IS - 1) * LDVR + o_vr]) + Math.Abs(VR[K+IS * LDVR + o_vr]));
883: }
884: // *
885: REMAX = ONE / EMAX;
886: this._dscal.Run(KI, REMAX, ref VR, 1+(IS - 1) * LDVR + o_vr, 1);
887: this._dscal.Run(KI, REMAX, ref VR, 1+IS * LDVR + o_vr, 1);
888: // *
889: for (K = KI + 1; K <= N; K++)
890: {
891: VR[K+(IS - 1) * LDVR + o_vr] = ZERO;
892: VR[K+IS * LDVR + o_vr] = ZERO;
893: }
894: // *
895: }
896: else
897: {
898: // *
899: if (KI > 2)
900: {
901: this._dgemv.Run("N", N, KI - 2, ONE, VR, offset_vr, LDVR
902: , WORK, 1 + N + o_work, 1, WORK[KI - 1 + N + o_work], ref VR, 1+(KI - 1) * LDVR + o_vr, 1);
903: this._dgemv.Run("N", N, KI - 2, ONE, VR, offset_vr, LDVR
904: , WORK, 1 + N2 + o_work, 1, WORK[KI + N2 + o_work], ref VR, 1+KI * LDVR + o_vr, 1);
905: }
906: else
907: {
908: this._dscal.Run(N, WORK[KI - 1 + N + o_work], ref VR, 1+(KI - 1) * LDVR + o_vr, 1);
909: this._dscal.Run(N, WORK[KI + N2 + o_work], ref VR, 1+KI * LDVR + o_vr, 1);
910: }
911: // *
912: EMAX = ZERO;
913: for (K = 1; K <= N; K++)
914: {
915: EMAX = Math.Max(EMAX, Math.Abs(VR[K+(KI - 1) * LDVR + o_vr]) + Math.Abs(VR[K+KI * LDVR + o_vr]));
916: }
917: REMAX = ONE / EMAX;
918: this._dscal.Run(N, REMAX, ref VR, 1+(KI - 1) * LDVR + o_vr, 1);
919: this._dscal.Run(N, REMAX, ref VR, 1+KI * LDVR + o_vr, 1);
920: }
921: }
922: // *
923: IS = IS - 1;
924: if (IP != 0) IS = IS - 1;
925: LABEL130:;
926: if (IP == 1) IP = 0;
927: if (IP == - 1) IP = 1;
928: }
929: }
930: // *
931: if (LEFTV)
932: {
933: // *
934: // * Compute left eigenvectors.
935: // *
936: IP = 0;
937: IS = 1;
938: for (KI = 1; KI <= N; KI++)
939: {
940: // *
941: if (IP == - 1) goto LABEL250;
942: if (KI == N) goto LABEL150;
943: if (T[KI + 1+KI * LDT + o_t] == ZERO) goto LABEL150;
944: IP = 1;
945: // *
946: LABEL150:;
947: if (SOMEV)
948: {
949: if (!SELECT[KI + o_select]) goto LABEL250;
950: }
951: // *
952: // * Compute the KI-th eigenvalue (WR,WI).
953: // *
954: WR = T[KI+KI * LDT + o_t];
955: WI = ZERO;
956: if (IP != 0) WI = Math.Sqrt(Math.Abs(T[KI+(KI + 1) * LDT + o_t])) * Math.Sqrt(Math.Abs(T[KI + 1+KI * LDT + o_t]));
957: SMIN = Math.Max(ULP * (Math.Abs(WR) + Math.Abs(WI)), SMLNUM);
958: // *
959: if (IP == 0)
960: {
961: // *
962: // * Real left eigenvector.
963: // *
964: WORK[KI + N + o_work] = ONE;
965: // *
966: // * Form right-hand side
967: // *
968: for (K = KI + 1; K <= N; K++)
969: {
970: WORK[K + N + o_work] = - T[KI+K * LDT + o_t];
971: }
972: // *
973: // * Solve the quasi-triangular system:
974: // * (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
975: // *
976: VMAX = ONE;
977: VCRIT = BIGNUM;
978: // *
979: JNXT = KI + 1;
980: for (J = KI + 1; J <= N; J++)
981: {
982: if (J < JNXT) goto LABEL170;
983: J1 = J;
984: J2 = J;
985: JNXT = J + 1;
986: if (J < N)
987: {
988: if (T[J + 1+J * LDT + o_t] != ZERO)
989: {
990: J2 = J + 1;
991: JNXT = J + 2;
992: }
993: }
994: // *
995: if (J1 == J2)
996: {
997: // *
998: // * 1-by-1 diagonal block
999: // *
1000: // * Scale if necessary to avoid overflow when forming
1001: // * the right-hand side.
1002: // *
1003: if (WORK[J + o_work] > VCRIT)
1004: {
1005: REC = ONE / VMAX;
1006: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);
1007: VMAX = ONE;
1008: VCRIT = BIGNUM;
1009: }
1010: // *
1011: WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+J * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);
1012: // *
1013: // * Solve (T(J,J)-WR)'*X = WORK
1014: // *
1015: this._dlaln2.Run(false, 1, 1, SMIN, ONE, T, J+J * LDT + o_t
1016: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
1017: , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
1018: // *
1019: // * Scale if necessary
1020: // *
1021: if (SCALE != ONE) this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);
1022: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
1023: VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), VMAX);
1024: VCRIT = BIGNUM / VMAX;
1025: // *
1026: }
1027: else
1028: {
1029: // *
1030: // * 2-by-2 diagonal block
1031: // *
1032: // * Scale if necessary to avoid overflow when forming
1033: // * the right-hand side.
1034: // *
1035: BETA = Math.Max(WORK[J + o_work], WORK[J + 1 + o_work]);
1036: if (BETA > VCRIT)
1037: {
1038: REC = ONE / VMAX;
1039: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);
1040: VMAX = ONE;
1041: VCRIT = BIGNUM;
1042: }
1043: // *
1044: WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+J * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);
1045: // *
1046: WORK[J + 1 + N + o_work] = WORK[J + 1 + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+(J + 1) * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);
1047: // *
1048: // * Solve
1049: // * [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 )
1050: // * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
1051: // *
1052: this._dlaln2.Run(true, 2, 1, SMIN, ONE, T, J+J * LDT + o_t
1053: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
1054: , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
1055: // *
1056: // * Scale if necessary
1057: // *
1058: if (SCALE != ONE) this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);
1059: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
1060: WORK[J + 1 + N + o_work] = X[2+1 * 2 + o_x];
1061: // *
1062: VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), Math.Max(Math.Abs(WORK[J + 1 + N + o_work]), VMAX));
1063: VCRIT = BIGNUM / VMAX;
1064: // *
1065: }
1066: LABEL170:;
1067: }
1068: // *
1069: // * Copy the vector x or Q*x to VL and normalize.
1070: // *
1071: if (!OVER)
1072: {
1073: this._dcopy.Run(N - KI + 1, WORK, KI + N + o_work, 1, ref VL, KI+IS * LDVL + o_vl, 1);
1074: // *
1075: II = this._idamax.Run(N - KI + 1, VL, KI+IS * LDVL + o_vl, 1) + KI - 1;
1076: REMAX = ONE / Math.Abs(VL[II+IS * LDVL + o_vl]);
1077: this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+IS * LDVL + o_vl, 1);
1078: // *
1079: for (K = 1; K <= KI - 1; K++)
1080: {
1081: VL[K+IS * LDVL + o_vl] = ZERO;
1082: }
1083: // *
1084: }
1085: else
1086: {
1087: // *
1088: if (KI < N)
1089: {
1090: this._dgemv.Run("N", N, N - KI, ONE, VL, 1+(KI + 1) * LDVL + o_vl, LDVL
1091: , WORK, KI + 1 + N + o_work, 1, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);
1092: }
1093: // *
1094: II = this._idamax.Run(N, VL, 1+KI * LDVL + o_vl, 1);
1095: REMAX = ONE / Math.Abs(VL[II+KI * LDVL + o_vl]);
1096: this._dscal.Run(N, REMAX, ref VL, 1+KI * LDVL + o_vl, 1);
1097: // *
1098: }
1099: // *
1100: }
1101: else
1102: {
1103: // *
1104: // * Complex left eigenvector.
1105: // *
1106: // * Initial solve:
1107: // * ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0.
1108: // * ((T(KI+1,KI) T(KI+1,KI+1)) )
1109: // *
1110: if (Math.Abs(T[KI+(KI + 1) * LDT + o_t]) >= Math.Abs(T[KI + 1+KI * LDT + o_t]))
1111: {
1112: WORK[KI + N + o_work] = WI / T[KI+(KI + 1) * LDT + o_t];
1113: WORK[KI + 1 + N2 + o_work] = ONE;
1114: }
1115: else
1116: {
1117: WORK[KI + N + o_work] = ONE;
1118: WORK[KI + 1 + N2 + o_work] = - WI / T[KI + 1+KI * LDT + o_t];
1119: }
1120: WORK[KI + 1 + N + o_work] = ZERO;
1121: WORK[KI + N2 + o_work] = ZERO;
1122: // *
1123: // * Form right-hand side
1124: // *
1125: for (K = KI + 2; K <= N; K++)
1126: {
1127: WORK[K + N + o_work] = - WORK[KI + N + o_work] * T[KI+K * LDT + o_t];
1128: WORK[K + N2 + o_work] = - WORK[KI + 1 + N2 + o_work] * T[KI + 1+K * LDT + o_t];
1129: }
1130: // *
1131: // * Solve complex quasi-triangular system:
1132: // * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
1133: // *
1134: VMAX = ONE;
1135: VCRIT = BIGNUM;
1136: // *
1137: JNXT = KI + 2;
1138: for (J = KI + 2; J <= N; J++)
1139: {
1140: if (J < JNXT) goto LABEL200;
1141: J1 = J;
1142: J2 = J;
1143: JNXT = J + 1;
1144: if (J < N)
1145: {
1146: if (T[J + 1+J * LDT + o_t] != ZERO)
1147: {
1148: J2 = J + 1;
1149: JNXT = J + 2;
1150: }
1151: }
1152: // *
1153: if (J1 == J2)
1154: {
1155: // *
1156: // * 1-by-1 diagonal block
1157: // *
1158: // * Scale if necessary to avoid overflow when
1159: // * forming the right-hand side elements.
1160: // *
1161: if (WORK[J + o_work] > VCRIT)
1162: {
1163: REC = ONE / VMAX;
1164: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);
1165: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N2 + o_work, 1);
1166: VMAX = ONE;
1167: VCRIT = BIGNUM;
1168: }
1169: // *
1170: WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);
1171: WORK[J + N2 + o_work] = WORK[J + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);
1172: // *
1173: // * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
1174: // *
1175: this._dlaln2.Run(false, 1, 2, SMIN, ONE, T, J+J * LDT + o_t
1176: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
1177: , - WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
1178: // *
1179: // * Scale if necessary
1180: // *
1181: if (SCALE != ONE)
1182: {
1183: this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);
1184: this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N2 + o_work, 1);
1185: }
1186: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
1187: WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];
1188: VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), Math.Max(Math.Abs(WORK[J + N2 + o_work]), VMAX));
1189: VCRIT = BIGNUM / VMAX;
1190: // *
1191: }
1192: else
1193: {
1194: // *
1195: // * 2-by-2 diagonal block
1196: // *
1197: // * Scale if necessary to avoid overflow when forming
1198: // * the right-hand side elements.
1199: // *
1200: BETA = Math.Max(WORK[J + o_work], WORK[J + 1 + o_work]);
1201: if (BETA > VCRIT)
1202: {
1203: REC = ONE / VMAX;
1204: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);
1205: this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N2 + o_work, 1);
1206: VMAX = ONE;
1207: VCRIT = BIGNUM;
1208: }
1209: // *
1210: WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);
1211: // *
1212: WORK[J + N2 + o_work] = WORK[J + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);
1213: // *
1214: WORK[J + 1 + N + o_work] = WORK[J + 1 + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+(J + 1) * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);
1215: // *
1216: WORK[J + 1 + N2 + o_work] = WORK[J + 1 + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+(J + 1) * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);
1217: // *
1218: // * Solve 2-by-2 complex linear equation
1219: // * ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B
1220: // * ([T(j+1,j) T(j+1,j+1)] )
1221: // *
1222: this._dlaln2.Run(true, 2, 2, SMIN, ONE, T, J+J * LDT + o_t
1223: , LDT, ONE, ONE, WORK, J + N + o_work, N, WR
1224: , - WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);
1225: // *
1226: // * Scale if necessary
1227: // *
1228: if (SCALE != ONE)
1229: {
1230: this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);
1231: this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N2 + o_work, 1);
1232: }
1233: WORK[J + N + o_work] = X[1+1 * 2 + o_x];
1234: WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];
1235: WORK[J + 1 + N + o_work] = X[2+1 * 2 + o_x];
1236: WORK[J + 1 + N2 + o_work] = X[2+2 * 2 + o_x];
1237: VMAX = Math.Max(Math.Abs(X[1+1 * 2 + o_x]), Math.Max(Math.Abs(X[1+2 * 2 + o_x]), Math.Max(Math.Abs(X[2+1 * 2 + o_x]), Math.Max(Math.Abs(X[2+2 * 2 + o_x]), VMAX))));
1238: VCRIT = BIGNUM / VMAX;
1239: // *
1240: }
1241: LABEL200:;
1242: }
1243: // *
1244: // * Copy the vector x or Q*x to VL and normalize.
1245: // *
1246: if (!OVER)
1247: {
1248: this._dcopy.Run(N - KI + 1, WORK, KI + N + o_work, 1, ref VL, KI+IS * LDVL + o_vl, 1);
1249: this._dcopy.Run(N - KI + 1, WORK, KI + N2 + o_work, 1, ref VL, KI+(IS + 1) * LDVL + o_vl, 1);
1250: // *
1251: EMAX = ZERO;
1252: for (K = KI; K <= N; K++)
1253: {
1254: EMAX = Math.Max(EMAX, Math.Abs(VL[K+IS * LDVL + o_vl]) + Math.Abs(VL[K+(IS + 1) * LDVL + o_vl]));
1255: }
1256: REMAX = ONE / EMAX;
1257: this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+IS * LDVL + o_vl, 1);
1258: this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+(IS + 1) * LDVL + o_vl, 1);
1259: // *
1260: for (K = 1; K <= KI - 1; K++)
1261: {
1262: VL[K+IS * LDVL + o_vl] = ZERO;
1263: VL[K+(IS + 1) * LDVL + o_vl] = ZERO;
1264: }
1265: }
1266: else
1267: {
1268: if (KI < N - 1)
1269: {
1270: this._dgemv.Run("N", N, N - KI - 1, ONE, VL, 1+(KI + 2) * LDVL + o_vl, LDVL
1271: , WORK, KI + 2 + N + o_work, 1, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);
1272: this._dgemv.Run("N", N, N - KI - 1, ONE, VL, 1+(KI + 2) * LDVL + o_vl, LDVL
1273: , WORK, KI + 2 + N2 + o_work, 1, WORK[KI + 1 + N2 + o_work], ref VL, 1+(KI + 1) * LDVL + o_vl, 1);
1274: }
1275: else
1276: {
1277: this._dscal.Run(N, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);
1278: this._dscal.Run(N, WORK[KI + 1 + N2 + o_work], ref VL, 1+(KI + 1) * LDVL + o_vl, 1);
1279: }
1280: // *
1281: EMAX = ZERO;
1282: for (K = 1; K <= N; K++)
1283: {
1284: EMAX = Math.Max(EMAX, Math.Abs(VL[K+KI * LDVL + o_vl]) + Math.Abs(VL[K+(KI + 1) * LDVL + o_vl]));
1285: }
1286: REMAX = ONE / EMAX;
1287: this._dscal.Run(N, REMAX, ref VL, 1+KI * LDVL + o_vl, 1);
1288: this._dscal.Run(N, REMAX, ref VL, 1+(KI + 1) * LDVL + o_vl, 1);
1289: // *
1290: }
1291: // *
1292: }
1293: // *
1294: IS = IS + 1;
1295: if (IP != 0) IS = IS + 1;
1296: LABEL250:;
1297: if (IP == - 1) IP = 0;
1298: if (IP == 1) IP = - 1;
1299: // *
1300: }
1301: // *
1302: }
1303: // *
1304: return;
1305: // *
1306: // * End of DTREVC
1307: // *
1308:
1309: #endregion
1310:
1311: }
1312: }
1313: }