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