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 auxiliary routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLAHQR is an auxiliary routine called by DHSEQR to update the
27: /// eigenvalues and Schur decomposition already computed by DHSEQR, by
28: /// dealing with the Hessenberg submatrix in rows and columns ILO to
29: /// IHI.
30: ///
31: ///</summary>
32: public class DLAHQR
33: {
34:
35:
36: #region Dependencies
37:
38: DLAMCH _dlamch; DCOPY _dcopy; DLABAD _dlabad; DLANV2 _dlanv2; DLARFG _dlarfg; DROT _drot;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const int ITMAX = 30; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0;
46: const double DAT1 = 3.0E0 / 4.0E0;const double DAT2 = - 0.4375E0; double AA = 0; double AB = 0; double BA = 0;
47: double BB = 0;double CS = 0; double DET = 0; double H11 = 0; double H12 = 0; double H21 = 0; double H21S = 0;
48: double H22 = 0;double RT1I = 0; double RT1R = 0; double RT2I = 0; double RT2R = 0; double RTDISC = 0; double S = 0;
49: double SAFMAX = 0;double SAFMIN = 0; double SMLNUM = 0; double SN = 0; double SUM = 0; double T1 = 0; double T2 = 0;
50: double T3 = 0;double TR = 0; double TST = 0; double ULP = 0; double V2 = 0; double V3 = 0; int I = 0; int I1 = 0;
51: int I2 = 0;int ITS = 0; int J = 0; int K = 0; int L = 0; int M = 0; int NH = 0; int NR = 0; int NZ = 0;
52: double[] V = new double[3]; int offset_v = 0; int o_v = -1;
53:
54: #endregion
55:
56: public DLAHQR(DLAMCH dlamch, DCOPY dcopy, DLABAD dlabad, DLANV2 dlanv2, DLARFG dlarfg, DROT drot)
57: {
58:
59:
60: #region Set Dependencies
61:
62: this._dlamch = dlamch; this._dcopy = dcopy; this._dlabad = dlabad; this._dlanv2 = dlanv2; this._dlarfg = dlarfg;
63: this._drot = drot;
64:
65: #endregion
66:
67: }
68:
69: public DLAHQR()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: LSAME lsame = new LSAME();
76: DLAMC3 dlamc3 = new DLAMC3();
77: DCOPY dcopy = new DCOPY();
78: DLABAD dlabad = new DLABAD();
79: DLAPY2 dlapy2 = new DLAPY2();
80: DNRM2 dnrm2 = new DNRM2();
81: DSCAL dscal = new DSCAL();
82: DROT drot = new DROT();
83: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
84: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
85: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
86: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
87: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
88: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
89: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
90:
91: #endregion
92:
93:
94: #region Set Dependencies
95:
96: this._dlamch = dlamch; this._dcopy = dcopy; this._dlabad = dlabad; this._dlanv2 = dlanv2; this._dlarfg = dlarfg;
97: this._drot = drot;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DLAHQR is an auxiliary routine called by DHSEQR to update the
107: /// eigenvalues and Schur decomposition already computed by DHSEQR, by
108: /// dealing with the Hessenberg submatrix in rows and columns ILO to
109: /// IHI.
110: ///
111: ///</summary>
112: /// <param name="WANTT">
113: /// (input) LOGICAL
114: /// = .TRUE. : the full Schur form T is required;
115: /// = .FALSE.: only eigenvalues are required.
116: ///</param>
117: /// <param name="WANTZ">
118: /// (input) LOGICAL
119: /// = .TRUE. : the matrix of Schur vectors Z is required;
120: /// = .FALSE.: Schur vectors are not required.
121: ///</param>
122: /// <param name="N">
123: /// (input) INTEGER
124: /// The order of the matrix H. N .GE. 0.
125: ///</param>
126: /// <param name="ILO">
127: /// (input) INTEGER
128: ///</param>
129: /// <param name="IHI">
130: /// (input) INTEGER
131: /// It is assumed that H is already upper quasi-triangular in
132: /// rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
133: /// ILO = 1). DLAHQR works primarily with the Hessenberg
134: /// submatrix in rows and columns ILO to IHI, but applies
135: /// transformations to all of H if WANTT is .TRUE..
136: /// 1 .LE. ILO .LE. max(1,IHI); IHI .LE. N.
137: ///</param>
138: /// <param name="H">
139: /// (input/output) DOUBLE PRECISION array, dimension (LDH,N)
140: /// On entry, the upper Hessenberg matrix H.
141: /// On exit, if INFO is zero and if WANTT is .TRUE., H is upper
142: /// quasi-triangular in rows and columns ILO:IHI, with any
143: /// 2-by-2 diagonal blocks in standard form. If INFO is zero
144: /// and WANTT is .FALSE., the contents of H are unspecified on
145: /// exit. The output state of H if INFO is nonzero is given
146: /// below under the description of INFO.
147: ///</param>
148: /// <param name="LDH">
149: /// (input) INTEGER
150: /// The leading dimension of the array H. LDH .GE. max(1,N).
151: ///</param>
152: /// <param name="WR">
153: /// (output) DOUBLE PRECISION array, dimension (N)
154: ///</param>
155: /// <param name="WI">
156: /// (output) DOUBLE PRECISION array, dimension (N)
157: /// The real and imaginary parts, respectively, of the computed
158: /// eigenvalues ILO to IHI are stored in the corresponding
159: /// elements of WR and WI. If two eigenvalues are computed as a
160: /// complex conjugate pair, they are stored in consecutive
161: /// elements of WR and WI, say the i-th and (i+1)th, with
162: /// WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., the
163: /// eigenvalues are stored in the same order as on the diagonal
164: /// of the Schur form returned in H, with WR(i) = H(i,i), and, if
165: /// H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
166: /// WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
167: ///</param>
168: /// <param name="ILOZ">
169: /// (input) INTEGER
170: ///</param>
171: /// <param name="IHIZ">
172: /// (input) INTEGER
173: /// Specify the rows of Z to which transformations must be
174: /// applied if WANTZ is .TRUE..
175: /// 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
176: ///</param>
177: /// <param name="Z">
178: /// (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
179: /// If WANTZ is .TRUE., on entry Z must contain the current
180: /// matrix Z of transformations accumulated by DHSEQR, and on
181: /// exit Z has been updated; transformations are applied only to
182: /// the submatrix Z(ILOZ:IHIZ,ILO:IHI).
183: /// If WANTZ is .FALSE., Z is not referenced.
184: ///</param>
185: /// <param name="LDZ">
186: /// (input) INTEGER
187: /// The leading dimension of the array Z. LDZ .GE. max(1,N).
188: ///</param>
189: /// <param name="INFO">
190: /// (output) INTEGER
191: /// = 0: successful exit
192: /// .GT. 0: If INFO = i, DLAHQR failed to compute all the
193: /// eigenvalues ILO to IHI in a total of 30 iterations
194: /// per eigenvalue; elements i+1:ihi of WR and WI
195: /// contain those eigenvalues which have been
196: /// successfully computed.
197: ///
198: /// If INFO .GT. 0 and WANTT is .FALSE., then on exit,
199: /// the remaining unconverged eigenvalues are the
200: /// eigenvalues of the upper Hessenberg matrix rows
201: /// and columns ILO thorugh INFO of the final, output
202: /// value of H.
203: ///
204: /// If INFO .GT. 0 and WANTT is .TRUE., then on exit
205: /// (*) (initial value of H)*U = U*(final value of H)
206: /// where U is an orthognal matrix. The final
207: /// value of H is upper Hessenberg and triangular in
208: /// rows and columns INFO+1 through IHI.
209: ///
210: /// If INFO .GT. 0 and WANTZ is .TRUE., then on exit
211: /// (final value of Z) = (initial value of Z)*U
212: /// where U is the orthogonal matrix in (*)
213: /// (regardless of the value of WANTT.)
214: ///</param>
215: public void Run(bool WANTT, bool WANTZ, int N, int ILO, int IHI, ref double[] H, int offset_h
216: , int LDH, ref double[] WR, int offset_wr, ref double[] WI, int offset_wi, int ILOZ, int IHIZ, ref double[] Z, int offset_z
217: , int LDZ, ref int INFO)
218: {
219:
220: #region Array Index Correction
221:
222: int o_h = -1 - LDH + offset_h; int o_wr = -1 + offset_wr; int o_wi = -1 + offset_wi;
223: int o_z = -1 - LDZ + offset_z;
224:
225: #endregion
226:
227:
228: #region Prolog
229:
230: // *
231: // * -- LAPACK auxiliary routine (version 3.1) --
232: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
233: // * November 2006
234: // *
235: // * .. Scalar Arguments ..
236: // * ..
237: // * .. Array Arguments ..
238: // * ..
239: // *
240: // * Purpose
241: // * =======
242: // *
243: // * DLAHQR is an auxiliary routine called by DHSEQR to update the
244: // * eigenvalues and Schur decomposition already computed by DHSEQR, by
245: // * dealing with the Hessenberg submatrix in rows and columns ILO to
246: // * IHI.
247: // *
248: // * Arguments
249: // * =========
250: // *
251: // * WANTT (input) LOGICAL
252: // * = .TRUE. : the full Schur form T is required;
253: // * = .FALSE.: only eigenvalues are required.
254: // *
255: // * WANTZ (input) LOGICAL
256: // * = .TRUE. : the matrix of Schur vectors Z is required;
257: // * = .FALSE.: Schur vectors are not required.
258: // *
259: // * N (input) INTEGER
260: // * The order of the matrix H. N >= 0.
261: // *
262: // * ILO (input) INTEGER
263: // * IHI (input) INTEGER
264: // * It is assumed that H is already upper quasi-triangular in
265: // * rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
266: // * ILO = 1). DLAHQR works primarily with the Hessenberg
267: // * submatrix in rows and columns ILO to IHI, but applies
268: // * transformations to all of H if WANTT is .TRUE..
269: // * 1 <= ILO <= max(1,IHI); IHI <= N.
270: // *
271: // * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
272: // * On entry, the upper Hessenberg matrix H.
273: // * On exit, if INFO is zero and if WANTT is .TRUE., H is upper
274: // * quasi-triangular in rows and columns ILO:IHI, with any
275: // * 2-by-2 diagonal blocks in standard form. If INFO is zero
276: // * and WANTT is .FALSE., the contents of H are unspecified on
277: // * exit. The output state of H if INFO is nonzero is given
278: // * below under the description of INFO.
279: // *
280: // * LDH (input) INTEGER
281: // * The leading dimension of the array H. LDH >= max(1,N).
282: // *
283: // * WR (output) DOUBLE PRECISION array, dimension (N)
284: // * WI (output) DOUBLE PRECISION array, dimension (N)
285: // * The real and imaginary parts, respectively, of the computed
286: // * eigenvalues ILO to IHI are stored in the corresponding
287: // * elements of WR and WI. If two eigenvalues are computed as a
288: // * complex conjugate pair, they are stored in consecutive
289: // * elements of WR and WI, say the i-th and (i+1)th, with
290: // * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
291: // * eigenvalues are stored in the same order as on the diagonal
292: // * of the Schur form returned in H, with WR(i) = H(i,i), and, if
293: // * H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
294: // * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
295: // *
296: // * ILOZ (input) INTEGER
297: // * IHIZ (input) INTEGER
298: // * Specify the rows of Z to which transformations must be
299: // * applied if WANTZ is .TRUE..
300: // * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
301: // *
302: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
303: // * If WANTZ is .TRUE., on entry Z must contain the current
304: // * matrix Z of transformations accumulated by DHSEQR, and on
305: // * exit Z has been updated; transformations are applied only to
306: // * the submatrix Z(ILOZ:IHIZ,ILO:IHI).
307: // * If WANTZ is .FALSE., Z is not referenced.
308: // *
309: // * LDZ (input) INTEGER
310: // * The leading dimension of the array Z. LDZ >= max(1,N).
311: // *
312: // * INFO (output) INTEGER
313: // * = 0: successful exit
314: // * .GT. 0: If INFO = i, DLAHQR failed to compute all the
315: // * eigenvalues ILO to IHI in a total of 30 iterations
316: // * per eigenvalue; elements i+1:ihi of WR and WI
317: // * contain those eigenvalues which have been
318: // * successfully computed.
319: // *
320: // * If INFO .GT. 0 and WANTT is .FALSE., then on exit,
321: // * the remaining unconverged eigenvalues are the
322: // * eigenvalues of the upper Hessenberg matrix rows
323: // * and columns ILO thorugh INFO of the final, output
324: // * value of H.
325: // *
326: // * If INFO .GT. 0 and WANTT is .TRUE., then on exit
327: // * (*) (initial value of H)*U = U*(final value of H)
328: // * where U is an orthognal matrix. The final
329: // * value of H is upper Hessenberg and triangular in
330: // * rows and columns INFO+1 through IHI.
331: // *
332: // * If INFO .GT. 0 and WANTZ is .TRUE., then on exit
333: // * (final value of Z) = (initial value of Z)*U
334: // * where U is the orthogonal matrix in (*)
335: // * (regardless of the value of WANTT.)
336: // *
337: // * Further Details
338: // * ===============
339: // *
340: // * 02-96 Based on modifications by
341: // * David Day, Sandia National Laboratory, USA
342: // *
343: // * 12-04 Further modifications by
344: // * Ralph Byers, University of Kansas, USA
345: // *
346: // * This is a modified version of DLAHQR from LAPACK version 3.0.
347: // * It is (1) more robust against overflow and underflow and
348: // * (2) adopts the more conservative Ahues & Tisseur stopping
349: // * criterion (LAWN 122, 1997).
350: // *
351: // * =========================================================
352: // *
353: // * .. Parameters ..
354: // * ..
355: // * .. Local Scalars ..
356: // * ..
357: // * .. Local Arrays ..
358: // * ..
359: // * .. External Functions ..
360: // * ..
361: // * .. External Subroutines ..
362: // * ..
363: // * .. Intrinsic Functions ..
364: // INTRINSIC ABS, DBLE, MAX, MIN, SQRT;
365: // * ..
366: // * .. Executable Statements ..
367: // *
368:
369: #endregion
370:
371:
372: #region Body
373:
374: INFO = 0;
375: // *
376: // * Quick return if possible
377: // *
378: if (N == 0) return;
379: if (ILO == IHI)
380: {
381: WR[ILO + o_wr] = H[ILO+ILO * LDH + o_h];
382: WI[ILO + o_wi] = ZERO;
383: return;
384: }
385: // *
386: // * ==== clear out the trash ====
387: for (J = ILO; J <= IHI - 3; J++)
388: {
389: H[J + 2+J * LDH + o_h] = ZERO;
390: H[J + 3+J * LDH + o_h] = ZERO;
391: }
392: if (ILO <= IHI - 2) H[IHI+(IHI - 2) * LDH + o_h] = ZERO;
393: // *
394: NH = IHI - ILO + 1;
395: NZ = IHIZ - ILOZ + 1;
396: // *
397: // * Set machine-dependent constants for the stopping criterion.
398: // *
399: SAFMIN = this._dlamch.Run("SAFE MINIMUM");
400: SAFMAX = ONE / SAFMIN;
401: this._dlabad.Run(ref SAFMIN, ref SAFMAX);
402: ULP = this._dlamch.Run("PRECISION");
403: SMLNUM = SAFMIN * (Convert.ToDouble(NH) / ULP);
404: // *
405: // * I1 and I2 are the indices of the first row and last column of H
406: // * to which transformations must be applied. If eigenvalues only are
407: // * being computed, I1 and I2 are set inside the main loop.
408: // *
409: if (WANTT)
410: {
411: I1 = 1;
412: I2 = N;
413: }
414: // *
415: // * The main loop begins here. I is the loop index and decreases from
416: // * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
417: // * with the active submatrix in rows and columns L to I.
418: // * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
419: // * H(L,L-1) is negligible so that the matrix splits.
420: // *
421: I = IHI;
422: LABEL20:;
423: L = ILO;
424: if (I < ILO) goto LABEL160;
425: // *
426: // * Perform QR iterations on rows and columns ILO to I until a
427: // * submatrix of order 1 or 2 splits off at the bottom because a
428: // * subdiagonal element has become negligible.
429: // *
430: for (ITS = 0; ITS <= ITMAX; ITS++)
431: {
432: // *
433: // * Look for a single small subdiagonal element.
434: // *
435: for (K = I; K >= L + 1; K += - 1)
436: {
437: if (Math.Abs(H[K+(K - 1) * LDH + o_h]) <= SMLNUM) goto LABEL40;
438: TST = Math.Abs(H[K - 1+(K - 1) * LDH + o_h]) + Math.Abs(H[K+K * LDH + o_h]);
439: if (TST == ZERO)
440: {
441: if (K - 2 >= ILO) TST = TST + Math.Abs(H[K - 1+(K - 2) * LDH + o_h]);
442: if (K + 1 <= IHI) TST = TST + Math.Abs(H[K + 1+K * LDH + o_h]);
443: }
444: // * ==== The following is a conservative small subdiagonal
445: // * . deflation criterion due to Ahues & Tisseur (LAWN 122,
446: // * . 1997). It has better mathematical foundation and
447: // * . improves accuracy in some cases. ====
448: if (Math.Abs(H[K+(K - 1) * LDH + o_h]) <= ULP * TST)
449: {
450: AB = Math.Max(Math.Abs(H[K+(K - 1) * LDH + o_h]), Math.Abs(H[K - 1+K * LDH + o_h]));
451: BA = Math.Min(Math.Abs(H[K+(K - 1) * LDH + o_h]), Math.Abs(H[K - 1+K * LDH + o_h]));
452: AA = Math.Max(Math.Abs(H[K+K * LDH + o_h]), Math.Abs(H[K - 1+(K - 1) * LDH + o_h] - H[K+K * LDH + o_h]));
453: BB = Math.Min(Math.Abs(H[K+K * LDH + o_h]), Math.Abs(H[K - 1+(K - 1) * LDH + o_h] - H[K+K * LDH + o_h]));
454: S = AA + AB;
455: if (BA * (AB / S) <= Math.Max(SMLNUM, ULP * (BB * (AA / S)))) goto LABEL40;
456: }
457: }
458: LABEL40:;
459: L = K;
460: if (L > ILO)
461: {
462: // *
463: // * H(L,L-1) is negligible
464: // *
465: H[L+(L - 1) * LDH + o_h] = ZERO;
466: }
467: // *
468: // * Exit from loop if a submatrix of order 1 or 2 has split off.
469: // *
470: if (L >= I - 1) goto LABEL150;
471: // *
472: // * Now the active submatrix is in rows and columns L to I. If
473: // * eigenvalues only are being computed, only the active submatrix
474: // * need be transformed.
475: // *
476: if (!WANTT)
477: {
478: I1 = L;
479: I2 = I;
480: }
481: // *
482: if (ITS == 10 || ITS == 20)
483: {
484: // *
485: // * Exceptional shift.
486: // *
487: H11 = DAT1 * S + H[I+I * LDH + o_h];
488: H12 = DAT2 * S;
489: H21 = S;
490: H22 = H11;
491: }
492: else
493: {
494: // *
495: // * Prepare to use Francis' double shift
496: // * (i.e. 2nd degree generalized Rayleigh quotient)
497: // *
498: H11 = H[I - 1+(I - 1) * LDH + o_h];
499: H21 = H[I+(I - 1) * LDH + o_h];
500: H12 = H[I - 1+I * LDH + o_h];
501: H22 = H[I+I * LDH + o_h];
502: }
503: S = Math.Abs(H11) + Math.Abs(H12) + Math.Abs(H21) + Math.Abs(H22);
504: if (S == ZERO)
505: {
506: RT1R = ZERO;
507: RT1I = ZERO;
508: RT2R = ZERO;
509: RT2I = ZERO;
510: }
511: else
512: {
513: H11 = H11 / S;
514: H21 = H21 / S;
515: H12 = H12 / S;
516: H22 = H22 / S;
517: TR = (H11 + H22) / TWO;
518: DET = (H11 - TR) * (H22 - TR) - H12 * H21;
519: RTDISC = Math.Sqrt(Math.Abs(DET));
520: if (DET >= ZERO)
521: {
522: // *
523: // * ==== complex conjugate shifts ====
524: // *
525: RT1R = TR * S;
526: RT2R = RT1R;
527: RT1I = RTDISC * S;
528: RT2I = - RT1I;
529: }
530: else
531: {
532: // *
533: // * ==== real shifts (use only one of them) ====
534: // *
535: RT1R = TR + RTDISC;
536: RT2R = TR - RTDISC;
537: if (Math.Abs(RT1R - H22) <= Math.Abs(RT2R - H22))
538: {
539: RT1R = RT1R * S;
540: RT2R = RT1R;
541: }
542: else
543: {
544: RT2R = RT2R * S;
545: RT1R = RT2R;
546: }
547: RT1I = ZERO;
548: RT2I = ZERO;
549: }
550: }
551: // *
552: // * Look for two consecutive small subdiagonal elements.
553: // *
554: for (M = I - 2; M >= L; M += - 1)
555: {
556: // * Determine the effect of starting the double-shift QR
557: // * iteration at row M, and see if this would make H(M,M-1)
558: // * negligible. (The following uses scaling to avoid
559: // * overflows and most underflows.)
560: // *
561: H21S = H[M + 1+M * LDH + o_h];
562: S = Math.Abs(H[M+M * LDH + o_h] - RT2R) + Math.Abs(RT2I) + Math.Abs(H21S);
563: H21S = H[M + 1+M * LDH + o_h] / S;
564: V[1 + o_v] = H21S * H[M+(M + 1) * LDH + o_h] + (H[M+M * LDH + o_h] - RT1R) * ((H[M+M * LDH + o_h] - RT2R) / S) - RT1I * (RT2I / S);
565: V[2 + o_v] = H21S * (H[M+M * LDH + o_h] + H[M + 1+(M + 1) * LDH + o_h] - RT1R - RT2R);
566: V[3 + o_v] = H21S * H[M + 2+(M + 1) * LDH + o_h];
567: S = Math.Abs(V[1 + o_v]) + Math.Abs(V[2 + o_v]) + Math.Abs(V[3 + o_v]);
568: V[1 + o_v] = V[1 + o_v] / S;
569: V[2 + o_v] = V[2 + o_v] / S;
570: V[3 + o_v] = V[3 + o_v] / S;
571: if (M == L) goto LABEL60;
572: if (Math.Abs(H[M+(M - 1) * LDH + o_h]) * (Math.Abs(V[2 + o_v]) + Math.Abs(V[3 + o_v])) <= ULP * Math.Abs(V[1 + o_v]) * (Math.Abs(H[M - 1+(M - 1) * LDH + o_h]) + Math.Abs(H[M+M * LDH + o_h]) + Math.Abs(H[M + 1+(M + 1) * LDH + o_h]))) goto LABEL60;
573: }
574: LABEL60:;
575: // *
576: // * Double-shift QR step
577: // *
578: for (K = M; K <= I - 1; K++)
579: {
580: // *
581: // * The first iteration of this loop determines a reflection G
582: // * from the vector V and applies it from left and right to H,
583: // * thus creating a nonzero bulge below the subdiagonal.
584: // *
585: // * Each subsequent iteration determines a reflection G to
586: // * restore the Hessenberg form in the (K-1)th column, and thus
587: // * chases the bulge one step toward the bottom of the active
588: // * submatrix. NR is the order of G.
589: // *
590: NR = Math.Min(3, I - K + 1);
591: if (K > M) this._dcopy.Run(NR, H, K+(K - 1) * LDH + o_h, 1, ref V, offset_v, 1);
592: this._dlarfg.Run(NR, ref V[1 + o_v], ref V, 2 + o_v, 1, ref T1);
593: if (K > M)
594: {
595: H[K+(K - 1) * LDH + o_h] = V[1 + o_v];
596: H[K + 1+(K - 1) * LDH + o_h] = ZERO;
597: if (K < I - 1) H[K + 2+(K - 1) * LDH + o_h] = ZERO;
598: }
599: else
600: {
601: if (M > L)
602: {
603: H[K+(K - 1) * LDH + o_h] = - H[K+(K - 1) * LDH + o_h];
604: }
605: }
606: V2 = V[2 + o_v];
607: T2 = T1 * V2;
608: if (NR == 3)
609: {
610: V3 = V[3 + o_v];
611: T3 = T1 * V3;
612: // *
613: // * Apply G from the left to transform the rows of the matrix
614: // * in columns K to I2.
615: // *
616: for (J = K; J <= I2; J++)
617: {
618: SUM = H[K+J * LDH + o_h] + V2 * H[K + 1+J * LDH + o_h] + V3 * H[K + 2+J * LDH + o_h];
619: H[K+J * LDH + o_h] = H[K+J * LDH + o_h] - SUM * T1;
620: H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - SUM * T2;
621: H[K + 2+J * LDH + o_h] = H[K + 2+J * LDH + o_h] - SUM * T3;
622: }
623: // *
624: // * Apply G from the right to transform the columns of the
625: // * matrix in rows I1 to min(K+3,I).
626: // *
627: for (J = I1; J <= Math.Min(K + 3, I); J++)
628: {
629: SUM = H[J+K * LDH + o_h] + V2 * H[J+(K + 1) * LDH + o_h] + V3 * H[J+(K + 2) * LDH + o_h];
630: H[J+K * LDH + o_h] = H[J+K * LDH + o_h] - SUM * T1;
631: H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - SUM * T2;
632: H[J+(K + 2) * LDH + o_h] = H[J+(K + 2) * LDH + o_h] - SUM * T3;
633: }
634: // *
635: if (WANTZ)
636: {
637: // *
638: // * Accumulate transformations in the matrix Z
639: // *
640: for (J = ILOZ; J <= IHIZ; J++)
641: {
642: SUM = Z[J+K * LDZ + o_z] + V2 * Z[J+(K + 1) * LDZ + o_z] + V3 * Z[J+(K + 2) * LDZ + o_z];
643: Z[J+K * LDZ + o_z] = Z[J+K * LDZ + o_z] - SUM * T1;
644: Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - SUM * T2;
645: Z[J+(K + 2) * LDZ + o_z] = Z[J+(K + 2) * LDZ + o_z] - SUM * T3;
646: }
647: }
648: }
649: else
650: {
651: if (NR == 2)
652: {
653: // *
654: // * Apply G from the left to transform the rows of the matrix
655: // * in columns K to I2.
656: // *
657: for (J = K; J <= I2; J++)
658: {
659: SUM = H[K+J * LDH + o_h] + V2 * H[K + 1+J * LDH + o_h];
660: H[K+J * LDH + o_h] = H[K+J * LDH + o_h] - SUM * T1;
661: H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - SUM * T2;
662: }
663: // *
664: // * Apply G from the right to transform the columns of the
665: // * matrix in rows I1 to min(K+3,I).
666: // *
667: for (J = I1; J <= I; J++)
668: {
669: SUM = H[J+K * LDH + o_h] + V2 * H[J+(K + 1) * LDH + o_h];
670: H[J+K * LDH + o_h] = H[J+K * LDH + o_h] - SUM * T1;
671: H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - SUM * T2;
672: }
673: // *
674: if (WANTZ)
675: {
676: // *
677: // * Accumulate transformations in the matrix Z
678: // *
679: for (J = ILOZ; J <= IHIZ; J++)
680: {
681: SUM = Z[J+K * LDZ + o_z] + V2 * Z[J+(K + 1) * LDZ + o_z];
682: Z[J+K * LDZ + o_z] = Z[J+K * LDZ + o_z] - SUM * T1;
683: Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - SUM * T2;
684: }
685: }
686: }
687: }
688: }
689: // *
690: }
691: // *
692: // * Failure to converge in remaining number of iterations
693: // *
694: INFO = I;
695: return;
696: // *
697: LABEL150:;
698: // *
699: if (L == I)
700: {
701: // *
702: // * H(I,I-1) is negligible: one eigenvalue has converged.
703: // *
704: WR[I + o_wr] = H[I+I * LDH + o_h];
705: WI[I + o_wi] = ZERO;
706: }
707: else
708: {
709: if (L == I - 1)
710: {
711: // *
712: // * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
713: // *
714: // * Transform the 2-by-2 submatrix to standard Schur form,
715: // * and compute and store the eigenvalues.
716: // *
717: this._dlanv2.Run(ref H[I - 1+(I - 1) * LDH + o_h], ref H[I - 1+I * LDH + o_h], ref H[I+(I - 1) * LDH + o_h], ref H[I+I * LDH + o_h], ref WR[I - 1 + o_wr], ref WI[I - 1 + o_wi]
718: , ref WR[I + o_wr], ref WI[I + o_wi], ref CS, ref SN);
719: // *
720: if (WANTT)
721: {
722: // *
723: // * Apply the transformation to the rest of H.
724: // *
725: if (I2 > I)
726: {
727: this._drot.Run(I2 - I, ref H, I - 1+(I + 1) * LDH + o_h, LDH, ref H, I+(I + 1) * LDH + o_h, LDH, CS
728: , SN);
729: }
730: this._drot.Run(I - I1 - 1, ref H, I1+(I - 1) * LDH + o_h, 1, ref H, I1+I * LDH + o_h, 1, CS
731: , SN);
732: }
733: if (WANTZ)
734: {
735: // *
736: // * Apply the transformation to Z.
737: // *
738: this._drot.Run(NZ, ref Z, ILOZ+(I - 1) * LDZ + o_z, 1, ref Z, ILOZ+I * LDZ + o_z, 1, CS
739: , SN);
740: }
741: }
742: }
743: // *
744: // * return to start of the main loop with new value of I.
745: // *
746: I = L - 1;
747: goto LABEL20;
748: // *
749: LABEL160:;
750: return;
751: // *
752: // * End of DLAHQR
753: // *
754:
755: #endregion
756:
757: }
758: }
759: }