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: /// DLAQR4 computes the eigenvalues of a Hessenberg matrix H
27: /// and, optionally, the matrices T and Z from the Schur decomposition
28: /// H = Z T Z**T, where T is an upper quasi-triangular matrix (the
29: /// Schur form), and Z is the orthogonal matrix of Schur vectors.
30: ///
31: /// Optionally Z may be postmultiplied into an input orthogonal
32: /// matrix Q so that this routine can give the Schur factorization
33: /// of a matrix A which has been reduced to the Hessenberg form H
34: /// by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
35: ///
36: ///</summary>
37: public class DLAQR4
38: {
39:
40:
41: #region Dependencies
42:
43: ILAENV _ilaenv; DLACPY _dlacpy; DLAHQR _dlahqr; DLANV2 _dlanv2; DLAQR2 _dlaqr2; DLAQR5 _dlaqr5;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const int NTINY = 11; const int KEXNW = 5; const int KEXSH = 6; const double WILK1 = 0.75E0;
51: const double WILK2 = - 0.4375E0;const double ZERO = 0.0E0; const double ONE = 1.0E0; double AA = 0; double BB = 0;
52: double CC = 0;double CS = 0; double DD = 0; double SN = 0; double SS = 0; double SWAP = 0; int I = 0; int INF = 0;
53: int IT = 0;int ITMAX = 0; int K = 0; int KACC22 = 0; int KBOT = 0; int KDU = 0; int KS = 0; int KT = 0; int KTOP = 0;
54: int KU = 0;int KV = 0; int KWH = 0; int KWTOP = 0; int KWV = 0; int LD = 0; int LS = 0; int LWKOPT = 0; int NDFL = 0;
55: int NH = 0;int NHO = 0; int NIBBLE = 0; int NMIN = 0; int NS = 0; int NSMAX = 0; int NSR = 0; int NVE = 0; int NW = 0;
56: int NWMAX = 0;int NWR = 0; bool NWINC = false; bool SORTED = false; string JBCMPZ = new string(' ', 2);
57: double[] ZDUM = new double[1 * 1]; int offset_zdum = 0;
58:
59: #endregion
60:
61: public DLAQR4(ILAENV ilaenv, DLACPY dlacpy, DLAHQR dlahqr, DLANV2 dlanv2, DLAQR2 dlaqr2, DLAQR5 dlaqr5)
62: {
63:
64:
65: #region Set Dependencies
66:
67: this._ilaenv = ilaenv; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlaqr2 = dlaqr2;
68: this._dlaqr5 = dlaqr5;
69:
70: #endregion
71:
72: }
73:
74: public DLAQR4()
75: {
76:
77:
78: #region Dependencies (Initialization)
79:
80: IEEECK ieeeck = new IEEECK();
81: IPARMQ iparmq = new IPARMQ();
82: LSAME lsame = new LSAME();
83: DLAMC3 dlamc3 = new DLAMC3();
84: DCOPY dcopy = new DCOPY();
85: DLABAD dlabad = new DLABAD();
86: DLAPY2 dlapy2 = new DLAPY2();
87: DNRM2 dnrm2 = new DNRM2();
88: DSCAL dscal = new DSCAL();
89: DROT drot = new DROT();
90: DAXPY daxpy = new DAXPY();
91: XERBLA xerbla = new XERBLA();
92: DLASSQ dlassq = new DLASSQ();
93: IDAMAX idamax = new IDAMAX();
94: DSWAP dswap = new DSWAP();
95: DLAQR1 dlaqr1 = new DLAQR1();
96: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
97: DLACPY dlacpy = new DLACPY(lsame);
98: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
99: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
100: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
101: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
102: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
103: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
104: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
105: DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
106: DGEMV dgemv = new DGEMV(lsame, xerbla);
107: DGER dger = new DGER(xerbla);
108: DLARF dlarf = new DLARF(dgemv, dger, lsame);
109: DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
110: DGEMM dgemm = new DGEMM(lsame, xerbla);
111: DTRMM dtrmm = new DTRMM(lsame, xerbla);
112: DTRMV dtrmv = new DTRMV(lsame, xerbla);
113: DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
114: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
115: DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
116: DLASET dlaset = new DLASET(lsame);
117: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
118: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
119: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
120: DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
121: DLANGE dlange = new DLANGE(dlassq, lsame);
122: DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
123: DLARTG dlartg = new DLARTG(dlamch);
124: DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
125: DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
126: DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
127: DLAQR2 dlaqr2 = new DLAQR2(dlamch, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlarf, dlarfg
128: , dlaset, dorghr, dtrexc);
129: DLAQR5 dlaqr5 = new DLAQR5(dlamch, dgemm, dlabad, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm);
130:
131: #endregion
132:
133:
134: #region Set Dependencies
135:
136: this._ilaenv = ilaenv; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlaqr2 = dlaqr2;
137: this._dlaqr5 = dlaqr5;
138:
139: #endregion
140:
141: }
142: /// <summary>
143: /// Purpose
144: /// =======
145: ///
146: /// DLAQR4 computes the eigenvalues of a Hessenberg matrix H
147: /// and, optionally, the matrices T and Z from the Schur decomposition
148: /// H = Z T Z**T, where T is an upper quasi-triangular matrix (the
149: /// Schur form), and Z is the orthogonal matrix of Schur vectors.
150: ///
151: /// Optionally Z may be postmultiplied into an input orthogonal
152: /// matrix Q so that this routine can give the Schur factorization
153: /// of a matrix A which has been reduced to the Hessenberg form H
154: /// by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
155: ///
156: ///</summary>
157: /// <param name="WANTT">
158: /// (input) LOGICAL
159: /// = .TRUE. : the full Schur form T is required;
160: /// = .FALSE.: only eigenvalues are required.
161: ///</param>
162: /// <param name="WANTZ">
163: /// (input) LOGICAL
164: /// = .TRUE. : the matrix of Schur vectors Z is required;
165: /// = .FALSE.: Schur vectors are not required.
166: ///</param>
167: /// <param name="N">
168: /// (input) INTEGER
169: /// The order of the matrix H. N .GE. 0.
170: ///</param>
171: /// <param name="ILO">
172: /// (input) INTEGER
173: ///</param>
174: /// <param name="IHI">
175: /// (input) INTEGER
176: /// It is assumed that H is already upper triangular in rows
177: /// and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
178: /// H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
179: /// previous call to DGEBAL, and then passed to DGEHRD when the
180: /// matrix output by DGEBAL is reduced to Hessenberg form.
181: /// Otherwise, ILO and IHI should be set to 1 and N,
182: /// respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
183: /// If N = 0, then ILO = 1 and IHI = 0.
184: ///</param>
185: /// <param name="H">
186: /// (input/output) DOUBLE PRECISION array, dimension (LDH,N)
187: /// On entry, the upper Hessenberg matrix H.
188: /// On exit, if INFO = 0 and WANTT is .TRUE., then H contains
189: /// the upper quasi-triangular matrix T from the Schur
190: /// decomposition (the Schur form); 2-by-2 diagonal blocks
191: /// (corresponding to complex conjugate pairs of eigenvalues)
192: /// are returned in standard form, with H(i,i) = H(i+1,i+1)
193: /// and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
194: /// .FALSE., then the contents of H are unspecified on exit.
195: /// (The output value of H when INFO.GT.0 is given under the
196: /// description of INFO below.)
197: ///
198: /// This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
199: /// j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
200: ///</param>
201: /// <param name="LDH">
202: /// (input) INTEGER
203: /// The leading dimension of the array H. LDH .GE. max(1,N).
204: ///</param>
205: /// <param name="WR">
206: /// (output) DOUBLE PRECISION array, dimension (IHI)
207: ///</param>
208: /// <param name="WI">
209: /// (output) DOUBLE PRECISION array, dimension (IHI)
210: /// The real and imaginary parts, respectively, of the computed
211: /// eigenvalues of H(ILO:IHI,ILO:IHI) are stored WR(ILO:IHI)
212: /// and WI(ILO:IHI). If two eigenvalues are computed as a
213: /// complex conjugate pair, they are stored in consecutive
214: /// elements of WR and WI, say the i-th and (i+1)th, with
215: /// WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
216: /// the eigenvalues are stored in the same order as on the
217: /// diagonal of the Schur form returned in H, with
218: /// WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
219: /// block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
220: /// WI(i+1) = -WI(i).
221: ///</param>
222: /// <param name="ILOZ">
223: /// (input) INTEGER
224: ///</param>
225: /// <param name="IHIZ">
226: /// (input) INTEGER
227: /// Specify the rows of Z to which transformations must be
228: /// applied if WANTZ is .TRUE..
229: /// 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
230: ///</param>
231: /// <param name="Z">
232: /// (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
233: /// If WANTZ is .FALSE., then Z is not referenced.
234: /// If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
235: /// replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
236: /// orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
237: /// (The output value of Z when INFO.GT.0 is given under
238: /// the description of INFO below.)
239: ///</param>
240: /// <param name="LDZ">
241: /// (input) INTEGER
242: /// The leading dimension of the array Z. if WANTZ is .TRUE.
243: /// then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
244: ///</param>
245: /// <param name="WORK">
246: /// (workspace/output) DOUBLE PRECISION array, dimension LWORK
247: /// On exit, if LWORK = -1, WORK(1) returns an estimate of
248: /// the optimal value for LWORK.
249: ///</param>
250: /// <param name="LWORK">
251: /// (input) INTEGER
252: /// The dimension of the array WORK. LWORK .GE. max(1,N)
253: /// is sufficient, but LWORK typically as large as 6*N may
254: /// be required for optimal performance. A workspace query
255: /// to determine the optimal workspace size is recommended.
256: ///
257: /// If LWORK = -1, then DLAQR4 does a workspace query.
258: /// In this case, DLAQR4 checks the input parameters and
259: /// estimates the optimal workspace size for the given
260: /// values of N, ILO and IHI. The estimate is returned
261: /// in WORK(1). No error message related to LWORK is
262: /// issued by XERBLA. Neither H nor Z are accessed.
263: ///
264: ///</param>
265: /// <param name="INFO">
266: /// (output) INTEGER
267: /// = 0: successful exit
268: /// .GT. 0: if INFO = i, DLAQR4 failed to compute all of
269: /// the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
270: /// and WI contain those eigenvalues which have been
271: /// successfully computed. (Failures are rare.)
272: ///
273: /// If INFO .GT. 0 and WANT is .FALSE., then on exit,
274: /// the remaining unconverged eigenvalues are the eigen-
275: /// values of the upper Hessenberg matrix rows and
276: /// columns ILO through INFO of the final, output
277: /// value of H.
278: ///
279: /// If INFO .GT. 0 and WANTT is .TRUE., then on exit
280: ///
281: /// (*) (initial value of H)*U = U*(final value of H)
282: ///
283: /// where U is an orthogonal matrix. The final
284: /// value of H is upper Hessenberg and quasi-triangular
285: /// in rows and columns INFO+1 through IHI.
286: ///
287: /// If INFO .GT. 0 and WANTZ is .TRUE., then on exit
288: ///
289: /// (final value of Z(ILO:IHI,ILOZ:IHIZ)
290: /// = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
291: ///
292: /// where U is the orthogonal matrix in (*) (regard-
293: /// less of the value of WANTT.)
294: ///
295: /// If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
296: /// accessed.
297: ///</param>
298: public void Run(bool WANTT, bool WANTZ, int N, int ILO, int IHI, ref double[] H, int offset_h
299: , int LDH, ref double[] WR, int offset_wr, ref double[] WI, int offset_wi, int ILOZ, int IHIZ, ref double[] Z, int offset_z
300: , int LDZ, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
301: {
302:
303: #region Array Index Correction
304:
305: int o_h = -1 - LDH + offset_h; int o_wr = -1 + offset_wr; int o_wi = -1 + offset_wi;
306: int o_z = -1 - LDZ + offset_z; int o_work = -1 + offset_work;
307:
308: #endregion
309:
310:
311: #region Prolog
312:
313: // *
314: // * -- LAPACK auxiliary routine (version 3.1) --
315: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
316: // * November 2006
317: // *
318: // * .. Scalar Arguments ..
319: // * ..
320: // * .. Array Arguments ..
321: // * ..
322: // *
323: // * This subroutine implements one level of recursion for DLAQR0.
324: // * It is a complete implementation of the small bulge multi-shift
325: // * QR algorithm. It may be called by DLAQR0 and, for large enough
326: // * deflation window size, it may be called by DLAQR3. This
327: // * subroutine is identical to DLAQR0 except that it calls DLAQR2
328: // * instead of DLAQR3.
329: // *
330: // * Purpose
331: // * =======
332: // *
333: // * DLAQR4 computes the eigenvalues of a Hessenberg matrix H
334: // * and, optionally, the matrices T and Z from the Schur decomposition
335: // * H = Z T Z**T, where T is an upper quasi-triangular matrix (the
336: // * Schur form), and Z is the orthogonal matrix of Schur vectors.
337: // *
338: // * Optionally Z may be postmultiplied into an input orthogonal
339: // * matrix Q so that this routine can give the Schur factorization
340: // * of a matrix A which has been reduced to the Hessenberg form H
341: // * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
342: // *
343: // * Arguments
344: // * =========
345: // *
346: // * WANTT (input) LOGICAL
347: // * = .TRUE. : the full Schur form T is required;
348: // * = .FALSE.: only eigenvalues are required.
349: // *
350: // * WANTZ (input) LOGICAL
351: // * = .TRUE. : the matrix of Schur vectors Z is required;
352: // * = .FALSE.: Schur vectors are not required.
353: // *
354: // * N (input) INTEGER
355: // * The order of the matrix H. N .GE. 0.
356: // *
357: // * ILO (input) INTEGER
358: // * IHI (input) INTEGER
359: // * It is assumed that H is already upper triangular in rows
360: // * and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
361: // * H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
362: // * previous call to DGEBAL, and then passed to DGEHRD when the
363: // * matrix output by DGEBAL is reduced to Hessenberg form.
364: // * Otherwise, ILO and IHI should be set to 1 and N,
365: // * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
366: // * If N = 0, then ILO = 1 and IHI = 0.
367: // *
368: // * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
369: // * On entry, the upper Hessenberg matrix H.
370: // * On exit, if INFO = 0 and WANTT is .TRUE., then H contains
371: // * the upper quasi-triangular matrix T from the Schur
372: // * decomposition (the Schur form); 2-by-2 diagonal blocks
373: // * (corresponding to complex conjugate pairs of eigenvalues)
374: // * are returned in standard form, with H(i,i) = H(i+1,i+1)
375: // * and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
376: // * .FALSE., then the contents of H are unspecified on exit.
377: // * (The output value of H when INFO.GT.0 is given under the
378: // * description of INFO below.)
379: // *
380: // * This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
381: // * j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
382: // *
383: // * LDH (input) INTEGER
384: // * The leading dimension of the array H. LDH .GE. max(1,N).
385: // *
386: // * WR (output) DOUBLE PRECISION array, dimension (IHI)
387: // * WI (output) DOUBLE PRECISION array, dimension (IHI)
388: // * The real and imaginary parts, respectively, of the computed
389: // * eigenvalues of H(ILO:IHI,ILO:IHI) are stored WR(ILO:IHI)
390: // * and WI(ILO:IHI). If two eigenvalues are computed as a
391: // * complex conjugate pair, they are stored in consecutive
392: // * elements of WR and WI, say the i-th and (i+1)th, with
393: // * WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
394: // * the eigenvalues are stored in the same order as on the
395: // * diagonal of the Schur form returned in H, with
396: // * WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
397: // * block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
398: // * WI(i+1) = -WI(i).
399: // *
400: // * ILOZ (input) INTEGER
401: // * IHIZ (input) INTEGER
402: // * Specify the rows of Z to which transformations must be
403: // * applied if WANTZ is .TRUE..
404: // * 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
405: // *
406: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
407: // * If WANTZ is .FALSE., then Z is not referenced.
408: // * If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
409: // * replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
410: // * orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
411: // * (The output value of Z when INFO.GT.0 is given under
412: // * the description of INFO below.)
413: // *
414: // * LDZ (input) INTEGER
415: // * The leading dimension of the array Z. if WANTZ is .TRUE.
416: // * then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
417: // *
418: // * WORK (workspace/output) DOUBLE PRECISION array, dimension LWORK
419: // * On exit, if LWORK = -1, WORK(1) returns an estimate of
420: // * the optimal value for LWORK.
421: // *
422: // * LWORK (input) INTEGER
423: // * The dimension of the array WORK. LWORK .GE. max(1,N)
424: // * is sufficient, but LWORK typically as large as 6*N may
425: // * be required for optimal performance. A workspace query
426: // * to determine the optimal workspace size is recommended.
427: // *
428: // * If LWORK = -1, then DLAQR4 does a workspace query.
429: // * In this case, DLAQR4 checks the input parameters and
430: // * estimates the optimal workspace size for the given
431: // * values of N, ILO and IHI. The estimate is returned
432: // * in WORK(1). No error message related to LWORK is
433: // * issued by XERBLA. Neither H nor Z are accessed.
434: // *
435: // *
436: // * INFO (output) INTEGER
437: // * = 0: successful exit
438: // * .GT. 0: if INFO = i, DLAQR4 failed to compute all of
439: // * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
440: // * and WI contain those eigenvalues which have been
441: // * successfully computed. (Failures are rare.)
442: // *
443: // * If INFO .GT. 0 and WANT is .FALSE., then on exit,
444: // * the remaining unconverged eigenvalues are the eigen-
445: // * values of the upper Hessenberg matrix rows and
446: // * columns ILO through INFO of the final, output
447: // * value of H.
448: // *
449: // * If INFO .GT. 0 and WANTT is .TRUE., then on exit
450: // *
451: // * (*) (initial value of H)*U = U*(final value of H)
452: // *
453: // * where U is an orthogonal matrix. The final
454: // * value of H is upper Hessenberg and quasi-triangular
455: // * in rows and columns INFO+1 through IHI.
456: // *
457: // * If INFO .GT. 0 and WANTZ is .TRUE., then on exit
458: // *
459: // * (final value of Z(ILO:IHI,ILOZ:IHIZ)
460: // * = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
461: // *
462: // * where U is the orthogonal matrix in (*) (regard-
463: // * less of the value of WANTT.)
464: // *
465: // * If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
466: // * accessed.
467: // *
468: // * ================================================================
469: // * Based on contributions by
470: // * Karen Braman and Ralph Byers, Department of Mathematics,
471: // * University of Kansas, USA
472: // *
473: // * ================================================================
474: // * References:
475: // * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
476: // * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
477: // * Performance, SIAM Journal of Matrix Analysis, volume 23, pages
478: // * 929--947, 2002.
479: // *
480: // * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
481: // * Algorithm Part II: Aggressive Early Deflation, SIAM Journal
482: // * of Matrix Analysis, volume 23, pages 948--973, 2002.
483: // *
484: // * ================================================================
485: // * .. Parameters ..
486: // *
487: // * ==== Matrices of order NTINY or smaller must be processed by
488: // * . DLAHQR because of insufficient subdiagonal scratch space.
489: // * . (This is a hard limit.) ====
490: // *
491: // * ==== Exceptional deflation windows: try to cure rare
492: // * . slow convergence by increasing the size of the
493: // * . deflation window after KEXNW iterations. =====
494: // *
495: // * ==== Exceptional shifts: try to cure rare slow convergence
496: // * . with ad-hoc exceptional shifts every KEXSH iterations.
497: // * . The constants WILK1 and WILK2 are used to form the
498: // * . exceptional shifts. ====
499: // *
500: // * ..
501: // * .. Local Scalars ..
502: // * ..
503: // * .. External Functions ..
504: // * ..
505: // * .. Local Arrays ..
506: // * ..
507: // * .. External Subroutines ..
508: // * ..
509: // * .. Intrinsic Functions ..
510: // INTRINSIC ABS, DBLE, INT, MAX, MIN, MOD;
511: // * ..
512: // * .. Executable Statements ..
513:
514: #endregion
515:
516:
517: #region Body
518:
519: INFO = 0;
520: // *
521: // * ==== Quick return for N = 0: nothing to do. ====
522: // *
523: if (N == 0)
524: {
525: WORK[1 + o_work] = ONE;
526: return;
527: }
528: // *
529: // * ==== Set up job flags for ILAENV. ====
530: // *
531: if (WANTT)
532: {
533: FortranLib.Copy(ref JBCMPZ, 1, 1, "S");
534: }
535: else
536: {
537: FortranLib.Copy(ref JBCMPZ, 1, 1, "E");
538: }
539: if (WANTZ)
540: {
541: FortranLib.Copy(ref JBCMPZ, 2, 2, "V");
542: }
543: else
544: {
545: FortranLib.Copy(ref JBCMPZ, 2, 2, "N");
546: }
547: // *
548: // * ==== Tiny matrices must use DLAHQR. ====
549: // *
550: if (N <= NTINY)
551: {
552: // *
553: // * ==== Estimate optimal workspace. ====
554: // *
555: LWKOPT = 1;
556: if (LWORK != - 1)
557: {
558: this._dlahqr.Run(WANTT, WANTZ, N, ILO, IHI, ref H, offset_h
559: , LDH, ref WR, offset_wr, ref WI, offset_wi, ILOZ, IHIZ, ref Z, offset_z
560: , LDZ, ref INFO);
561: }
562: }
563: else
564: {
565: // *
566: // * ==== Use small bulge multi-shift QR with aggressive early
567: // * . deflation on larger-than-tiny matrices. ====
568: // *
569: // * ==== Hope for the best. ====
570: // *
571: INFO = 0;
572: // *
573: // * ==== NWR = recommended deflation window size. At this
574: // * . point, N .GT. NTINY = 11, so there is enough
575: // * . subdiagonal workspace for NWR.GE.2 as required.
576: // * . (In fact, there is enough subdiagonal space for
577: // * . NWR.GE.3.) ====
578: // *
579: NWR = this._ilaenv.Run(13, "DLAQR4", JBCMPZ, N, ILO, IHI, LWORK);
580: NWR = Math.Max(2, NWR);
581: NWR = Math.Min(IHI - ILO + 1, Math.Min((N - 1) / 3, NWR));
582: NW = NWR;
583: // *
584: // * ==== NSR = recommended number of simultaneous shifts.
585: // * . At this point N .GT. NTINY = 11, so there is at
586: // * . enough subdiagonal workspace for NSR to be even
587: // * . and greater than or equal to two as required. ====
588: // *
589: NSR = this._ilaenv.Run(15, "DLAQR4", JBCMPZ, N, ILO, IHI, LWORK);
590: NSR = Math.Min(NSR, Math.Min((N + 6) / 9, IHI - ILO));
591: NSR = Math.Max(2, NSR - FortranLib.Mod(NSR,2));
592: // *
593: // * ==== Estimate optimal workspace ====
594: // *
595: // * ==== Workspace query call to DLAQR2 ====
596: // *
597: this._dlaqr2.Run(WANTT, WANTZ, N, ILO, IHI, NWR + 1
598: , ref H, offset_h, LDH, ILOZ, IHIZ, ref Z, offset_z, LDZ
599: , ref LS, ref LD, ref WR, offset_wr, ref WI, offset_wi, ref H, offset_h, LDH
600: , N, ref H, offset_h, LDH, N, ref H, offset_h, LDH
601: , ref WORK, offset_work, - 1);
602: // *
603: // * ==== Optimal workspace = MAX(DLAQR5, DLAQR2) ====
604: // *
605: LWKOPT = Math.Max(3 * NSR / 2, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
606: // *
607: // * ==== Quick return in case of workspace query. ====
608: // *
609: if (LWORK == - 1)
610: {
611: WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
612: return;
613: }
614: // *
615: // * ==== DLAHQR/DLAQR0 crossover point ====
616: // *
617: NMIN = this._ilaenv.Run(12, "DLAQR4", JBCMPZ, N, ILO, IHI, LWORK);
618: NMIN = Math.Max(NTINY, NMIN);
619: // *
620: // * ==== Nibble crossover point ====
621: // *
622: NIBBLE = this._ilaenv.Run(14, "DLAQR4", JBCMPZ, N, ILO, IHI, LWORK);
623: NIBBLE = Math.Max(0, NIBBLE);
624: // *
625: // * ==== Accumulate reflections during ttswp? Use block
626: // * . 2-by-2 structure during matrix-matrix multiply? ====
627: // *
628: KACC22 = this._ilaenv.Run(16, "DLAQR4", JBCMPZ, N, ILO, IHI, LWORK);
629: KACC22 = Math.Max(0, KACC22);
630: KACC22 = Math.Min(2, KACC22);
631: // *
632: // * ==== NWMAX = the largest possible deflation window for
633: // * . which there is sufficient workspace. ====
634: // *
635: NWMAX = Math.Min((N - 1) / 3, LWORK / 2);
636: // *
637: // * ==== NSMAX = the Largest number of simultaneous shifts
638: // * . for which there is sufficient workspace. ====
639: // *
640: NSMAX = Math.Min((N + 6) / 9, 2 * LWORK / 3);
641: NSMAX = NSMAX - FortranLib.Mod(NSMAX,2);
642: // *
643: // * ==== NDFL: an iteration count restarted at deflation. ====
644: // *
645: NDFL = 1;
646: // *
647: // * ==== ITMAX = iteration limit ====
648: // *
649: ITMAX = Math.Max(30, 2 * KEXSH) * Math.Max(10, (IHI - ILO + 1));
650: // *
651: // * ==== Last row and column in the active block ====
652: // *
653: KBOT = IHI;
654: // *
655: // * ==== Main Loop ====
656: // *
657: for (IT = 1; IT <= ITMAX; IT++)
658: {
659: // *
660: // * ==== Done when KBOT falls below ILO ====
661: // *
662: if (KBOT < ILO) goto LABEL90;
663: // *
664: // * ==== Locate active block ====
665: // *
666: for (K = KBOT; K >= ILO + 1; K += - 1)
667: {
668: if (H[K+(K - 1) * LDH + o_h] == ZERO) goto LABEL20;
669: }
670: K = ILO;
671: LABEL20:;
672: KTOP = K;
673: // *
674: // * ==== Select deflation window size ====
675: // *
676: NH = KBOT - KTOP + 1;
677: if (NDFL < KEXNW || NH < NW)
678: {
679: // *
680: // * ==== Typical deflation window. If possible and
681: // * . advisable, nibble the entire active block.
682: // * . If not, use size NWR or NWR+1 depending upon
683: // * . which has the smaller corresponding subdiagonal
684: // * . entry (a heuristic). ====
685: // *
686: NWINC = true;
687: if (NH <= Math.Min(NMIN, NWMAX))
688: {
689: NW = NH;
690: }
691: else
692: {
693: NW = Math.Min(NWR, Math.Min(NH, NWMAX));
694: if (NW < NWMAX)
695: {
696: if (NW >= NH - 1)
697: {
698: NW = NH;
699: }
700: else
701: {
702: KWTOP = KBOT - NW + 1;
703: if (Math.Abs(H[KWTOP+(KWTOP - 1) * LDH + o_h]) > Math.Abs(H[KWTOP - 1+(KWTOP - 2) * LDH + o_h])) NW = NW + 1;
704: }
705: }
706: }
707: }
708: else
709: {
710: // *
711: // * ==== Exceptional deflation window. If there have
712: // * . been no deflations in KEXNW or more iterations,
713: // * . then vary the deflation window size. At first,
714: // * . because, larger windows are, in general, more
715: // * . powerful than smaller ones, rapidly increase the
716: // * . window up to the maximum reasonable and possible.
717: // * . Then maybe try a slightly smaller window. ====
718: // *
719: if (NWINC && NW < Math.Min(NWMAX, NH))
720: {
721: NW = Math.Min(NWMAX, Math.Min(NH, 2 * NW));
722: }
723: else
724: {
725: NWINC = false;
726: if (NW == NH && NH > 2) NW = NH - 1;
727: }
728: }
729: // *
730: // * ==== Aggressive early deflation:
731: // * . split workspace under the subdiagonal into
732: // * . - an nw-by-nw work array V in the lower
733: // * . left-hand-corner,
734: // * . - an NW-by-at-least-NW-but-more-is-better
735: // * . (NW-by-NHO) horizontal work array along
736: // * . the bottom edge,
737: // * . - an at-least-NW-but-more-is-better (NHV-by-NW)
738: // * . vertical work array along the left-hand-edge.
739: // * . ====
740: // *
741: KV = N - NW + 1;
742: KT = NW + 1;
743: NHO = (N - NW - 1) - KT + 1;
744: KWV = NW + 2;
745: NVE = (N - NW) - KWV + 1;
746: // *
747: // * ==== Aggressive early deflation ====
748: // *
749: this._dlaqr2.Run(WANTT, WANTZ, N, KTOP, KBOT, NW
750: , ref H, offset_h, LDH, ILOZ, IHIZ, ref Z, offset_z, LDZ
751: , ref LS, ref LD, ref WR, offset_wr, ref WI, offset_wi, ref H, KV+1 * LDH + o_h, LDH
752: , NHO, ref H, KV+KT * LDH + o_h, LDH, NVE, ref H, KWV+1 * LDH + o_h, LDH
753: , ref WORK, offset_work, LWORK);
754: // *
755: // * ==== Adjust KBOT accounting for new deflations. ====
756: // *
757: KBOT = KBOT - LD;
758: // *
759: // * ==== KS points to the shifts. ====
760: // *
761: KS = KBOT - LS + 1;
762: // *
763: // * ==== Skip an expensive QR sweep if there is a (partly
764: // * . heuristic) reason to expect that many eigenvalues
765: // * . will deflate without it. Here, the QR sweep is
766: // * . skipped if many eigenvalues have just been deflated
767: // * . or if the remaining active block is small.
768: // *
769: if ((LD == 0) || ((100 * LD <= NW * NIBBLE) && (KBOT - KTOP + 1 > Math.Min(NMIN, NWMAX))))
770: {
771: // *
772: // * ==== NS = nominal number of simultaneous shifts.
773: // * . This may be lowered (slightly) if DLAQR2
774: // * . did not provide that many shifts. ====
775: // *
776: NS = Math.Min(NSMAX, Math.Min(NSR, Math.Max(2, KBOT - KTOP)));
777: NS = NS - FortranLib.Mod(NS,2);
778: // *
779: // * ==== If there have been no deflations
780: // * . in a multiple of KEXSH iterations,
781: // * . then try exceptional shifts.
782: // * . Otherwise use shifts provided by
783: // * . DLAQR2 above or from the eigenvalues
784: // * . of a trailing principal submatrix. ====
785: // *
786: if (FortranLib.Mod(NDFL,KEXSH) == 0)
787: {
788: KS = KBOT - NS + 1;
789: for (I = KBOT; I >= Math.Max(KS + 1, KTOP + 2); I += - 2)
790: {
791: SS = Math.Abs(H[I+(I - 1) * LDH + o_h]) + Math.Abs(H[I - 1+(I - 2) * LDH + o_h]);
792: AA = WILK1 * SS + H[I+I * LDH + o_h];
793: BB = SS;
794: CC = WILK2 * SS;
795: DD = AA;
796: this._dlanv2.Run(ref AA, ref BB, ref CC, ref DD, ref WR[I - 1 + o_wr], ref WI[I - 1 + o_wi]
797: , ref WR[I + o_wr], ref WI[I + o_wi], ref CS, ref SN);
798: }
799: if (KS == KTOP)
800: {
801: WR[KS + 1 + o_wr] = H[KS + 1+(KS + 1) * LDH + o_h];
802: WI[KS + 1 + o_wi] = ZERO;
803: WR[KS + o_wr] = WR[KS + 1 + o_wr];
804: WI[KS + o_wi] = WI[KS + 1 + o_wi];
805: }
806: }
807: else
808: {
809: // *
810: // * ==== Got NS/2 or fewer shifts? Use DLAHQR
811: // * . on a trailing principal submatrix to
812: // * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
813: // * . there is enough space below the subdiagonal
814: // * . to fit an NS-by-NS scratch array.) ====
815: // *
816: if (KBOT - KS + 1 <= NS / 2)
817: {
818: KS = KBOT - NS + 1;
819: KT = N - NS + 1;
820: this._dlacpy.Run("A", NS, NS, H, KS+KS * LDH + o_h, LDH, ref H, KT+1 * LDH + o_h
821: , LDH);
822: this._dlahqr.Run(false, false, NS, 1, NS, ref H, KT+1 * LDH + o_h
823: , LDH, ref WR, KS + o_wr, ref WI, KS + o_wi, 1, 1, ref ZDUM, offset_zdum
824: , 1, ref INF);
825: KS = KS + INF;
826: // *
827: // * ==== In case of a rare QR failure use
828: // * . eigenvalues of the trailing 2-by-2
829: // * . principal submatrix. ====
830: // *
831: if (KS >= KBOT)
832: {
833: AA = H[KBOT - 1+(KBOT - 1) * LDH + o_h];
834: CC = H[KBOT+(KBOT - 1) * LDH + o_h];
835: BB = H[KBOT - 1+KBOT * LDH + o_h];
836: DD = H[KBOT+KBOT * LDH + o_h];
837: this._dlanv2.Run(ref AA, ref BB, ref CC, ref DD, ref WR[KBOT - 1 + o_wr], ref WI[KBOT - 1 + o_wi]
838: , ref WR[KBOT + o_wr], ref WI[KBOT + o_wi], ref CS, ref SN);
839: KS = KBOT - 1;
840: }
841: }
842: // *
843: if (KBOT - KS + 1 > NS)
844: {
845: // *
846: // * ==== Sort the shifts (Helps a little)
847: // * . Bubble sort keeps complex conjugate
848: // * . pairs together. ====
849: // *
850: SORTED = false;
851: for (K = KBOT; K >= KS + 1; K += - 1)
852: {
853: if (SORTED) goto LABEL60;
854: SORTED = true;
855: for (I = KS; I <= K - 1; I++)
856: {
857: if (Math.Abs(WR[I + o_wr]) + Math.Abs(WI[I + o_wi]) < Math.Abs(WR[I + 1 + o_wr]) + Math.Abs(WI[I + 1 + o_wi]))
858: {
859: SORTED = false;
860: // *
861: SWAP = WR[I + o_wr];
862: WR[I + o_wr] = WR[I + 1 + o_wr];
863: WR[I + 1 + o_wr] = SWAP;
864: // *
865: SWAP = WI[I + o_wi];
866: WI[I + o_wi] = WI[I + 1 + o_wi];
867: WI[I + 1 + o_wi] = SWAP;
868: }
869: }
870: }
871: LABEL60:;
872: }
873: // *
874: // * ==== Shuffle shifts into pairs of real shifts
875: // * . and pairs of complex conjugate shifts
876: // * . assuming complex conjugate shifts are
877: // * . already adjacent to one another. (Yes,
878: // * . they are.) ====
879: // *
880: for (I = KBOT; I >= KS + 2; I += - 2)
881: {
882: if (WI[I + o_wi] != - WI[I - 1 + o_wi])
883: {
884: // *
885: SWAP = WR[I + o_wr];
886: WR[I + o_wr] = WR[I - 1 + o_wr];
887: WR[I - 1 + o_wr] = WR[I - 2 + o_wr];
888: WR[I - 2 + o_wr] = SWAP;
889: // *
890: SWAP = WI[I + o_wi];
891: WI[I + o_wi] = WI[I - 1 + o_wi];
892: WI[I - 1 + o_wi] = WI[I - 2 + o_wi];
893: WI[I - 2 + o_wi] = SWAP;
894: }
895: }
896: }
897: // *
898: // * ==== If there are only two shifts and both are
899: // * . real, then use only one. ====
900: // *
901: if (KBOT - KS + 1 == 2)
902: {
903: if (WI[KBOT + o_wi] == ZERO)
904: {
905: if (Math.Abs(WR[KBOT + o_wr] - H[KBOT+KBOT * LDH + o_h]) < Math.Abs(WR[KBOT - 1 + o_wr] - H[KBOT+KBOT * LDH + o_h]))
906: {
907: WR[KBOT - 1 + o_wr] = WR[KBOT + o_wr];
908: }
909: else
910: {
911: WR[KBOT + o_wr] = WR[KBOT - 1 + o_wr];
912: }
913: }
914: }
915: // *
916: // * ==== Use up to NS of the the smallest magnatiude
917: // * . shifts. If there aren't NS shifts available,
918: // * . then use them all, possibly dropping one to
919: // * . make the number of shifts even. ====
920: // *
921: NS = Math.Min(NS, KBOT - KS + 1);
922: NS = NS - FortranLib.Mod(NS,2);
923: KS = KBOT - NS + 1;
924: // *
925: // * ==== Small-bulge multi-shift QR sweep:
926: // * . split workspace under the subdiagonal into
927: // * . - a KDU-by-KDU work array U in the lower
928: // * . left-hand-corner,
929: // * . - a KDU-by-at-least-KDU-but-more-is-better
930: // * . (KDU-by-NHo) horizontal work array WH along
931: // * . the bottom edge,
932: // * . - and an at-least-KDU-but-more-is-better-by-KDU
933: // * . (NVE-by-KDU) vertical work WV arrow along
934: // * . the left-hand-edge. ====
935: // *
936: KDU = 3 * NS - 3;
937: KU = N - KDU + 1;
938: KWH = KDU + 1;
939: NHO = (N - KDU + 1 - 4) - (KDU + 1) + 1;
940: KWV = KDU + 4;
941: NVE = N - KDU - KWV + 1;
942: // *
943: // * ==== Small-bulge multi-shift QR sweep ====
944: // *
945: this._dlaqr5.Run(WANTT, WANTZ, KACC22, N, KTOP, KBOT
946: , NS, ref WR, KS + o_wr, ref WI, KS + o_wi, ref H, offset_h, LDH, ILOZ
947: , IHIZ, ref Z, offset_z, LDZ, ref WORK, offset_work, 3, ref H, KU+1 * LDH + o_h
948: , LDH, NVE, ref H, KWV+1 * LDH + o_h, LDH, NHO, ref H, KU+KWH * LDH + o_h
949: , LDH);
950: }
951: // *
952: // * ==== Note progress (or the lack of it). ====
953: // *
954: if (LD > 0)
955: {
956: NDFL = 1;
957: }
958: else
959: {
960: NDFL = NDFL + 1;
961: }
962: // *
963: // * ==== End of main loop ====
964: }
965: // *
966: // * ==== Iteration limit exceeded. Set INFO to show where
967: // * . the problem occurred and exit. ====
968: // *
969: INFO = KBOT;
970: LABEL90:;
971: }
972: // *
973: // * ==== Return the optimal value of LWORK. ====
974: // *
975: WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
976: // *
977: // * ==== End of DLAQR4 ====
978: // *
979:
980: #endregion
981:
982: }
983: }
984: }