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 driver routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DHSEQR 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 DHSEQR
38: {
39:
40:
41: #region Dependencies
42:
43: ILAENV _ilaenv; LSAME _lsame; DLACPY _dlacpy; DLAHQR _dlahqr; DLAQR0 _dlaqr0; DLASET _dlaset; XERBLA _xerbla;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const int NTINY = 11; const int NL = 49; const double ZERO = 0.0E0; const double ONE = 1.0E0;
51: double[] HL = new double[NL * NL]; int offset_hl = 0; int o_hl = -1 - NL;
52: double[] WORKL = new double[NL]; int offset_workl = 0; int I = 0; int KBOT = 0; int NMIN = 0;
53: bool INITZ = false;bool LQUERY = false; bool WANTT = false; bool WANTZ = false;
54:
55: #endregion
56:
57: public DHSEQR(ILAENV ilaenv, LSAME lsame, DLACPY dlacpy, DLAHQR dlahqr, DLAQR0 dlaqr0, DLASET dlaset, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._ilaenv = ilaenv; this._lsame = lsame; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlaqr0 = dlaqr0;
64: this._dlaset = dlaset;this._xerbla = xerbla;
65:
66: #endregion
67:
68: }
69:
70: public DHSEQR()
71: {
72:
73:
74: #region Dependencies (Initialization)
75:
76: IEEECK ieeeck = new IEEECK();
77: IPARMQ iparmq = new IPARMQ();
78: LSAME lsame = new LSAME();
79: DLAMC3 dlamc3 = new DLAMC3();
80: DCOPY dcopy = new DCOPY();
81: DLABAD dlabad = new DLABAD();
82: DLAPY2 dlapy2 = new DLAPY2();
83: DNRM2 dnrm2 = new DNRM2();
84: DSCAL dscal = new DSCAL();
85: DROT drot = new DROT();
86: DAXPY daxpy = new DAXPY();
87: XERBLA xerbla = new XERBLA();
88: DLASSQ dlassq = new DLASSQ();
89: IDAMAX idamax = new IDAMAX();
90: DSWAP dswap = new DSWAP();
91: DLAQR1 dlaqr1 = new DLAQR1();
92: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
93: DLACPY dlacpy = new DLACPY(lsame);
94: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
95: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
96: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
97: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
98: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
99: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
100: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
101: DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
102: DGEMV dgemv = new DGEMV(lsame, xerbla);
103: DGER dger = new DGER(xerbla);
104: DLARF dlarf = new DLARF(dgemv, dger, lsame);
105: DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
106: DGEMM dgemm = new DGEMM(lsame, xerbla);
107: DTRMM dtrmm = new DTRMM(lsame, xerbla);
108: DTRMV dtrmv = new DTRMV(lsame, xerbla);
109: DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
110: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
111: DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
112: DLASET dlaset = new DLASET(lsame);
113: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
114: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
115: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
116: DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
117: DLANGE dlange = new DLANGE(dlassq, lsame);
118: DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
119: DLARTG dlartg = new DLARTG(dlamch);
120: DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
121: DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
122: DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
123: DLAQR2 dlaqr2 = new DLAQR2(dlamch, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlarf, dlarfg
124: , dlaset, dorghr, dtrexc);
125: DLAQR5 dlaqr5 = new DLAQR5(dlamch, dgemm, dlabad, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm);
126: DLAQR4 dlaqr4 = new DLAQR4(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5);
127: DLAQR3 dlaqr3 = new DLAQR3(dlamch, ilaenv, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlaqr4
128: , dlarf, dlarfg, dlaset, dorghr, dtrexc);
129: DLAQR0 dlaqr0 = new DLAQR0(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr3, dlaqr4, dlaqr5);
130:
131: #endregion
132:
133:
134: #region Set Dependencies
135:
136: this._ilaenv = ilaenv; this._lsame = lsame; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlaqr0 = dlaqr0;
137: this._dlaset = dlaset;this._xerbla = xerbla;
138:
139: #endregion
140:
141: }
142: /// <summary>
143: /// Purpose
144: /// =======
145: ///
146: /// DHSEQR 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="JOB">
158: /// (input) CHARACTER*1
159: /// = 'E': compute eigenvalues only;
160: /// = 'S': compute eigenvalues and the Schur form T.
161: ///</param>
162: /// <param name="COMPZ">
163: /// (input) CHARACTER*1
164: /// = 'N': no Schur vectors are computed;
165: /// = 'I': Z is initialized to the unit matrix and the matrix Z
166: /// of Schur vectors of H is returned;
167: /// = 'V': Z must contain an orthogonal matrix Q on entry, and
168: /// the product Q*Z is returned.
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. ILO and IHI are normally
181: /// set by a previous call to DGEBAL, and then passed to DGEHRD
182: /// when the matrix output by DGEBAL is reduced to Hessenberg
183: /// form. Otherwise ILO and IHI should be set to 1 and N
184: /// respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
185: /// If N = 0, then ILO = 1 and IHI = 0.
186: ///</param>
187: /// <param name="H">
188: /// (input/output) DOUBLE PRECISION array, dimension (LDH,N)
189: /// On entry, the upper Hessenberg matrix H.
190: /// On exit, if INFO = 0 and JOB = 'S', then H contains the
191: /// upper quasi-triangular matrix T from the Schur decomposition
192: /// (the Schur form); 2-by-2 diagonal blocks (corresponding to
193: /// complex conjugate pairs of eigenvalues) are returned in
194: /// standard form, with H(i,i) = H(i+1,i+1) and
195: /// H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
196: /// contents of H are unspecified on exit. (The output value of
197: /// H when INFO.GT.0 is given under the description of INFO
198: /// below.)
199: ///
200: /// Unlike earlier versions of DHSEQR, this subroutine may
201: /// explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
202: /// 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 (N)
210: ///</param>
211: /// <param name="WI">
212: /// (output) DOUBLE PRECISION array, dimension (N)
213: /// The real and imaginary parts, respectively, of the computed
214: /// eigenvalues. If two eigenvalues are computed as a complex
215: /// conjugate pair, they are stored in consecutive elements of
216: /// WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
217: /// WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
218: /// the same order as on the diagonal of the Schur form returned
219: /// in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
220: /// diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
221: /// WI(i+1) = -WI(i).
222: ///</param>
223: /// <param name="Z">
224: /// (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
225: /// If COMPZ = 'N', Z is not referenced.
226: /// If COMPZ = 'I', on entry Z need not be set and on exit,
227: /// if INFO = 0, Z contains the orthogonal matrix Z of the Schur
228: /// vectors of H. If COMPZ = 'V', on entry Z must contain an
229: /// N-by-N matrix Q, which is assumed to be equal to the unit
230: /// matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
231: /// if INFO = 0, Z contains Q*Z.
232: /// Normally Q is the orthogonal matrix generated by DORGHR
233: /// after the call to DGEHRD which formed the Hessenberg matrix
234: /// H. (The output value of Z when INFO.GT.0 is given under
235: /// the description of INFO below.)
236: ///</param>
237: /// <param name="LDZ">
238: /// (input) INTEGER
239: /// The leading dimension of the array Z. if COMPZ = 'I' or
240: /// COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
241: ///</param>
242: /// <param name="WORK">
243: /// (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
244: /// On exit, if INFO = 0, WORK(1) returns an estimate of
245: /// the optimal value for LWORK.
246: ///</param>
247: /// <param name="LWORK">
248: /// (input) INTEGER
249: /// The dimension of the array WORK. LWORK .GE. max(1,N)
250: /// is sufficient, but LWORK typically as large as 6*N may
251: /// be required for optimal performance. A workspace query
252: /// to determine the optimal workspace size is recommended.
253: ///
254: /// If LWORK = -1, then DHSEQR does a workspace query.
255: /// In this case, DHSEQR checks the input parameters and
256: /// estimates the optimal workspace size for the given
257: /// values of N, ILO and IHI. The estimate is returned
258: /// in WORK(1). No error message related to LWORK is
259: /// issued by XERBLA. Neither H nor Z are accessed.
260: ///
261: ///</param>
262: /// <param name="INFO">
263: /// (output) INTEGER
264: /// = 0: successful exit
265: /// .LT. 0: if INFO = -i, the i-th argument had an illegal
266: /// value
267: /// .GT. 0: if INFO = i, DHSEQR failed to compute all of
268: /// the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
269: /// and WI contain those eigenvalues which have been
270: /// successfully computed. (Failures are rare.)
271: ///
272: /// If INFO .GT. 0 and JOB = 'E', then on exit, the
273: /// remaining unconverged eigenvalues are the eigen-
274: /// values of the upper Hessenberg matrix rows and
275: /// columns ILO through INFO of the final, output
276: /// value of H.
277: ///
278: /// If INFO .GT. 0 and JOB = 'S', then on exit
279: ///
280: /// (*) (initial value of H)*U = U*(final value of H)
281: ///
282: /// where U is an orthogonal matrix. The final
283: /// value of H is upper Hessenberg and quasi-triangular
284: /// in rows and columns INFO+1 through IHI.
285: ///
286: /// If INFO .GT. 0 and COMPZ = 'V', then on exit
287: ///
288: /// (final value of Z) = (initial value of Z)*U
289: ///
290: /// where U is the orthogonal matrix in (*) (regard-
291: /// less of the value of JOB.)
292: ///
293: /// If INFO .GT. 0 and COMPZ = 'I', then on exit
294: /// (final value of Z) = U
295: /// where U is the orthogonal matrix in (*) (regard-
296: /// less of the value of JOB.)
297: ///
298: /// If INFO .GT. 0 and COMPZ = 'N', then Z is not
299: /// accessed.
300: ///</param>
301: public void Run(string JOB, string COMPZ, int N, int ILO, int IHI, ref double[] H, int offset_h
302: , int LDH, ref double[] WR, int offset_wr, ref double[] WI, int offset_wi, ref double[] Z, int offset_z, int LDZ, ref double[] WORK, int offset_work
303: , int LWORK, ref int INFO)
304: {
305:
306: #region Array Index Correction
307:
308: int o_h = -1 - LDH + offset_h; int o_wr = -1 + offset_wr; int o_wi = -1 + offset_wi;
309: int o_z = -1 - LDZ + offset_z; int o_work = -1 + offset_work;
310:
311: #endregion
312:
313:
314: #region Strings
315:
316: JOB = JOB.Substring(0, 1); COMPZ = COMPZ.Substring(0, 1);
317:
318: #endregion
319:
320:
321: #region Prolog
322:
323: // *
324: // * -- LAPACK driver routine (version 3.1) --
325: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
326: // * November 2006
327: // *
328: // * .. Scalar Arguments ..
329: // * ..
330: // * .. Array Arguments ..
331: // * ..
332: // * Purpose
333: // * =======
334: // *
335: // * DHSEQR computes the eigenvalues of a Hessenberg matrix H
336: // * and, optionally, the matrices T and Z from the Schur decomposition
337: // * H = Z T Z**T, where T is an upper quasi-triangular matrix (the
338: // * Schur form), and Z is the orthogonal matrix of Schur vectors.
339: // *
340: // * Optionally Z may be postmultiplied into an input orthogonal
341: // * matrix Q so that this routine can give the Schur factorization
342: // * of a matrix A which has been reduced to the Hessenberg form H
343: // * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
344: // *
345: // * Arguments
346: // * =========
347: // *
348: // * JOB (input) CHARACTER*1
349: // * = 'E': compute eigenvalues only;
350: // * = 'S': compute eigenvalues and the Schur form T.
351: // *
352: // * COMPZ (input) CHARACTER*1
353: // * = 'N': no Schur vectors are computed;
354: // * = 'I': Z is initialized to the unit matrix and the matrix Z
355: // * of Schur vectors of H is returned;
356: // * = 'V': Z must contain an orthogonal matrix Q on entry, and
357: // * the product Q*Z is returned.
358: // *
359: // * N (input) INTEGER
360: // * The order of the matrix H. N .GE. 0.
361: // *
362: // * ILO (input) INTEGER
363: // * IHI (input) INTEGER
364: // * It is assumed that H is already upper triangular in rows
365: // * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
366: // * set by a previous call to DGEBAL, and then passed to DGEHRD
367: // * when the matrix output by DGEBAL is reduced to Hessenberg
368: // * form. Otherwise ILO and IHI should be set to 1 and N
369: // * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
370: // * If N = 0, then ILO = 1 and IHI = 0.
371: // *
372: // * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
373: // * On entry, the upper Hessenberg matrix H.
374: // * On exit, if INFO = 0 and JOB = 'S', then H contains the
375: // * upper quasi-triangular matrix T from the Schur decomposition
376: // * (the Schur form); 2-by-2 diagonal blocks (corresponding to
377: // * complex conjugate pairs of eigenvalues) are returned in
378: // * standard form, with H(i,i) = H(i+1,i+1) and
379: // * H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
380: // * contents of H are unspecified on exit. (The output value of
381: // * H when INFO.GT.0 is given under the description of INFO
382: // * below.)
383: // *
384: // * Unlike earlier versions of DHSEQR, this subroutine may
385: // * explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
386: // * or j = IHI+1, IHI+2, ... N.
387: // *
388: // * LDH (input) INTEGER
389: // * The leading dimension of the array H. LDH .GE. max(1,N).
390: // *
391: // * WR (output) DOUBLE PRECISION array, dimension (N)
392: // * WI (output) DOUBLE PRECISION array, dimension (N)
393: // * The real and imaginary parts, respectively, of the computed
394: // * eigenvalues. If two eigenvalues are computed as a complex
395: // * conjugate pair, they are stored in consecutive elements of
396: // * WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
397: // * WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
398: // * the same order as on the diagonal of the Schur form returned
399: // * in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
400: // * diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
401: // * WI(i+1) = -WI(i).
402: // *
403: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
404: // * If COMPZ = 'N', Z is not referenced.
405: // * If COMPZ = 'I', on entry Z need not be set and on exit,
406: // * if INFO = 0, Z contains the orthogonal matrix Z of the Schur
407: // * vectors of H. If COMPZ = 'V', on entry Z must contain an
408: // * N-by-N matrix Q, which is assumed to be equal to the unit
409: // * matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
410: // * if INFO = 0, Z contains Q*Z.
411: // * Normally Q is the orthogonal matrix generated by DORGHR
412: // * after the call to DGEHRD which formed the Hessenberg matrix
413: // * H. (The output value of Z when INFO.GT.0 is given under
414: // * the description of INFO below.)
415: // *
416: // * LDZ (input) INTEGER
417: // * The leading dimension of the array Z. if COMPZ = 'I' or
418: // * COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
419: // *
420: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
421: // * On exit, if INFO = 0, WORK(1) returns an estimate of
422: // * the optimal value for LWORK.
423: // *
424: // * LWORK (input) INTEGER
425: // * The dimension of the array WORK. LWORK .GE. max(1,N)
426: // * is sufficient, but LWORK typically as large as 6*N may
427: // * be required for optimal performance. A workspace query
428: // * to determine the optimal workspace size is recommended.
429: // *
430: // * If LWORK = -1, then DHSEQR does a workspace query.
431: // * In this case, DHSEQR checks the input parameters and
432: // * estimates the optimal workspace size for the given
433: // * values of N, ILO and IHI. The estimate is returned
434: // * in WORK(1). No error message related to LWORK is
435: // * issued by XERBLA. Neither H nor Z are accessed.
436: // *
437: // *
438: // * INFO (output) INTEGER
439: // * = 0: successful exit
440: // * .LT. 0: if INFO = -i, the i-th argument had an illegal
441: // * value
442: // * .GT. 0: if INFO = i, DHSEQR failed to compute all of
443: // * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
444: // * and WI contain those eigenvalues which have been
445: // * successfully computed. (Failures are rare.)
446: // *
447: // * If INFO .GT. 0 and JOB = 'E', then on exit, the
448: // * remaining unconverged eigenvalues are the eigen-
449: // * values of the upper Hessenberg matrix rows and
450: // * columns ILO through INFO of the final, output
451: // * value of H.
452: // *
453: // * If INFO .GT. 0 and JOB = 'S', then on exit
454: // *
455: // * (*) (initial value of H)*U = U*(final value of H)
456: // *
457: // * where U is an orthogonal matrix. The final
458: // * value of H is upper Hessenberg and quasi-triangular
459: // * in rows and columns INFO+1 through IHI.
460: // *
461: // * If INFO .GT. 0 and COMPZ = 'V', then on exit
462: // *
463: // * (final value of Z) = (initial value of Z)*U
464: // *
465: // * where U is the orthogonal matrix in (*) (regard-
466: // * less of the value of JOB.)
467: // *
468: // * If INFO .GT. 0 and COMPZ = 'I', then on exit
469: // * (final value of Z) = U
470: // * where U is the orthogonal matrix in (*) (regard-
471: // * less of the value of JOB.)
472: // *
473: // * If INFO .GT. 0 and COMPZ = 'N', then Z is not
474: // * accessed.
475: // *
476: // * ================================================================
477: // * Default values supplied by
478: // * ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
479: // * It is suggested that these defaults be adjusted in order
480: // * to attain best performance in each particular
481: // * computational environment.
482: // *
483: // * ISPEC=1: The DLAHQR vs DLAQR0 crossover point.
484: // * Default: 75. (Must be at least 11.)
485: // *
486: // * ISPEC=2: Recommended deflation window size.
487: // * This depends on ILO, IHI and NS. NS is the
488: // * number of simultaneous shifts returned
489: // * by ILAENV(ISPEC=4). (See ISPEC=4 below.)
490: // * The default for (IHI-ILO+1).LE.500 is NS.
491: // * The default for (IHI-ILO+1).GT.500 is 3*NS/2.
492: // *
493: // * ISPEC=3: Nibble crossover point. (See ILAENV for
494: // * details.) Default: 14% of deflation window
495: // * size.
496: // *
497: // * ISPEC=4: Number of simultaneous shifts, NS, in
498: // * a multi-shift QR iteration.
499: // *
500: // * If IHI-ILO+1 is ...
501: // *
502: // * greater than ...but less ... the
503: // * or equal to ... than default is
504: // *
505: // * 1 30 NS - 2(+)
506: // * 30 60 NS - 4(+)
507: // * 60 150 NS = 10(+)
508: // * 150 590 NS = **
509: // * 590 3000 NS = 64
510: // * 3000 6000 NS = 128
511: // * 6000 infinity NS = 256
512: // *
513: // * (+) By default some or all matrices of this order
514: // * are passed to the implicit double shift routine
515: // * DLAHQR and NS is ignored. See ISPEC=1 above
516: // * and comments in IPARM for details.
517: // *
518: // * The asterisks (**) indicate an ad-hoc
519: // * function of N increasing from 10 to 64.
520: // *
521: // * ISPEC=5: Select structured matrix multiply.
522: // * (See ILAENV for details.) Default: 3.
523: // *
524: // * ================================================================
525: // * Based on contributions by
526: // * Karen Braman and Ralph Byers, Department of Mathematics,
527: // * University of Kansas, USA
528: // *
529: // * ================================================================
530: // * References:
531: // * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
532: // * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
533: // * Performance, SIAM Journal of Matrix Analysis, volume 23, pages
534: // * 929--947, 2002.
535: // *
536: // * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
537: // * Algorithm Part II: Aggressive Early Deflation, SIAM Journal
538: // * of Matrix Analysis, volume 23, pages 948--973, 2002.
539: // *
540: // * ================================================================
541: // * .. Parameters ..
542: // *
543: // * ==== Matrices of order NTINY or smaller must be processed by
544: // * . DLAHQR because of insufficient subdiagonal scratch space.
545: // * . (This is a hard limit.) ====
546: // *
547: // * ==== NL allocates some local workspace to help small matrices
548: // * . through a rare DLAHQR failure. NL .GT. NTINY = 11 is
549: // * . required and NL .LE. NMIN = ILAENV(ISPEC=1,...) is recom-
550: // * . mended. (The default value of NMIN is 75.) Using NL = 49
551: // * . allows up to six simultaneous shifts and a 16-by-16
552: // * . deflation window. ====
553: // *
554: // * ..
555: // * .. Local Arrays ..
556: // * ..
557: // * .. Local Scalars ..
558: // * ..
559: // * .. External Functions ..
560: // * ..
561: // * .. External Subroutines ..
562: // * ..
563: // * .. Intrinsic Functions ..
564: // INTRINSIC DBLE, MAX, MIN;
565: // * ..
566: // * .. Executable Statements ..
567: // *
568: // * ==== Decode and check the input parameters. ====
569: // *
570:
571: #endregion
572:
573:
574: #region Body
575:
576: WANTT = this._lsame.Run(JOB, "S");
577: INITZ = this._lsame.Run(COMPZ, "I");
578: WANTZ = INITZ || this._lsame.Run(COMPZ, "V");
579: WORK[1 + o_work] = Convert.ToDouble(Math.Max(1, N));
580: LQUERY = LWORK == - 1;
581: // *
582: INFO = 0;
583: if (!this._lsame.Run(JOB, "E") && !WANTT)
584: {
585: INFO = - 1;
586: }
587: else
588: {
589: if (!this._lsame.Run(COMPZ, "N") && !WANTZ)
590: {
591: INFO = - 2;
592: }
593: else
594: {
595: if (N < 0)
596: {
597: INFO = - 3;
598: }
599: else
600: {
601: if (ILO < 1 || ILO > Math.Max(1, N))
602: {
603: INFO = - 4;
604: }
605: else
606: {
607: if (IHI < Math.Min(ILO, N) || IHI > N)
608: {
609: INFO = - 5;
610: }
611: else
612: {
613: if (LDH < Math.Max(1, N))
614: {
615: INFO = - 7;
616: }
617: else
618: {
619: if (LDZ < 1 || (WANTZ && LDZ < Math.Max(1, N)))
620: {
621: INFO = - 11;
622: }
623: else
624: {
625: if (LWORK < Math.Max(1, N) && !LQUERY)
626: {
627: INFO = - 13;
628: }
629: }
630: }
631: }
632: }
633: }
634: }
635: }
636: // *
637: if (INFO != 0)
638: {
639: // *
640: // * ==== Quick return in case of invalid argument. ====
641: // *
642: this._xerbla.Run("DHSEQR", - INFO);
643: return;
644: // *
645: }
646: else
647: {
648: if (N == 0)
649: {
650: // *
651: // * ==== Quick return in case N = 0; nothing to do. ====
652: // *
653: return;
654: // *
655: }
656: else
657: {
658: if (LQUERY)
659: {
660: // *
661: // * ==== Quick return in case of a workspace query ====
662: // *
663: this._dlaqr0.Run(WANTT, WANTZ, N, ILO, IHI, ref H, offset_h
664: , LDH, ref WR, offset_wr, ref WI, offset_wi, ILO, IHI, ref Z, offset_z
665: , LDZ, ref WORK, offset_work, LWORK, ref INFO);
666: // * ==== Ensure reported workspace size is backward-compatible with
667: // * . previous LAPACK versions. ====
668: WORK[1 + o_work] = Math.Max(Convert.ToDouble(Math.Max(1, N)), WORK[1 + o_work]);
669: return;
670: // *
671: }
672: else
673: {
674: // *
675: // * ==== copy eigenvalues isolated by DGEBAL ====
676: // *
677: for (I = 1; I <= ILO - 1; I++)
678: {
679: WR[I + o_wr] = H[I+I * LDH + o_h];
680: WI[I + o_wi] = ZERO;
681: }
682: for (I = IHI + 1; I <= N; I++)
683: {
684: WR[I + o_wr] = H[I+I * LDH + o_h];
685: WI[I + o_wi] = ZERO;
686: }
687: // *
688: // * ==== Initialize Z, if requested ====
689: // *
690: if (INITZ)
691: {
692: this._dlaset.Run("A", N, N, ZERO, ONE, ref Z, offset_z
693: , LDZ);
694: }
695: // *
696: // * ==== Quick return if possible ====
697: // *
698: if (ILO == IHI)
699: {
700: WR[ILO + o_wr] = H[ILO+ILO * LDH + o_h];
701: WI[ILO + o_wi] = ZERO;
702: return;
703: }
704: // *
705: // * ==== DLAHQR/DLAQR0 crossover point ====
706: // *
707: NMIN = this._ilaenv.Run(12, "DHSEQR", FortranLib.Substring(JOB, 1, 1) + FortranLib.Substring(COMPZ, 1, 1), N, ILO, IHI, LWORK);
708: NMIN = Math.Max(NTINY, NMIN);
709: // *
710: // * ==== DLAQR0 for big matrices; DLAHQR for small ones ====
711: // *
712: if (N > NMIN)
713: {
714: this._dlaqr0.Run(WANTT, WANTZ, N, ILO, IHI, ref H, offset_h
715: , LDH, ref WR, offset_wr, ref WI, offset_wi, ILO, IHI, ref Z, offset_z
716: , LDZ, ref WORK, offset_work, LWORK, ref INFO);
717: }
718: else
719: {
720: // *
721: // * ==== Small matrix ====
722: // *
723: this._dlahqr.Run(WANTT, WANTZ, N, ILO, IHI, ref H, offset_h
724: , LDH, ref WR, offset_wr, ref WI, offset_wi, ILO, IHI, ref Z, offset_z
725: , LDZ, ref INFO);
726: // *
727: if (INFO > 0)
728: {
729: // *
730: // * ==== A rare DLAHQR failure! DLAQR0 sometimes succeeds
731: // * . when DLAHQR fails. ====
732: // *
733: KBOT = INFO;
734: // *
735: if (N >= NL)
736: {
737: // *
738: // * ==== Larger matrices have enough subdiagonal scratch
739: // * . space to call DLAQR0 directly. ====
740: // *
741: this._dlaqr0.Run(WANTT, WANTZ, N, ILO, KBOT, ref H, offset_h
742: , LDH, ref WR, offset_wr, ref WI, offset_wi, ILO, IHI, ref Z, offset_z
743: , LDZ, ref WORK, offset_work, LWORK, ref INFO);
744: // *
745: }
746: else
747: {
748: // *
749: // * ==== Tiny matrices don't have enough subdiagonal
750: // * . scratch space to benefit from DLAQR0. Hence,
751: // * . tiny matrices must be copied into a larger
752: // * . array before calling DLAQR0. ====
753: // *
754: this._dlacpy.Run("A", N, N, H, offset_h, LDH, ref HL, offset_hl
755: , NL);
756: HL[N + 1+N * NL + o_hl] = ZERO;
757: this._dlaset.Run("A", NL, NL - N, ZERO, ZERO, ref HL, 1+(N + 1) * NL + o_hl
758: , NL);
759: this._dlaqr0.Run(WANTT, WANTZ, NL, ILO, KBOT, ref HL, offset_hl
760: , NL, ref WR, offset_wr, ref WI, offset_wi, ILO, IHI, ref Z, offset_z
761: , LDZ, ref WORKL, offset_workl, NL, ref INFO);
762: if (WANTT || INFO != 0)
763: {
764: this._dlacpy.Run("A", N, N, HL, offset_hl, NL, ref H, offset_h
765: , LDH);
766: }
767: }
768: }
769: }
770: // *
771: // * ==== Clear out the trash, if necessary. ====
772: // *
773: if ((WANTT || INFO != 0) && N > 2)
774: {
775: this._dlaset.Run("L", N - 2, N - 2, ZERO, ZERO, ref H, 3+1 * LDH + o_h
776: , LDH);
777: }
778: // *
779: // * ==== Ensure reported workspace size is backward-compatible with
780: // * . previous LAPACK versions. ====
781: // *
782: WORK[1 + o_work] = Math.Max(Convert.ToDouble(Math.Max(1, N)), WORK[1 + o_work]);
783: }
784: }
785: }
786: // *
787: // * ==== End of DHSEQR ====
788: // *
789:
790: #endregion
791:
792: }
793: }
794: }