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