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: ///</summary>
24: public class DLAQR2
25: {
26:
27:
28: #region Dependencies
29:
30: DLAMCH _dlamch; DCOPY _dcopy; DGEHRD _dgehrd; DGEMM _dgemm; DLABAD _dlabad; DLACPY _dlacpy; DLAHQR _dlahqr;
31: DLANV2 _dlanv2;DLARF _dlarf; DLARFG _dlarfg; DLASET _dlaset; DORGHR _dorghr; DTREXC _dtrexc;
32:
33: #endregion
34:
35:
36: #region Fields
37:
38: const double ZERO = 0.0E0; const double ONE = 1.0E0; double AA = 0; double BB = 0; double BETA = 0; double CC = 0;
39: double CS = 0;double DD = 0; double EVI = 0; double EVK = 0; double FOO = 0; double S = 0; double SAFMAX = 0;
40: double SAFMIN = 0;double SMLNUM = 0; double SN = 0; double TAU = 0; double ULP = 0; int I = 0; int IFST = 0; int ILST = 0;
41: int INFO = 0;int INFQR = 0; int J = 0; int JW = 0; int K = 0; int KCOL = 0; int KEND = 0; int KLN = 0; int KROW = 0;
42: int KWTOP = 0;int LTOP = 0; int LWK1 = 0; int LWK2 = 0; int LWKOPT = 0; bool BULGE = false; bool SORTED = false;
43:
44: #endregion
45:
46: public DLAQR2(DLAMCH dlamch, DCOPY dcopy, DGEHRD dgehrd, DGEMM dgemm, DLABAD dlabad, DLACPY dlacpy, DLAHQR dlahqr, DLANV2 dlanv2, DLARF dlarf, DLARFG dlarfg
47: , DLASET dlaset, DORGHR dorghr, DTREXC dtrexc)
48: {
49:
50:
51: #region Set Dependencies
52:
53: this._dlamch = dlamch; this._dcopy = dcopy; this._dgehrd = dgehrd; this._dgemm = dgemm; this._dlabad = dlabad;
54: this._dlacpy = dlacpy;this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlarf = dlarf; this._dlarfg = dlarfg;
55: this._dlaset = dlaset;this._dorghr = dorghr; this._dtrexc = dtrexc;
56:
57: #endregion
58:
59: }
60:
61: public DLAQR2()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: LSAME lsame = new LSAME();
68: DLAMC3 dlamc3 = new DLAMC3();
69: DCOPY dcopy = new DCOPY();
70: DAXPY daxpy = new DAXPY();
71: XERBLA xerbla = new XERBLA();
72: DLAPY2 dlapy2 = new DLAPY2();
73: DNRM2 dnrm2 = new DNRM2();
74: DSCAL dscal = new DSCAL();
75: IEEECK ieeeck = new IEEECK();
76: IPARMQ iparmq = new IPARMQ();
77: DLABAD dlabad = new DLABAD();
78: DROT drot = new DROT();
79: DLASSQ dlassq = new DLASSQ();
80: IDAMAX idamax = new IDAMAX();
81: DSWAP dswap = new DSWAP();
82: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
83: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
84: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
85: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
86: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
87: DGEMV dgemv = new DGEMV(lsame, xerbla);
88: DGER dger = new DGER(xerbla);
89: DLARF dlarf = new DLARF(dgemv, dger, lsame);
90: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
91: DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
92: DGEMM dgemm = new DGEMM(lsame, xerbla);
93: DLACPY dlacpy = new DLACPY(lsame);
94: DTRMM dtrmm = new DTRMM(lsame, xerbla);
95: DTRMV dtrmv = new DTRMV(lsame, xerbla);
96: DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
97: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
98: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
99: DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
100: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
101: DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
102: DLASET dlaset = new DLASET(lsame);
103: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
104: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
105: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
106: DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
107: DLANGE dlange = new DLANGE(dlassq, lsame);
108: DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
109: DLARTG dlartg = new DLARTG(dlamch);
110: DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
111: DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
112: DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
113:
114: #endregion
115:
116:
117: #region Set Dependencies
118:
119: this._dlamch = dlamch; this._dcopy = dcopy; this._dgehrd = dgehrd; this._dgemm = dgemm; this._dlabad = dlabad;
120: this._dlacpy = dlacpy;this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlarf = dlarf; this._dlarfg = dlarfg;
121: this._dlaset = dlaset;this._dorghr = dorghr; this._dtrexc = dtrexc;
122:
123: #endregion
124:
125: }
126: /// <param name="WANTT">
127: /// (input) LOGICAL
128: /// If .TRUE., then the Hessenberg matrix H is fully updated
129: /// so that the quasi-triangular Schur factor may be
130: /// computed (in cooperation with the calling subroutine).
131: /// If .FALSE., then only enough of H is updated to preserve
132: /// the eigenvalues.
133: ///</param>
134: /// <param name="WANTZ">
135: /// (input) LOGICAL
136: /// If .TRUE., then the orthogonal matrix Z is updated so
137: /// so that the orthogonal Schur factor may be computed
138: /// (in cooperation with the calling subroutine).
139: /// If .FALSE., then Z is not referenced.
140: ///</param>
141: /// <param name="N">
142: /// (input) INTEGER
143: /// The order of the matrix H and (if WANTZ is .TRUE.) the
144: /// order of the orthogonal matrix Z.
145: ///</param>
146: /// <param name="KTOP">
147: /// (input) INTEGER
148: /// It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
149: /// KBOT and KTOP together determine an isolated block
150: /// along the diagonal of the Hessenberg matrix.
151: ///</param>
152: /// <param name="KBOT">
153: /// (input) INTEGER
154: /// It is assumed without a check that either
155: /// KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
156: /// determine an isolated block along the diagonal of the
157: /// Hessenberg matrix.
158: ///</param>
159: /// <param name="NW">
160: /// (input) INTEGER
161: /// Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
162: ///</param>
163: /// <param name="H">
164: /// (input/output) DOUBLE PRECISION array, dimension (LDH,N)
165: /// On input the initial N-by-N section of H stores the
166: /// Hessenberg matrix undergoing aggressive early deflation.
167: /// On output H has been transformed by an orthogonal
168: /// similarity transformation, perturbed, and the returned
169: /// to Hessenberg form that (it is to be hoped) has some
170: /// zero subdiagonal entries.
171: ///</param>
172: /// <param name="LDH">
173: /// (input) integer
174: /// Leading dimension of H just as declared in the calling
175: /// subroutine. N .LE. LDH
176: ///</param>
177: /// <param name="ILOZ">
178: /// (input) INTEGER
179: ///</param>
180: /// <param name="IHIZ">
181: /// (input) INTEGER
182: /// Specify the rows of Z to which transformations must be
183: /// applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
184: ///</param>
185: /// <param name="Z">
186: /// (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
187: /// IF WANTZ is .TRUE., then on output, the orthogonal
188: /// similarity transformation mentioned above has been
189: /// accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
190: /// If WANTZ is .FALSE., then Z is unreferenced.
191: ///</param>
192: /// <param name="LDZ">
193: /// (input) integer
194: /// The leading dimension of Z just as declared in the
195: /// calling subroutine. 1 .LE. LDZ.
196: ///</param>
197: /// <param name="NS">
198: /// (output) integer
199: /// The number of unconverged (ie approximate) eigenvalues
200: /// returned in SR and SI that may be used as shifts by the
201: /// calling subroutine.
202: ///</param>
203: /// <param name="ND">
204: /// (output) integer
205: /// The number of converged eigenvalues uncovered by this
206: /// subroutine.
207: ///</param>
208: /// <param name="SR">
209: /// (output) DOUBLE PRECISION array, dimension KBOT
210: ///</param>
211: /// <param name="SI">
212: /// (output) DOUBLE PRECISION array, dimension KBOT
213: /// On output, the real and imaginary parts of approximate
214: /// eigenvalues that may be used for shifts are stored in
215: /// SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
216: /// SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
217: /// The real and imaginary parts of converged eigenvalues
218: /// are stored in SR(KBOT-ND+1) through SR(KBOT) and
219: /// SI(KBOT-ND+1) through SI(KBOT), respectively.
220: ///</param>
221: /// <param name="V">
222: /// (workspace) DOUBLE PRECISION array, dimension (LDV,NW)
223: /// An NW-by-NW work array.
224: ///</param>
225: /// <param name="LDV">
226: /// (input) integer scalar
227: /// The leading dimension of V just as declared in the
228: /// calling subroutine. NW .LE. LDV
229: ///</param>
230: /// <param name="NH">
231: /// (input) integer scalar
232: /// The number of columns of T. NH.GE.NW.
233: ///</param>
234: /// <param name="T">
235: /// (workspace) DOUBLE PRECISION array, dimension (LDT,NW)
236: ///</param>
237: /// <param name="LDT">
238: /// (input) integer
239: /// The leading dimension of T just as declared in the
240: /// calling subroutine. NW .LE. LDT
241: ///</param>
242: /// <param name="NV">
243: /// (input) integer
244: /// The number of rows of work array WV available for
245: /// workspace. NV.GE.NW.
246: ///</param>
247: /// <param name="WV">
248: /// (workspace) DOUBLE PRECISION array, dimension (LDWV,NW)
249: ///</param>
250: /// <param name="LDWV">
251: /// (input) integer
252: /// The leading dimension of W just as declared in the
253: /// calling subroutine. NW .LE. LDV
254: ///</param>
255: /// <param name="WORK">
256: /// (workspace) DOUBLE PRECISION array, dimension LWORK.
257: /// On exit, WORK(1) is set to an estimate of the optimal value
258: /// of LWORK for the given values of N, NW, KTOP and KBOT.
259: ///</param>
260: /// <param name="LWORK">
261: /// (input) integer
262: /// The dimension of the work array WORK. LWORK = 2*NW
263: /// suffices, but greater efficiency may result from larger
264: /// values of LWORK.
265: ///
266: /// If LWORK = -1, then a workspace query is assumed; DLAQR2
267: /// only estimates the optimal workspace size for the given
268: /// values of N, NW, KTOP and KBOT. The estimate is returned
269: /// in WORK(1). No error message related to LWORK is issued
270: /// by XERBLA. Neither H nor Z are accessed.
271: ///</param>
272: public void Run(bool WANTT, bool WANTZ, int N, int KTOP, int KBOT, int NW
273: , ref double[] H, int offset_h, int LDH, int ILOZ, int IHIZ, ref double[] Z, int offset_z, int LDZ
274: , ref int NS, ref int ND, ref double[] SR, int offset_sr, ref double[] SI, int offset_si, ref double[] V, int offset_v, int LDV
275: , int NH, ref double[] T, int offset_t, int LDT, int NV, ref double[] WV, int offset_wv, int LDWV
276: , ref double[] WORK, int offset_work, int LWORK)
277: {
278:
279: #region Array Index Correction
280:
281: int o_h = -1 - LDH + offset_h; int o_z = -1 - LDZ + offset_z; int o_sr = -1 + offset_sr;
282: int o_si = -1 + offset_si; int o_v = -1 - LDV + offset_v; int o_t = -1 - LDT + offset_t;
283: int o_wv = -1 - LDWV + offset_wv; int o_work = -1 + offset_work;
284:
285: #endregion
286:
287:
288: #region Prolog
289:
290: // *
291: // * -- LAPACK auxiliary routine (version 3.1) --
292: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
293: // * November 2006
294: // *
295: // * .. Scalar Arguments ..
296: // * ..
297: // * .. Array Arguments ..
298: // * ..
299: // *
300: // * This subroutine is identical to DLAQR3 except that it avoids
301: // * recursion by calling DLAHQR instead of DLAQR4.
302: // *
303: // *
304: // * ******************************************************************
305: // * Aggressive early deflation:
306: // *
307: // * This subroutine accepts as input an upper Hessenberg matrix
308: // * H and performs an orthogonal similarity transformation
309: // * designed to detect and deflate fully converged eigenvalues from
310: // * a trailing principal submatrix. On output H has been over-
311: // * written by a new Hessenberg matrix that is a perturbation of
312: // * an orthogonal similarity transformation of H. It is to be
313: // * hoped that the final version of H has many zero subdiagonal
314: // * entries.
315: // *
316: // * ******************************************************************
317: // * WANTT (input) LOGICAL
318: // * If .TRUE., then the Hessenberg matrix H is fully updated
319: // * so that the quasi-triangular Schur factor may be
320: // * computed (in cooperation with the calling subroutine).
321: // * If .FALSE., then only enough of H is updated to preserve
322: // * the eigenvalues.
323: // *
324: // * WANTZ (input) LOGICAL
325: // * If .TRUE., then the orthogonal matrix Z is updated so
326: // * so that the orthogonal Schur factor may be computed
327: // * (in cooperation with the calling subroutine).
328: // * If .FALSE., then Z is not referenced.
329: // *
330: // * N (input) INTEGER
331: // * The order of the matrix H and (if WANTZ is .TRUE.) the
332: // * order of the orthogonal matrix Z.
333: // *
334: // * KTOP (input) INTEGER
335: // * It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
336: // * KBOT and KTOP together determine an isolated block
337: // * along the diagonal of the Hessenberg matrix.
338: // *
339: // * KBOT (input) INTEGER
340: // * It is assumed without a check that either
341: // * KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
342: // * determine an isolated block along the diagonal of the
343: // * Hessenberg matrix.
344: // *
345: // * NW (input) INTEGER
346: // * Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
347: // *
348: // * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
349: // * On input the initial N-by-N section of H stores the
350: // * Hessenberg matrix undergoing aggressive early deflation.
351: // * On output H has been transformed by an orthogonal
352: // * similarity transformation, perturbed, and the returned
353: // * to Hessenberg form that (it is to be hoped) has some
354: // * zero subdiagonal entries.
355: // *
356: // * LDH (input) integer
357: // * Leading dimension of H just as declared in the calling
358: // * subroutine. N .LE. LDH
359: // *
360: // * ILOZ (input) INTEGER
361: // * IHIZ (input) INTEGER
362: // * Specify the rows of Z to which transformations must be
363: // * applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
364: // *
365: // * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
366: // * IF WANTZ is .TRUE., then on output, the orthogonal
367: // * similarity transformation mentioned above has been
368: // * accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
369: // * If WANTZ is .FALSE., then Z is unreferenced.
370: // *
371: // * LDZ (input) integer
372: // * The leading dimension of Z just as declared in the
373: // * calling subroutine. 1 .LE. LDZ.
374: // *
375: // * NS (output) integer
376: // * The number of unconverged (ie approximate) eigenvalues
377: // * returned in SR and SI that may be used as shifts by the
378: // * calling subroutine.
379: // *
380: // * ND (output) integer
381: // * The number of converged eigenvalues uncovered by this
382: // * subroutine.
383: // *
384: // * SR (output) DOUBLE PRECISION array, dimension KBOT
385: // * SI (output) DOUBLE PRECISION array, dimension KBOT
386: // * On output, the real and imaginary parts of approximate
387: // * eigenvalues that may be used for shifts are stored in
388: // * SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
389: // * SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
390: // * The real and imaginary parts of converged eigenvalues
391: // * are stored in SR(KBOT-ND+1) through SR(KBOT) and
392: // * SI(KBOT-ND+1) through SI(KBOT), respectively.
393: // *
394: // * V (workspace) DOUBLE PRECISION array, dimension (LDV,NW)
395: // * An NW-by-NW work array.
396: // *
397: // * LDV (input) integer scalar
398: // * The leading dimension of V just as declared in the
399: // * calling subroutine. NW .LE. LDV
400: // *
401: // * NH (input) integer scalar
402: // * The number of columns of T. NH.GE.NW.
403: // *
404: // * T (workspace) DOUBLE PRECISION array, dimension (LDT,NW)
405: // *
406: // * LDT (input) integer
407: // * The leading dimension of T just as declared in the
408: // * calling subroutine. NW .LE. LDT
409: // *
410: // * NV (input) integer
411: // * The number of rows of work array WV available for
412: // * workspace. NV.GE.NW.
413: // *
414: // * WV (workspace) DOUBLE PRECISION array, dimension (LDWV,NW)
415: // *
416: // * LDWV (input) integer
417: // * The leading dimension of W just as declared in the
418: // * calling subroutine. NW .LE. LDV
419: // *
420: // * WORK (workspace) DOUBLE PRECISION array, dimension LWORK.
421: // * On exit, WORK(1) is set to an estimate of the optimal value
422: // * of LWORK for the given values of N, NW, KTOP and KBOT.
423: // *
424: // * LWORK (input) integer
425: // * The dimension of the work array WORK. LWORK = 2*NW
426: // * suffices, but greater efficiency may result from larger
427: // * values of LWORK.
428: // *
429: // * If LWORK = -1, then a workspace query is assumed; DLAQR2
430: // * only estimates the optimal workspace size for the given
431: // * values of N, NW, KTOP and KBOT. The estimate is returned
432: // * in WORK(1). No error message related to LWORK is issued
433: // * by XERBLA. Neither H nor Z are accessed.
434: // *
435: // * ================================================================
436: // * Based on contributions by
437: // * Karen Braman and Ralph Byers, Department of Mathematics,
438: // * University of Kansas, USA
439: // *
440: // * ================================================================
441: // * .. Parameters ..
442: // * ..
443: // * .. Local Scalars ..
444: // * ..
445: // * .. External Functions ..
446: // * ..
447: // * .. External Subroutines ..
448: // * ..
449: // * .. Intrinsic Functions ..
450: // INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT;
451: // * ..
452: // * .. Executable Statements ..
453: // *
454: // * ==== Estimate optimal workspace. ====
455: // *
456:
457: #endregion
458:
459:
460: #region Body
461:
462: JW = Math.Min(NW, KBOT - KTOP + 1);
463: if (JW <= 2)
464: {
465: LWKOPT = 1;
466: }
467: else
468: {
469: // *
470: // * ==== Workspace query call to DGEHRD ====
471: // *
472: this._dgehrd.Run(JW, 1, JW - 1, ref T, offset_t, LDT, ref WORK, offset_work
473: , ref WORK, offset_work, - 1, ref INFO);
474: LWK1 = Convert.ToInt32(Math.Truncate(WORK[1 + o_work]));
475: // *
476: // * ==== Workspace query call to DORGHR ====
477: // *
478: this._dorghr.Run(JW, 1, JW - 1, ref T, offset_t, LDT, WORK, offset_work
479: , ref WORK, offset_work, - 1, ref INFO);
480: LWK2 = Convert.ToInt32(Math.Truncate(WORK[1 + o_work]));
481: // *
482: // * ==== Optimal workspace ====
483: // *
484: LWKOPT = JW + Math.Max(LWK1, LWK2);
485: }
486: // *
487: // * ==== Quick return in case of workspace query. ====
488: // *
489: if (LWORK == - 1)
490: {
491: WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
492: return;
493: }
494: // *
495: // * ==== Nothing to do ...
496: // * ... for an empty active block ... ====
497: NS = 0;
498: ND = 0;
499: if (KTOP > KBOT) return;
500: // * ... nor for an empty deflation window. ====
501: if (NW < 1) return;
502: // *
503: // * ==== Machine constants ====
504: // *
505: SAFMIN = this._dlamch.Run("SAFE MINIMUM");
506: SAFMAX = ONE / SAFMIN;
507: this._dlabad.Run(ref SAFMIN, ref SAFMAX);
508: ULP = this._dlamch.Run("PRECISION");
509: SMLNUM = SAFMIN * (Convert.ToDouble(N) / ULP);
510: // *
511: // * ==== Setup deflation window ====
512: // *
513: JW = Math.Min(NW, KBOT - KTOP + 1);
514: KWTOP = KBOT - JW + 1;
515: if (KWTOP == KTOP)
516: {
517: S = ZERO;
518: }
519: else
520: {
521: S = H[KWTOP+(KWTOP - 1) * LDH + o_h];
522: }
523: // *
524: if (KBOT == KWTOP)
525: {
526: // *
527: // * ==== 1-by-1 deflation window: not much to do ====
528: // *
529: SR[KWTOP + o_sr] = H[KWTOP+KWTOP * LDH + o_h];
530: SI[KWTOP + o_si] = ZERO;
531: NS = 1;
532: ND = 0;
533: if (Math.Abs(S) <= Math.Max(SMLNUM, ULP * Math.Abs(H[KWTOP+KWTOP * LDH + o_h])))
534: {
535: NS = 0;
536: ND = 1;
537: if (KWTOP > KTOP) H[KWTOP+(KWTOP - 1) * LDH + o_h] = ZERO;
538: }
539: return;
540: }
541: // *
542: // * ==== Convert to spike-triangular form. (In case of a
543: // * . rare QR failure, this routine continues to do
544: // * . aggressive early deflation using that part of
545: // * . the deflation window that converged using INFQR
546: // * . here and there to keep track.) ====
547: // *
548: this._dlacpy.Run("U", JW, JW, H, KWTOP+KWTOP * LDH + o_h, LDH, ref T, offset_t
549: , LDT);
550: this._dcopy.Run(JW - 1, H, KWTOP + 1+KWTOP * LDH + o_h, LDH + 1, ref T, 2+1 * LDT + o_t, LDT + 1);
551: // *
552: this._dlaset.Run("A", JW, JW, ZERO, ONE, ref V, offset_v
553: , LDV);
554: this._dlahqr.Run(true, true, JW, 1, JW, ref T, offset_t
555: , LDT, ref SR, KWTOP + o_sr, ref SI, KWTOP + o_si, 1, JW, ref V, offset_v
556: , LDV, ref INFQR);
557: // *
558: // * ==== DTREXC needs a clean margin near the diagonal ====
559: // *
560: for (J = 1; J <= JW - 3; J++)
561: {
562: T[J + 2+J * LDT + o_t] = ZERO;
563: T[J + 3+J * LDT + o_t] = ZERO;
564: }
565: if (JW > 2) T[JW+(JW - 2) * LDT + o_t] = ZERO;
566: // *
567: // * ==== Deflation detection loop ====
568: // *
569: NS = JW;
570: ILST = INFQR + 1;
571: LABEL20:;
572: if (ILST <= NS)
573: {
574: if (NS == 1)
575: {
576: BULGE = false;
577: }
578: else
579: {
580: BULGE = T[NS+(NS - 1) * LDT + o_t] != ZERO;
581: }
582: // *
583: // * ==== Small spike tip test for deflation ====
584: // *
585: if (!BULGE)
586: {
587: // *
588: // * ==== Real eigenvalue ====
589: // *
590: FOO = Math.Abs(T[NS+NS * LDT + o_t]);
591: if (FOO == ZERO) FOO = Math.Abs(S);
592: if (Math.Abs(S * V[1+NS * LDV + o_v]) <= Math.Max(SMLNUM, ULP * FOO))
593: {
594: // *
595: // * ==== Deflatable ====
596: // *
597: NS = NS - 1;
598: }
599: else
600: {
601: // *
602: // * ==== Undeflatable. Move it up out of the way.
603: // * . (DTREXC can not fail in this case.) ====
604: // *
605: IFST = NS;
606: this._dtrexc.Run("V", JW, ref T, offset_t, LDT, ref V, offset_v, LDV
607: , ref IFST, ref ILST, ref WORK, offset_work, ref INFO);
608: ILST = ILST + 1;
609: }
610: }
611: else
612: {
613: // *
614: // * ==== Complex conjugate pair ====
615: // *
616: FOO = Math.Abs(T[NS+NS * LDT + o_t]) + Math.Sqrt(Math.Abs(T[NS+(NS - 1) * LDT + o_t])) * Math.Sqrt(Math.Abs(T[NS - 1+NS * LDT + o_t]));
617: if (FOO == ZERO) FOO = Math.Abs(S);
618: if (Math.Max(Math.Abs(S * V[1+NS * LDV + o_v]), Math.Abs(S * V[1+(NS - 1) * LDV + o_v])) <= Math.Max(SMLNUM, ULP * FOO))
619: {
620: // *
621: // * ==== Deflatable ====
622: // *
623: NS = NS - 2;
624: }
625: else
626: {
627: // *
628: // * ==== Undflatable. Move them up out of the way.
629: // * . Fortunately, DTREXC does the right thing with
630: // * . ILST in case of a rare exchange failure. ====
631: // *
632: IFST = NS;
633: this._dtrexc.Run("V", JW, ref T, offset_t, LDT, ref V, offset_v, LDV
634: , ref IFST, ref ILST, ref WORK, offset_work, ref INFO);
635: ILST = ILST + 2;
636: }
637: }
638: // *
639: // * ==== End deflation detection loop ====
640: // *
641: goto LABEL20;
642: }
643: // *
644: // * ==== Return to Hessenberg form ====
645: // *
646: if (NS == 0) S = ZERO;
647: // *
648: if (NS < JW)
649: {
650: // *
651: // * ==== sorting diagonal blocks of T improves accuracy for
652: // * . graded matrices. Bubble sort deals well with
653: // * . exchange failures. ====
654: // *
655: SORTED = false;
656: I = NS + 1;
657: LABEL30:;
658: if (SORTED) goto LABEL50;
659: SORTED = true;
660: // *
661: KEND = I - 1;
662: I = INFQR + 1;
663: if (I == NS)
664: {
665: K = I + 1;
666: }
667: else
668: {
669: if (T[I + 1+I * LDT + o_t] == ZERO)
670: {
671: K = I + 1;
672: }
673: else
674: {
675: K = I + 2;
676: }
677: }
678: LABEL40:;
679: if (K <= KEND)
680: {
681: if (K == I + 1)
682: {
683: EVI = Math.Abs(T[I+I * LDT + o_t]);
684: }
685: else
686: {
687: EVI = Math.Abs(T[I+I * LDT + o_t]) + Math.Sqrt(Math.Abs(T[I + 1+I * LDT + o_t])) * Math.Sqrt(Math.Abs(T[I+(I + 1) * LDT + o_t]));
688: }
689: // *
690: if (K == KEND)
691: {
692: EVK = Math.Abs(T[K+K * LDT + o_t]);
693: }
694: else
695: {
696: if (T[K + 1+K * LDT + o_t] == ZERO)
697: {
698: EVK = Math.Abs(T[K+K * LDT + o_t]);
699: }
700: else
701: {
702: EVK = Math.Abs(T[K+K * LDT + o_t]) + Math.Sqrt(Math.Abs(T[K + 1+K * LDT + o_t])) * Math.Sqrt(Math.Abs(T[K+(K + 1) * LDT + o_t]));
703: }
704: }
705: // *
706: if (EVI >= EVK)
707: {
708: I = K;
709: }
710: else
711: {
712: SORTED = false;
713: IFST = I;
714: ILST = K;
715: this._dtrexc.Run("V", JW, ref T, offset_t, LDT, ref V, offset_v, LDV
716: , ref IFST, ref ILST, ref WORK, offset_work, ref INFO);
717: if (INFO == 0)
718: {
719: I = ILST;
720: }
721: else
722: {
723: I = K;
724: }
725: }
726: if (I == KEND)
727: {
728: K = I + 1;
729: }
730: else
731: {
732: if (T[I + 1+I * LDT + o_t] == ZERO)
733: {
734: K = I + 1;
735: }
736: else
737: {
738: K = I + 2;
739: }
740: }
741: goto LABEL40;
742: }
743: goto LABEL30;
744: LABEL50:;
745: }
746: // *
747: // * ==== Restore shift/eigenvalue array from T ====
748: // *
749: I = JW;
750: LABEL60:;
751: if (I >= INFQR + 1)
752: {
753: if (I == INFQR + 1)
754: {
755: SR[KWTOP + I - 1 + o_sr] = T[I+I * LDT + o_t];
756: SI[KWTOP + I - 1 + o_si] = ZERO;
757: I = I - 1;
758: }
759: else
760: {
761: if (T[I+(I - 1) * LDT + o_t] == ZERO)
762: {
763: SR[KWTOP + I - 1 + o_sr] = T[I+I * LDT + o_t];
764: SI[KWTOP + I - 1 + o_si] = ZERO;
765: I = I - 1;
766: }
767: else
768: {
769: AA = T[I - 1+(I - 1) * LDT + o_t];
770: CC = T[I+(I - 1) * LDT + o_t];
771: BB = T[I - 1+I * LDT + o_t];
772: DD = T[I+I * LDT + o_t];
773: this._dlanv2.Run(ref AA, ref BB, ref CC, ref DD, ref SR[KWTOP + I - 2 + o_sr], ref SI[KWTOP + I - 2 + o_si]
774: , ref SR[KWTOP + I - 1 + o_sr], ref SI[KWTOP + I - 1 + o_si], ref CS, ref SN);
775: I = I - 2;
776: }
777: }
778: goto LABEL60;
779: }
780: // *
781: if (NS < JW || S == ZERO)
782: {
783: if (NS > 1 && S != ZERO)
784: {
785: // *
786: // * ==== Reflect spike back into lower triangle ====
787: // *
788: this._dcopy.Run(NS, V, offset_v, LDV, ref WORK, offset_work, 1);
789: BETA = WORK[1 + o_work];
790: this._dlarfg.Run(NS, ref BETA, ref WORK, 2 + o_work, 1, ref TAU);
791: WORK[1 + o_work] = ONE;
792: // *
793: this._dlaset.Run("L", JW - 2, JW - 2, ZERO, ZERO, ref T, 3+1 * LDT + o_t
794: , LDT);
795: // *
796: this._dlarf.Run("L", NS, JW, WORK, offset_work, 1, TAU
797: , ref T, offset_t, LDT, ref WORK, JW + 1 + o_work);
798: this._dlarf.Run("R", NS, NS, WORK, offset_work, 1, TAU
799: , ref T, offset_t, LDT, ref WORK, JW + 1 + o_work);
800: this._dlarf.Run("R", JW, NS, WORK, offset_work, 1, TAU
801: , ref V, offset_v, LDV, ref WORK, JW + 1 + o_work);
802: // *
803: this._dgehrd.Run(JW, 1, NS, ref T, offset_t, LDT, ref WORK, offset_work
804: , ref WORK, JW + 1 + o_work, LWORK - JW, ref INFO);
805: }
806: // *
807: // * ==== Copy updated reduced window into place ====
808: // *
809: if (KWTOP > 1) H[KWTOP+(KWTOP - 1) * LDH + o_h] = S * V[1+1 * LDV + o_v];
810: this._dlacpy.Run("U", JW, JW, T, offset_t, LDT, ref H, KWTOP+KWTOP * LDH + o_h
811: , LDH);
812: this._dcopy.Run(JW - 1, T, 2+1 * LDT + o_t, LDT + 1, ref H, KWTOP + 1+KWTOP * LDH + o_h, LDH + 1);
813: // *
814: // * ==== Accumulate orthogonal matrix in order update
815: // * . H and Z, if requested. (A modified version
816: // * . of DORGHR that accumulates block Householder
817: // * . transformations into V directly might be
818: // * . marginally more efficient than the following.) ====
819: // *
820: if (NS > 1 && S != ZERO)
821: {
822: this._dorghr.Run(JW, 1, NS, ref T, offset_t, LDT, WORK, offset_work
823: , ref WORK, JW + 1 + o_work, LWORK - JW, ref INFO);
824: this._dgemm.Run("N", "N", JW, NS, NS, ONE
825: , V, offset_v, LDV, T, offset_t, LDT, ZERO, ref WV, offset_wv
826: , LDWV);
827: this._dlacpy.Run("A", JW, NS, WV, offset_wv, LDWV, ref V, offset_v
828: , LDV);
829: }
830: // *
831: // * ==== Update vertical slab in H ====
832: // *
833: if (WANTT)
834: {
835: LTOP = 1;
836: }
837: else
838: {
839: LTOP = KTOP;
840: }
841: for (KROW = LTOP; (NV >= 0) ? (KROW <= KWTOP - 1) : (KROW >= KWTOP - 1); KROW += NV)
842: {
843: KLN = Math.Min(NV, KWTOP - KROW);
844: this._dgemm.Run("N", "N", KLN, JW, JW, ONE
845: , H, KROW+KWTOP * LDH + o_h, LDH, V, offset_v, LDV, ZERO, ref WV, offset_wv
846: , LDWV);
847: this._dlacpy.Run("A", KLN, JW, WV, offset_wv, LDWV, ref H, KROW+KWTOP * LDH + o_h
848: , LDH);
849: }
850: // *
851: // * ==== Update horizontal slab in H ====
852: // *
853: if (WANTT)
854: {
855: for (KCOL = KBOT + 1; (NH >= 0) ? (KCOL <= N) : (KCOL >= N); KCOL += NH)
856: {
857: KLN = Math.Min(NH, N - KCOL + 1);
858: this._dgemm.Run("C", "N", JW, KLN, JW, ONE
859: , V, offset_v, LDV, H, KWTOP+KCOL * LDH + o_h, LDH, ZERO, ref T, offset_t
860: , LDT);
861: this._dlacpy.Run("A", JW, KLN, T, offset_t, LDT, ref H, KWTOP+KCOL * LDH + o_h
862: , LDH);
863: }
864: }
865: // *
866: // * ==== Update vertical slab in Z ====
867: // *
868: if (WANTZ)
869: {
870: for (KROW = ILOZ; (NV >= 0) ? (KROW <= IHIZ) : (KROW >= IHIZ); KROW += NV)
871: {
872: KLN = Math.Min(NV, IHIZ - KROW + 1);
873: this._dgemm.Run("N", "N", KLN, JW, JW, ONE
874: , Z, KROW+KWTOP * LDZ + o_z, LDZ, V, offset_v, LDV, ZERO, ref WV, offset_wv
875: , LDWV);
876: this._dlacpy.Run("A", KLN, JW, WV, offset_wv, LDWV, ref Z, KROW+KWTOP * LDZ + o_z
877: , LDZ);
878: }
879: }
880: }
881: // *
882: // * ==== Return the number of deflations ... ====
883: // *
884: ND = JW - NS;
885: // *
886: // * ==== ... and the number of shifts. (Subtracting
887: // * . INFQR from the spike length takes care
888: // * . of the case of a rare QR failure while
889: // * . calculating eigenvalues of the deflation
890: // * . window.) ====
891: // *
892: NS = NS - INFQR;
893: // *
894: // * ==== Return optimal workspace. ====
895: // *
896: WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
897: // *
898: // * ==== End of DLAQR2 ====
899: // *
900:
901: #endregion
902:
903: }
904: }
905: }