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 DLAQR5
25: {
26:
27:
28: #region Dependencies
29:
30: DLAMCH _dlamch; DGEMM _dgemm; DLABAD _dlabad; DLACPY _dlacpy; DLAQR1 _dlaqr1; DLARFG _dlarfg; DLASET _dlaset;
31: DTRMM _dtrmm;
32:
33: #endregion
34:
35:
36: #region Fields
37:
38: const double ZERO = 0.0E0; const double ONE = 1.0E0; double ALPHA = 0; double BETA = 0; double H11 = 0; double H12 = 0;
39: double H21 = 0;double H22 = 0; double REFSUM = 0; double SAFMAX = 0; double SAFMIN = 0; double SCL = 0; double SMLNUM = 0;
40: double SWAP = 0;double TST1 = 0; double TST2 = 0; double ULP = 0; int I = 0; int I2 = 0; int I4 = 0; int INCOL = 0;
41: int J = 0;int J2 = 0; int J4 = 0; int JBOT = 0; int JCOL = 0; int JLEN = 0; int JROW = 0; int JTOP = 0; int K = 0;
42: int K1 = 0;int KDU = 0; int KMS = 0; int KNZ = 0; int KRCOL = 0; int KZS = 0; int M = 0; int M22 = 0; int MBOT = 0;
43: int MEND = 0;int MSTART = 0; int MTOP = 0; int NBMPS = 0; int NDCOL = 0; int NS = 0; int NU = 0; bool ACCUM = false;
44: bool BLK22 = false;bool BMP22 = false; double[] VT = new double[3]; int offset_vt = 0; int o_vt = -1;
45:
46: #endregion
47:
48: public DLAQR5(DLAMCH dlamch, DGEMM dgemm, DLABAD dlabad, DLACPY dlacpy, DLAQR1 dlaqr1, DLARFG dlarfg, DLASET dlaset, DTRMM dtrmm)
49: {
50:
51:
52: #region Set Dependencies
53:
54: this._dlamch = dlamch; this._dgemm = dgemm; this._dlabad = dlabad; this._dlacpy = dlacpy; this._dlaqr1 = dlaqr1;
55: this._dlarfg = dlarfg;this._dlaset = dlaset; this._dtrmm = dtrmm;
56:
57: #endregion
58:
59: }
60:
61: public DLAQR5()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: LSAME lsame = new LSAME();
68: DLAMC3 dlamc3 = new DLAMC3();
69: XERBLA xerbla = new XERBLA();
70: DLABAD dlabad = new DLABAD();
71: DLAQR1 dlaqr1 = new DLAQR1();
72: DLAPY2 dlapy2 = new DLAPY2();
73: DNRM2 dnrm2 = new DNRM2();
74: DSCAL dscal = new DSCAL();
75: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
76: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
77: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
78: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
79: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
80: DGEMM dgemm = new DGEMM(lsame, xerbla);
81: DLACPY dlacpy = new DLACPY(lsame);
82: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
83: DLASET dlaset = new DLASET(lsame);
84: DTRMM dtrmm = new DTRMM(lsame, xerbla);
85:
86: #endregion
87:
88:
89: #region Set Dependencies
90:
91: this._dlamch = dlamch; this._dgemm = dgemm; this._dlabad = dlabad; this._dlacpy = dlacpy; this._dlaqr1 = dlaqr1;
92: this._dlarfg = dlarfg;this._dlaset = dlaset; this._dtrmm = dtrmm;
93:
94: #endregion
95:
96: }
97: /// <param name="WANTT">
98: /// (input) logical scalar
99: /// WANTT = .true. if the quasi-triangular Schur factor
100: /// is being computed. WANTT is set to .false. otherwise.
101: ///</param>
102: /// <param name="WANTZ">
103: /// (input) logical scalar
104: /// WANTZ = .true. if the orthogonal Schur factor is being
105: /// computed. WANTZ is set to .false. otherwise.
106: ///</param>
107: /// <param name="KACC22">
108: /// (input) integer with value 0, 1, or 2.
109: /// Specifies the computation mode of far-from-diagonal
110: /// orthogonal updates.
111: /// = 0: DLAQR5 does not accumulate reflections and does not
112: /// use matrix-matrix multiply to update far-from-diagonal
113: /// matrix entries.
114: /// = 1: DLAQR5 accumulates reflections and uses matrix-matrix
115: /// multiply to update the far-from-diagonal matrix entries.
116: /// = 2: DLAQR5 accumulates reflections, uses matrix-matrix
117: /// multiply to update the far-from-diagonal matrix entries,
118: /// and takes advantage of 2-by-2 block structure during
119: /// matrix multiplies.
120: ///</param>
121: /// <param name="N">
122: /// (input) integer scalar
123: /// N is the order of the Hessenberg matrix H upon which this
124: /// subroutine operates.
125: ///</param>
126: /// <param name="KTOP">
127: /// (input) integer scalar
128: ///</param>
129: /// <param name="KBOT">
130: /// (input) integer scalar
131: /// These are the first and last rows and columns of an
132: /// isolated diagonal block upon which the QR sweep is to be
133: /// applied. It is assumed without a check that
134: /// either KTOP = 1 or H(KTOP,KTOP-1) = 0
135: /// and
136: /// either KBOT = N or H(KBOT+1,KBOT) = 0.
137: ///</param>
138: /// <param name="NSHFTS">
139: /// (input) integer scalar
140: /// NSHFTS gives the number of simultaneous shifts. NSHFTS
141: /// must be positive and even.
142: ///</param>
143: /// <param name="SR">
144: /// (input) DOUBLE PRECISION array of size (NSHFTS)
145: ///</param>
146: /// <param name="SI">
147: /// (input) DOUBLE PRECISION array of size (NSHFTS)
148: /// SR contains the real parts and SI contains the imaginary
149: /// parts of the NSHFTS shifts of origin that define the
150: /// multi-shift QR sweep.
151: ///</param>
152: /// <param name="H">
153: /// (input/output) DOUBLE PRECISION array of size (LDH,N)
154: /// On input H contains a Hessenberg matrix. On output a
155: /// multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
156: /// to the isolated diagonal block in rows and columns KTOP
157: /// through KBOT.
158: ///</param>
159: /// <param name="LDH">
160: /// (input) integer scalar
161: /// LDH is the leading dimension of H just as declared in the
162: /// calling procedure. LDH.GE.MAX(1,N).
163: ///</param>
164: /// <param name="ILOZ">
165: /// (input) INTEGER
166: ///</param>
167: /// <param name="IHIZ">
168: /// (input) INTEGER
169: /// Specify the rows of Z to which transformations must be
170: /// applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
171: ///</param>
172: /// <param name="Z">
173: /// (input/output) DOUBLE PRECISION array of size (LDZ,IHI)
174: /// If WANTZ = .TRUE., then the QR Sweep orthogonal
175: /// similarity transformation is accumulated into
176: /// Z(ILOZ:IHIZ,ILO:IHI) from the right.
177: /// If WANTZ = .FALSE., then Z is unreferenced.
178: ///</param>
179: /// <param name="LDZ">
180: /// (input) integer scalar
181: /// LDA is the leading dimension of Z just as declared in
182: /// the calling procedure. LDZ.GE.N.
183: ///</param>
184: /// <param name="V">
185: /// (workspace) DOUBLE PRECISION array of size (LDV,NSHFTS/2)
186: ///</param>
187: /// <param name="LDV">
188: /// (input) integer scalar
189: /// LDV is the leading dimension of V as declared in the
190: /// calling procedure. LDV.GE.3.
191: ///</param>
192: /// <param name="U">
193: /// (workspace) DOUBLE PRECISION array of size
194: /// (LDU,3*NSHFTS-3)
195: ///</param>
196: /// <param name="LDU">
197: /// (input) integer scalar
198: /// LDU is the leading dimension of U just as declared in the
199: /// in the calling subroutine. LDU.GE.3*NSHFTS-3.
200: ///</param>
201: /// <param name="NV">
202: /// (input) integer scalar
203: /// NV is the number of rows in WV agailable for workspace.
204: /// NV.GE.1.
205: ///</param>
206: /// <param name="WV">
207: /// (workspace) DOUBLE PRECISION array of size
208: /// (LDWV,3*NSHFTS-3)
209: ///</param>
210: /// <param name="LDWV">
211: /// (input) integer scalar
212: /// LDWV is the leading dimension of WV as declared in the
213: /// in the calling subroutine. LDWV.GE.NV.
214: ///
215: ///</param>
216: /// <param name="NH">
217: /// (input) integer scalar
218: /// NH is the number of columns in array WH available for
219: /// workspace. NH.GE.1.
220: ///</param>
221: /// <param name="WH">
222: /// (workspace) DOUBLE PRECISION array of size (LDWH,NH)
223: ///</param>
224: /// <param name="LDWH">
225: /// (input) integer scalar
226: /// Leading dimension of WH just as declared in the
227: /// calling procedure. LDWH.GE.3*NSHFTS-3.
228: ///</param>
229: public void Run(bool WANTT, bool WANTZ, int KACC22, int N, int KTOP, int KBOT
230: , int NSHFTS, ref double[] SR, int offset_sr, ref double[] SI, int offset_si, ref double[] H, int offset_h, int LDH, int ILOZ
231: , int IHIZ, ref double[] Z, int offset_z, int LDZ, ref double[] V, int offset_v, int LDV, ref double[] U, int offset_u
232: , int LDU, int NV, ref double[] WV, int offset_wv, int LDWV, int NH, ref double[] WH, int offset_wh
233: , int LDWH)
234: {
235:
236: #region Array Index Correction
237:
238: int o_sr = -1 + offset_sr; int o_si = -1 + offset_si; int o_h = -1 - LDH + offset_h;
239: int o_z = -1 - LDZ + offset_z; int o_v = -1 - LDV + offset_v; int o_u = -1 - LDU + offset_u;
240: int o_wv = -1 - LDWV + offset_wv; int o_wh = -1 - LDWH + offset_wh;
241:
242: #endregion
243:
244:
245: #region Prolog
246:
247: // *
248: // * -- LAPACK auxiliary routine (version 3.1) --
249: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
250: // * November 2006
251: // *
252: // * .. Scalar Arguments ..
253: // * ..
254: // * .. Array Arguments ..
255: // * ..
256: // *
257: // * This auxiliary subroutine called by DLAQR0 performs a
258: // * single small-bulge multi-shift QR sweep.
259: // *
260: // * WANTT (input) logical scalar
261: // * WANTT = .true. if the quasi-triangular Schur factor
262: // * is being computed. WANTT is set to .false. otherwise.
263: // *
264: // * WANTZ (input) logical scalar
265: // * WANTZ = .true. if the orthogonal Schur factor is being
266: // * computed. WANTZ is set to .false. otherwise.
267: // *
268: // * KACC22 (input) integer with value 0, 1, or 2.
269: // * Specifies the computation mode of far-from-diagonal
270: // * orthogonal updates.
271: // * = 0: DLAQR5 does not accumulate reflections and does not
272: // * use matrix-matrix multiply to update far-from-diagonal
273: // * matrix entries.
274: // * = 1: DLAQR5 accumulates reflections and uses matrix-matrix
275: // * multiply to update the far-from-diagonal matrix entries.
276: // * = 2: DLAQR5 accumulates reflections, uses matrix-matrix
277: // * multiply to update the far-from-diagonal matrix entries,
278: // * and takes advantage of 2-by-2 block structure during
279: // * matrix multiplies.
280: // *
281: // * N (input) integer scalar
282: // * N is the order of the Hessenberg matrix H upon which this
283: // * subroutine operates.
284: // *
285: // * KTOP (input) integer scalar
286: // * KBOT (input) integer scalar
287: // * These are the first and last rows and columns of an
288: // * isolated diagonal block upon which the QR sweep is to be
289: // * applied. It is assumed without a check that
290: // * either KTOP = 1 or H(KTOP,KTOP-1) = 0
291: // * and
292: // * either KBOT = N or H(KBOT+1,KBOT) = 0.
293: // *
294: // * NSHFTS (input) integer scalar
295: // * NSHFTS gives the number of simultaneous shifts. NSHFTS
296: // * must be positive and even.
297: // *
298: // * SR (input) DOUBLE PRECISION array of size (NSHFTS)
299: // * SI (input) DOUBLE PRECISION array of size (NSHFTS)
300: // * SR contains the real parts and SI contains the imaginary
301: // * parts of the NSHFTS shifts of origin that define the
302: // * multi-shift QR sweep.
303: // *
304: // * H (input/output) DOUBLE PRECISION array of size (LDH,N)
305: // * On input H contains a Hessenberg matrix. On output a
306: // * multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
307: // * to the isolated diagonal block in rows and columns KTOP
308: // * through KBOT.
309: // *
310: // * LDH (input) integer scalar
311: // * LDH is the leading dimension of H just as declared in the
312: // * calling procedure. LDH.GE.MAX(1,N).
313: // *
314: // * ILOZ (input) INTEGER
315: // * IHIZ (input) INTEGER
316: // * Specify the rows of Z to which transformations must be
317: // * applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
318: // *
319: // * Z (input/output) DOUBLE PRECISION array of size (LDZ,IHI)
320: // * If WANTZ = .TRUE., then the QR Sweep orthogonal
321: // * similarity transformation is accumulated into
322: // * Z(ILOZ:IHIZ,ILO:IHI) from the right.
323: // * If WANTZ = .FALSE., then Z is unreferenced.
324: // *
325: // * LDZ (input) integer scalar
326: // * LDA is the leading dimension of Z just as declared in
327: // * the calling procedure. LDZ.GE.N.
328: // *
329: // * V (workspace) DOUBLE PRECISION array of size (LDV,NSHFTS/2)
330: // *
331: // * LDV (input) integer scalar
332: // * LDV is the leading dimension of V as declared in the
333: // * calling procedure. LDV.GE.3.
334: // *
335: // * U (workspace) DOUBLE PRECISION array of size
336: // * (LDU,3*NSHFTS-3)
337: // *
338: // * LDU (input) integer scalar
339: // * LDU is the leading dimension of U just as declared in the
340: // * in the calling subroutine. LDU.GE.3*NSHFTS-3.
341: // *
342: // * NH (input) integer scalar
343: // * NH is the number of columns in array WH available for
344: // * workspace. NH.GE.1.
345: // *
346: // * WH (workspace) DOUBLE PRECISION array of size (LDWH,NH)
347: // *
348: // * LDWH (input) integer scalar
349: // * Leading dimension of WH just as declared in the
350: // * calling procedure. LDWH.GE.3*NSHFTS-3.
351: // *
352: // * NV (input) integer scalar
353: // * NV is the number of rows in WV agailable for workspace.
354: // * NV.GE.1.
355: // *
356: // * WV (workspace) DOUBLE PRECISION array of size
357: // * (LDWV,3*NSHFTS-3)
358: // *
359: // * LDWV (input) integer scalar
360: // * LDWV is the leading dimension of WV as declared in the
361: // * in the calling subroutine. LDWV.GE.NV.
362: // *
363: // *
364: // * ================================================================
365: // * Based on contributions by
366: // * Karen Braman and Ralph Byers, Department of Mathematics,
367: // * University of Kansas, USA
368: // *
369: // * ============================================================
370: // * Reference:
371: // *
372: // * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
373: // * Algorithm Part I: Maintaining Well Focused Shifts, and
374: // * Level 3 Performance, SIAM Journal of Matrix Analysis,
375: // * volume 23, pages 929--947, 2002.
376: // *
377: // * ============================================================
378: // * .. Parameters ..
379: // * ..
380: // * .. Local Scalars ..
381: // * ..
382: // * .. External Functions ..
383: // * ..
384: // * .. Intrinsic Functions ..
385: // *
386: // INTRINSIC ABS, DBLE, MAX, MIN, MOD;
387: // * ..
388: // * .. Local Arrays ..
389: // * ..
390: // * .. External Subroutines ..
391: // * ..
392: // * .. Executable Statements ..
393: // *
394: // * ==== If there are no shifts, then there is nothing to do. ====
395: // *
396:
397: #endregion
398:
399:
400: #region Body
401:
402: if (NSHFTS < 2) return;
403: // *
404: // * ==== If the active block is empty or 1-by-1, then there
405: // * . is nothing to do. ====
406: // *
407: if (KTOP >= KBOT) return;
408: // *
409: // * ==== Shuffle shifts into pairs of real shifts and pairs
410: // * . of complex conjugate shifts assuming complex
411: // * . conjugate shifts are already adjacent to one
412: // * . another. ====
413: // *
414: for (I = 1; I <= NSHFTS - 2; I += 2)
415: {
416: if (SI[I + o_si] != - SI[I + 1 + o_si])
417: {
418: // *
419: SWAP = SR[I + o_sr];
420: SR[I + o_sr] = SR[I + 1 + o_sr];
421: SR[I + 1 + o_sr] = SR[I + 2 + o_sr];
422: SR[I + 2 + o_sr] = SWAP;
423: // *
424: SWAP = SI[I + o_si];
425: SI[I + o_si] = SI[I + 1 + o_si];
426: SI[I + 1 + o_si] = SI[I + 2 + o_si];
427: SI[I + 2 + o_si] = SWAP;
428: }
429: }
430: // *
431: // * ==== NSHFTS is supposed to be even, but if is odd,
432: // * . then simply reduce it by one. The shuffle above
433: // * . ensures that the dropped shift is real and that
434: // * . the remaining shifts are paired. ====
435: // *
436: NS = NSHFTS - FortranLib.Mod(NSHFTS,2);
437: // *
438: // * ==== Machine constants for deflation ====
439: // *
440: SAFMIN = this._dlamch.Run("SAFE MINIMUM");
441: SAFMAX = ONE / SAFMIN;
442: this._dlabad.Run(ref SAFMIN, ref SAFMAX);
443: ULP = this._dlamch.Run("PRECISION");
444: SMLNUM = SAFMIN * (Convert.ToDouble(N) / ULP);
445: // *
446: // * ==== Use accumulated reflections to update far-from-diagonal
447: // * . entries ? ====
448: // *
449: ACCUM = (KACC22 == 1) || (KACC22 == 2);
450: // *
451: // * ==== If so, exploit the 2-by-2 block structure? ====
452: // *
453: BLK22 = (NS > 2) && (KACC22 == 2);
454: // *
455: // * ==== clear trash ====
456: // *
457: if (KTOP + 2 <= KBOT) H[KTOP + 2+KTOP * LDH + o_h] = ZERO;
458: // *
459: // * ==== NBMPS = number of 2-shift bulges in the chain ====
460: // *
461: NBMPS = NS / 2;
462: // *
463: // * ==== KDU = width of slab ====
464: // *
465: KDU = 6 * NBMPS - 3;
466: // *
467: // * ==== Create and chase chains of NBMPS bulges ====
468: // *
469: for (INCOL = 3 * (1 - NBMPS) + KTOP - 1; (3 * NBMPS - 2 >= 0) ? (INCOL <= KBOT - 2) : (INCOL >= KBOT - 2); INCOL += 3 * NBMPS - 2)
470: {
471: NDCOL = INCOL + KDU;
472: if (ACCUM)
473: {
474: this._dlaset.Run("ALL", KDU, KDU, ZERO, ONE, ref U, offset_u
475: , LDU);
476: }
477: // *
478: // * ==== Near-the-diagonal bulge chase. The following loop
479: // * . performs the near-the-diagonal part of a small bulge
480: // * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
481: // * . chunk extends from column INCOL to column NDCOL
482: // * . (including both column INCOL and column NDCOL). The
483: // * . following loop chases a 3*NBMPS column long chain of
484: // * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
485: // * . may be less than KTOP and and NDCOL may be greater than
486: // * . KBOT indicating phantom columns from which to chase
487: // * . bulges before they are actually introduced or to which
488: // * . to chase bulges beyond column KBOT.) ====
489: // *
490: for (KRCOL = INCOL; KRCOL <= Math.Min(INCOL + 3 * NBMPS - 3, KBOT - 2); KRCOL++)
491: {
492: // *
493: // * ==== Bulges number MTOP to MBOT are active double implicit
494: // * . shift bulges. There may or may not also be small
495: // * . 2-by-2 bulge, if there is room. The inactive bulges
496: // * . (if any) must wait until the active bulges have moved
497: // * . down the diagonal to make room. The phantom matrix
498: // * . paradigm described above helps keep track. ====
499: // *
500: MTOP = Math.Max(1, ((KTOP - 1) - KRCOL + 2) / 3 + 1);
501: MBOT = Math.Min(NBMPS, (KBOT - KRCOL) / 3);
502: M22 = MBOT + 1;
503: BMP22 = (MBOT < NBMPS) && (KRCOL + 3 * (M22 - 1)) == (KBOT - 2);
504: // *
505: // * ==== Generate reflections to chase the chain right
506: // * . one column. (The minimum value of K is KTOP-1.) ====
507: // *
508: for (M = MTOP; M <= MBOT; M++)
509: {
510: K = KRCOL + 3 * (M - 1);
511: if (K == KTOP - 1)
512: {
513: this._dlaqr1.Run(3, H, KTOP+KTOP * LDH + o_h, LDH, SR[2 * M - 1 + o_sr], SI[2 * M - 1 + o_si], SR[2 * M + o_sr]
514: , SI[2 * M + o_si], ref V, 1+M * LDV + o_v);
515: ALPHA = V[1+M * LDV + o_v];
516: this._dlarfg.Run(3, ref ALPHA, ref V, 2+M * LDV + o_v, 1, ref V[1+M * LDV + o_v]);
517: }
518: else
519: {
520: BETA = H[K + 1+K * LDH + o_h];
521: V[2+M * LDV + o_v] = H[K + 2+K * LDH + o_h];
522: V[3+M * LDV + o_v] = H[K + 3+K * LDH + o_h];
523: this._dlarfg.Run(3, ref BETA, ref V, 2+M * LDV + o_v, 1, ref V[1+M * LDV + o_v]);
524: // *
525: // * ==== A Bulge may collapse because of vigilant
526: // * . deflation or destructive underflow. (The
527: // * . initial bulge is always collapsed.) Use
528: // * . the two-small-subdiagonals trick to try
529: // * . to get it started again. If V(2,M).NE.0 and
530: // * . V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
531: // * . this bulge is collapsing into a zero
532: // * . subdiagonal. It will be restarted next
533: // * . trip through the loop.)
534: // *
535: if (V[1+M * LDV + o_v] != ZERO && (V[3+M * LDV + o_v] != ZERO || (H[K + 3+(K + 1) * LDH + o_h] == ZERO && H[K + 3+(K + 2) * LDH + o_h] == ZERO)))
536: {
537: // *
538: // * ==== Typical case: not collapsed (yet). ====
539: // *
540: H[K + 1+K * LDH + o_h] = BETA;
541: H[K + 2+K * LDH + o_h] = ZERO;
542: H[K + 3+K * LDH + o_h] = ZERO;
543: }
544: else
545: {
546: // *
547: // * ==== Atypical case: collapsed. Attempt to
548: // * . reintroduce ignoring H(K+1,K). If the
549: // * . fill resulting from the new reflector
550: // * . is too large, then abandon it.
551: // * . Otherwise, use the new one. ====
552: // *
553: this._dlaqr1.Run(3, H, K + 1+(K + 1) * LDH + o_h, LDH, SR[2 * M - 1 + o_sr], SI[2 * M - 1 + o_si], SR[2 * M + o_sr]
554: , SI[2 * M + o_si], ref VT, offset_vt);
555: SCL = Math.Abs(VT[1 + o_vt]) + Math.Abs(VT[2 + o_vt]) + Math.Abs(VT[3 + o_vt]);
556: if (SCL != ZERO)
557: {
558: VT[1 + o_vt] = VT[1 + o_vt] / SCL;
559: VT[2 + o_vt] = VT[2 + o_vt] / SCL;
560: VT[3 + o_vt] = VT[3 + o_vt] / SCL;
561: }
562: // *
563: // * ==== The following is the traditional and
564: // * . conservative two-small-subdiagonals
565: // * . test. ====
566: // * .
567: if (Math.Abs(H[K + 1+K * LDH + o_h]) * (Math.Abs(VT[2 + o_vt]) + Math.Abs(VT[3 + o_vt])) > ULP * Math.Abs(VT[1 + o_vt]) * (Math.Abs(H[K+K * LDH + o_h]) + Math.Abs(H[K + 1+(K + 1) * LDH + o_h]) + Math.Abs(H[K + 2+(K + 2) * LDH + o_h])))
568: {
569: // *
570: // * ==== Starting a new bulge here would
571: // * . create non-negligible fill. If
572: // * . the old reflector is diagonal (only
573: // * . possible with underflows), then
574: // * . change it to I. Otherwise, use
575: // * . it with trepidation. ====
576: // *
577: if (V[2+M * LDV + o_v] == ZERO && V[3+M * LDV + o_v] == ZERO)
578: {
579: V[1+M * LDV + o_v] = ZERO;
580: }
581: else
582: {
583: H[K + 1+K * LDH + o_h] = BETA;
584: H[K + 2+K * LDH + o_h] = ZERO;
585: H[K + 3+K * LDH + o_h] = ZERO;
586: }
587: }
588: else
589: {
590: // *
591: // * ==== Stating a new bulge here would
592: // * . create only negligible fill.
593: // * . Replace the old reflector with
594: // * . the new one. ====
595: // *
596: ALPHA = VT[1 + o_vt];
597: this._dlarfg.Run(3, ref ALPHA, ref VT, 2 + o_vt, 1, ref VT[1 + o_vt]);
598: REFSUM = H[K + 1+K * LDH + o_h] + H[K + 2+K * LDH + o_h] * VT[2 + o_vt] + H[K + 3+K * LDH + o_h] * VT[3 + o_vt];
599: H[K + 1+K * LDH + o_h] = H[K + 1+K * LDH + o_h] - VT[1 + o_vt] * REFSUM;
600: H[K + 2+K * LDH + o_h] = ZERO;
601: H[K + 3+K * LDH + o_h] = ZERO;
602: V[1+M * LDV + o_v] = VT[1 + o_vt];
603: V[2+M * LDV + o_v] = VT[2 + o_vt];
604: V[3+M * LDV + o_v] = VT[3 + o_vt];
605: }
606: }
607: }
608: }
609: // *
610: // * ==== Generate a 2-by-2 reflection, if needed. ====
611: // *
612: K = KRCOL + 3 * (M22 - 1);
613: if (BMP22)
614: {
615: if (K == KTOP - 1)
616: {
617: this._dlaqr1.Run(2, H, K + 1+(K + 1) * LDH + o_h, LDH, SR[2 * M22 - 1 + o_sr], SI[2 * M22 - 1 + o_si], SR[2 * M22 + o_sr]
618: , SI[2 * M22 + o_si], ref V, 1+M22 * LDV + o_v);
619: BETA = V[1+M22 * LDV + o_v];
620: this._dlarfg.Run(2, ref BETA, ref V, 2+M22 * LDV + o_v, 1, ref V[1+M22 * LDV + o_v]);
621: }
622: else
623: {
624: BETA = H[K + 1+K * LDH + o_h];
625: V[2+M22 * LDV + o_v] = H[K + 2+K * LDH + o_h];
626: this._dlarfg.Run(2, ref BETA, ref V, 2+M22 * LDV + o_v, 1, ref V[1+M22 * LDV + o_v]);
627: H[K + 1+K * LDH + o_h] = BETA;
628: H[K + 2+K * LDH + o_h] = ZERO;
629: }
630: }
631: else
632: {
633: // *
634: // * ==== Initialize V(1,M22) here to avoid possible undefined
635: // * . variable problems later. ====
636: // *
637: V[1+M22 * LDV + o_v] = ZERO;
638: }
639: // *
640: // * ==== Multiply H by reflections from the left ====
641: // *
642: if (ACCUM)
643: {
644: JBOT = Math.Min(NDCOL, KBOT);
645: }
646: else
647: {
648: if (WANTT)
649: {
650: JBOT = N;
651: }
652: else
653: {
654: JBOT = KBOT;
655: }
656: }
657: for (J = Math.Max(KTOP, KRCOL); J <= JBOT; J++)
658: {
659: MEND = Math.Min(MBOT, (J - KRCOL + 2) / 3);
660: for (M = MTOP; M <= MEND; M++)
661: {
662: K = KRCOL + 3 * (M - 1);
663: REFSUM = V[1+M * LDV + o_v] * (H[K + 1+J * LDH + o_h] + V[2+M * LDV + o_v] * H[K + 2+J * LDH + o_h] + V[3+M * LDV + o_v] * H[K + 3+J * LDH + o_h]);
664: H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - REFSUM;
665: H[K + 2+J * LDH + o_h] = H[K + 2+J * LDH + o_h] - REFSUM * V[2+M * LDV + o_v];
666: H[K + 3+J * LDH + o_h] = H[K + 3+J * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
667: }
668: }
669: if (BMP22)
670: {
671: K = KRCOL + 3 * (M22 - 1);
672: for (J = Math.Max(K + 1, KTOP); J <= JBOT; J++)
673: {
674: REFSUM = V[1+M22 * LDV + o_v] * (H[K + 1+J * LDH + o_h] + V[2+M22 * LDV + o_v] * H[K + 2+J * LDH + o_h]);
675: H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - REFSUM;
676: H[K + 2+J * LDH + o_h] = H[K + 2+J * LDH + o_h] - REFSUM * V[2+M22 * LDV + o_v];
677: }
678: }
679: // *
680: // * ==== Multiply H by reflections from the right.
681: // * . Delay filling in the last row until the
682: // * . vigilant deflation check is complete. ====
683: // *
684: if (ACCUM)
685: {
686: JTOP = Math.Max(KTOP, INCOL);
687: }
688: else
689: {
690: if (WANTT)
691: {
692: JTOP = 1;
693: }
694: else
695: {
696: JTOP = KTOP;
697: }
698: }
699: for (M = MTOP; M <= MBOT; M++)
700: {
701: if (V[1+M * LDV + o_v] != ZERO)
702: {
703: K = KRCOL + 3 * (M - 1);
704: for (J = JTOP; J <= Math.Min(KBOT, K + 3); J++)
705: {
706: REFSUM = V[1+M * LDV + o_v] * (H[J+(K + 1) * LDH + o_h] + V[2+M * LDV + o_v] * H[J+(K + 2) * LDH + o_h] + V[3+M * LDV + o_v] * H[J+(K + 3) * LDH + o_h]);
707: H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - REFSUM;
708: H[J+(K + 2) * LDH + o_h] = H[J+(K + 2) * LDH + o_h] - REFSUM * V[2+M * LDV + o_v];
709: H[J+(K + 3) * LDH + o_h] = H[J+(K + 3) * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
710: }
711: // *
712: if (ACCUM)
713: {
714: // *
715: // * ==== Accumulate U. (If necessary, update Z later
716: // * . with with an efficient matrix-matrix
717: // * . multiply.) ====
718: // *
719: KMS = K - INCOL;
720: for (J = Math.Max(1, KTOP - INCOL); J <= KDU; J++)
721: {
722: REFSUM = V[1+M * LDV + o_v] * (U[J+(KMS + 1) * LDU + o_u] + V[2+M * LDV + o_v] * U[J+(KMS + 2) * LDU + o_u] + V[3+M * LDV + o_v] * U[J+(KMS + 3) * LDU + o_u]);
723: U[J+(KMS + 1) * LDU + o_u] = U[J+(KMS + 1) * LDU + o_u] - REFSUM;
724: U[J+(KMS + 2) * LDU + o_u] = U[J+(KMS + 2) * LDU + o_u] - REFSUM * V[2+M * LDV + o_v];
725: U[J+(KMS + 3) * LDU + o_u] = U[J+(KMS + 3) * LDU + o_u] - REFSUM * V[3+M * LDV + o_v];
726: }
727: }
728: else
729: {
730: if (WANTZ)
731: {
732: // *
733: // * ==== U is not accumulated, so update Z
734: // * . now by multiplying by reflections
735: // * . from the right. ====
736: // *
737: for (J = ILOZ; J <= IHIZ; J++)
738: {
739: REFSUM = V[1+M * LDV + o_v] * (Z[J+(K + 1) * LDZ + o_z] + V[2+M * LDV + o_v] * Z[J+(K + 2) * LDZ + o_z] + V[3+M * LDV + o_v] * Z[J+(K + 3) * LDZ + o_z]);
740: Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - REFSUM;
741: Z[J+(K + 2) * LDZ + o_z] = Z[J+(K + 2) * LDZ + o_z] - REFSUM * V[2+M * LDV + o_v];
742: Z[J+(K + 3) * LDZ + o_z] = Z[J+(K + 3) * LDZ + o_z] - REFSUM * V[3+M * LDV + o_v];
743: }
744: }
745: }
746: }
747: }
748: // *
749: // * ==== Special case: 2-by-2 reflection (if needed) ====
750: // *
751: K = KRCOL + 3 * (M22 - 1);
752: if (BMP22 && (V[1+M22 * LDV + o_v] != ZERO))
753: {
754: for (J = JTOP; J <= Math.Min(KBOT, K + 3); J++)
755: {
756: REFSUM = V[1+M22 * LDV + o_v] * (H[J+(K + 1) * LDH + o_h] + V[2+M22 * LDV + o_v] * H[J+(K + 2) * LDH + o_h]);
757: H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - REFSUM;
758: H[J+(K + 2) * LDH + o_h] = H[J+(K + 2) * LDH + o_h] - REFSUM * V[2+M22 * LDV + o_v];
759: }
760: // *
761: if (ACCUM)
762: {
763: KMS = K - INCOL;
764: for (J = Math.Max(1, KTOP - INCOL); J <= KDU; J++)
765: {
766: REFSUM = V[1+M22 * LDV + o_v] * (U[J+(KMS + 1) * LDU + o_u] + V[2+M22 * LDV + o_v] * U[J+(KMS + 2) * LDU + o_u]);
767: U[J+(KMS + 1) * LDU + o_u] = U[J+(KMS + 1) * LDU + o_u] - REFSUM;
768: U[J+(KMS + 2) * LDU + o_u] = U[J+(KMS + 2) * LDU + o_u] - REFSUM * V[2+M22 * LDV + o_v];
769: }
770: }
771: else
772: {
773: if (WANTZ)
774: {
775: for (J = ILOZ; J <= IHIZ; J++)
776: {
777: REFSUM = V[1+M22 * LDV + o_v] * (Z[J+(K + 1) * LDZ + o_z] + V[2+M22 * LDV + o_v] * Z[J+(K + 2) * LDZ + o_z]);
778: Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - REFSUM;
779: Z[J+(K + 2) * LDZ + o_z] = Z[J+(K + 2) * LDZ + o_z] - REFSUM * V[2+M22 * LDV + o_v];
780: }
781: }
782: }
783: }
784: // *
785: // * ==== Vigilant deflation check ====
786: // *
787: MSTART = MTOP;
788: if (KRCOL + 3 * (MSTART - 1) < KTOP) MSTART = MSTART + 1;
789: MEND = MBOT;
790: if (BMP22) MEND = MEND + 1;
791: if (KRCOL == KBOT - 2) MEND = MEND + 1;
792: for (M = MSTART; M <= MEND; M++)
793: {
794: K = Math.Min(KBOT - 1, KRCOL + 3 * (M - 1));
795: // *
796: // * ==== The following convergence test requires that
797: // * . the tradition small-compared-to-nearby-diagonals
798: // * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
799: // * . criteria both be satisfied. The latter improves
800: // * . accuracy in some examples. Falling back on an
801: // * . alternate convergence criterion when TST1 or TST2
802: // * . is zero (as done here) is traditional but probably
803: // * . unnecessary. ====
804: // *
805: if (H[K + 1+K * LDH + o_h] != ZERO)
806: {
807: TST1 = Math.Abs(H[K+K * LDH + o_h]) + Math.Abs(H[K + 1+(K + 1) * LDH + o_h]);
808: if (TST1 == ZERO)
809: {
810: if (K >= KTOP + 1) TST1 = TST1 + Math.Abs(H[K+(K - 1) * LDH + o_h]);
811: if (K >= KTOP + 2) TST1 = TST1 + Math.Abs(H[K+(K - 2) * LDH + o_h]);
812: if (K >= KTOP + 3) TST1 = TST1 + Math.Abs(H[K+(K - 3) * LDH + o_h]);
813: if (K <= KBOT - 2) TST1 = TST1 + Math.Abs(H[K + 2+(K + 1) * LDH + o_h]);
814: if (K <= KBOT - 3) TST1 = TST1 + Math.Abs(H[K + 3+(K + 1) * LDH + o_h]);
815: if (K <= KBOT - 4) TST1 = TST1 + Math.Abs(H[K + 4+(K + 1) * LDH + o_h]);
816: }
817: if (Math.Abs(H[K + 1+K * LDH + o_h]) <= Math.Max(SMLNUM, ULP * TST1))
818: {
819: H12 = Math.Max(Math.Abs(H[K + 1+K * LDH + o_h]), Math.Abs(H[K+(K + 1) * LDH + o_h]));
820: H21 = Math.Min(Math.Abs(H[K + 1+K * LDH + o_h]), Math.Abs(H[K+(K + 1) * LDH + o_h]));
821: H11 = Math.Max(Math.Abs(H[K + 1+(K + 1) * LDH + o_h]), Math.Abs(H[K+K * LDH + o_h] - H[K + 1+(K + 1) * LDH + o_h]));
822: H22 = Math.Min(Math.Abs(H[K + 1+(K + 1) * LDH + o_h]), Math.Abs(H[K+K * LDH + o_h] - H[K + 1+(K + 1) * LDH + o_h]));
823: SCL = H11 + H12;
824: TST2 = H22 * (H11 / SCL);
825: // *
826: if (TST2 == ZERO || H21 * (H12 / SCL) <= Math.Max(SMLNUM, ULP * TST2)) H[K + 1+K * LDH + o_h] = ZERO;
827: }
828: }
829: }
830: // *
831: // * ==== Fill in the last row of each bulge. ====
832: // *
833: MEND = Math.Min(NBMPS, (KBOT - KRCOL - 1) / 3);
834: for (M = MTOP; M <= MEND; M++)
835: {
836: K = KRCOL + 3 * (M - 1);
837: REFSUM = V[1+M * LDV + o_v] * V[3+M * LDV + o_v] * H[K + 4+(K + 3) * LDH + o_h];
838: H[K + 4+(K + 1) * LDH + o_h] = - REFSUM;
839: H[K + 4+(K + 2) * LDH + o_h] = - REFSUM * V[2+M * LDV + o_v];
840: H[K + 4+(K + 3) * LDH + o_h] = H[K + 4+(K + 3) * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
841: }
842: // *
843: // * ==== End of near-the-diagonal bulge chase. ====
844: // *
845: }
846: // *
847: // * ==== Use U (if accumulated) to update far-from-diagonal
848: // * . entries in H. If required, use U to update Z as
849: // * . well. ====
850: // *
851: if (ACCUM)
852: {
853: if (WANTT)
854: {
855: JTOP = 1;
856: JBOT = N;
857: }
858: else
859: {
860: JTOP = KTOP;
861: JBOT = KBOT;
862: }
863: if ((!BLK22) || (INCOL < KTOP) || (NDCOL > KBOT) || (NS <= 2))
864: {
865: // *
866: // * ==== Updates not exploiting the 2-by-2 block
867: // * . structure of U. K1 and NU keep track of
868: // * . the location and size of U in the special
869: // * . cases of introducing bulges and chasing
870: // * . bulges off the bottom. In these special
871: // * . cases and in case the number of shifts
872: // * . is NS = 2, there is no 2-by-2 block
873: // * . structure to exploit. ====
874: // *
875: K1 = Math.Max(1, KTOP - INCOL);
876: NU = (KDU - Math.Max(0, NDCOL - KBOT)) - K1 + 1;
877: // *
878: // * ==== Horizontal Multiply ====
879: // *
880: for (JCOL = Math.Min(NDCOL, KBOT) + 1; (NH >= 0) ? (JCOL <= JBOT) : (JCOL >= JBOT); JCOL += NH)
881: {
882: JLEN = Math.Min(NH, JBOT - JCOL + 1);
883: this._dgemm.Run("C", "N", NU, JLEN, NU, ONE
884: , U, K1+K1 * LDU + o_u, LDU, H, INCOL + K1+JCOL * LDH + o_h, LDH, ZERO, ref WH, offset_wh
885: , LDWH);
886: this._dlacpy.Run("ALL", NU, JLEN, WH, offset_wh, LDWH, ref H, INCOL + K1+JCOL * LDH + o_h
887: , LDH);
888: }
889: // *
890: // * ==== Vertical multiply ====
891: // *
892: for (JROW = JTOP; (NV >= 0) ? (JROW <= Math.Max(KTOP, INCOL) - 1) : (JROW >= Math.Max(KTOP, INCOL) - 1); JROW += NV)
893: {
894: JLEN = Math.Min(NV, Math.Max(KTOP, INCOL) - JROW);
895: this._dgemm.Run("N", "N", JLEN, NU, NU, ONE
896: , H, JROW+(INCOL + K1) * LDH + o_h, LDH, U, K1+K1 * LDU + o_u, LDU, ZERO, ref WV, offset_wv
897: , LDWV);
898: this._dlacpy.Run("ALL", JLEN, NU, WV, offset_wv, LDWV, ref H, JROW+(INCOL + K1) * LDH + o_h
899: , LDH);
900: }
901: // *
902: // * ==== Z multiply (also vertical) ====
903: // *
904: if (WANTZ)
905: {
906: for (JROW = ILOZ; (NV >= 0) ? (JROW <= IHIZ) : (JROW >= IHIZ); JROW += NV)
907: {
908: JLEN = Math.Min(NV, IHIZ - JROW + 1);
909: this._dgemm.Run("N", "N", JLEN, NU, NU, ONE
910: , Z, JROW+(INCOL + K1) * LDZ + o_z, LDZ, U, K1+K1 * LDU + o_u, LDU, ZERO, ref WV, offset_wv
911: , LDWV);
912: this._dlacpy.Run("ALL", JLEN, NU, WV, offset_wv, LDWV, ref Z, JROW+(INCOL + K1) * LDZ + o_z
913: , LDZ);
914: }
915: }
916: }
917: else
918: {
919: // *
920: // * ==== Updates exploiting U's 2-by-2 block structure.
921: // * . (I2, I4, J2, J4 are the last rows and columns
922: // * . of the blocks.) ====
923: // *
924: I2 = (KDU + 1) / 2;
925: I4 = KDU;
926: J2 = I4 - I2;
927: J4 = KDU;
928: // *
929: // * ==== KZS and KNZ deal with the band of zeros
930: // * . along the diagonal of one of the triangular
931: // * . blocks. ====
932: // *
933: KZS = (J4 - J2) - (NS + 1);
934: KNZ = NS + 1;
935: // *
936: // * ==== Horizontal multiply ====
937: // *
938: for (JCOL = Math.Min(NDCOL, KBOT) + 1; (NH >= 0) ? (JCOL <= JBOT) : (JCOL >= JBOT); JCOL += NH)
939: {
940: JLEN = Math.Min(NH, JBOT - JCOL + 1);
941: // *
942: // * ==== Copy bottom of H to top+KZS of scratch ====
943: // * (The first KZS rows get multiplied by zero.) ====
944: // *
945: this._dlacpy.Run("ALL", KNZ, JLEN, H, INCOL + 1 + J2+JCOL * LDH + o_h, LDH, ref WH, KZS + 1+1 * LDWH + o_wh
946: , LDWH);
947: // *
948: // * ==== Multiply by U21' ====
949: // *
950: this._dlaset.Run("ALL", KZS, JLEN, ZERO, ZERO, ref WH, offset_wh
951: , LDWH);
952: this._dtrmm.Run("L", "U", "C", "N", KNZ, JLEN
953: , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WH, KZS + 1+1 * LDWH + o_wh, LDWH);
954: // *
955: // * ==== Multiply top of H by U11' ====
956: // *
957: this._dgemm.Run("C", "N", I2, JLEN, J2, ONE
958: , U, offset_u, LDU, H, INCOL + 1+JCOL * LDH + o_h, LDH, ONE, ref WH, offset_wh
959: , LDWH);
960: // *
961: // * ==== Copy top of H bottom of WH ====
962: // *
963: this._dlacpy.Run("ALL", J2, JLEN, H, INCOL + 1+JCOL * LDH + o_h, LDH, ref WH, I2 + 1+1 * LDWH + o_wh
964: , LDWH);
965: // *
966: // * ==== Multiply by U21' ====
967: // *
968: this._dtrmm.Run("L", "L", "C", "N", J2, JLEN
969: , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WH, I2 + 1+1 * LDWH + o_wh, LDWH);
970: // *
971: // * ==== Multiply by U22 ====
972: // *
973: this._dgemm.Run("C", "N", I4 - I2, JLEN, J4 - J2, ONE
974: , U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, H, INCOL + 1 + J2+JCOL * LDH + o_h, LDH, ONE, ref WH, I2 + 1+1 * LDWH + o_wh
975: , LDWH);
976: // *
977: // * ==== Copy it back ====
978: // *
979: this._dlacpy.Run("ALL", KDU, JLEN, WH, offset_wh, LDWH, ref H, INCOL + 1+JCOL * LDH + o_h
980: , LDH);
981: }
982: // *
983: // * ==== Vertical multiply ====
984: // *
985: for (JROW = JTOP; (NV >= 0) ? (JROW <= Math.Max(INCOL, KTOP) - 1) : (JROW >= Math.Max(INCOL, KTOP) - 1); JROW += NV)
986: {
987: JLEN = Math.Min(NV, Math.Max(INCOL, KTOP) - JROW);
988: // *
989: // * ==== Copy right of H to scratch (the first KZS
990: // * . columns get multiplied by zero) ====
991: // *
992: this._dlacpy.Run("ALL", JLEN, KNZ, H, JROW+(INCOL + 1 + J2) * LDH + o_h, LDH, ref WV, 1+(1 + KZS) * LDWV + o_wv
993: , LDWV);
994: // *
995: // * ==== Multiply by U21 ====
996: // *
997: this._dlaset.Run("ALL", JLEN, KZS, ZERO, ZERO, ref WV, offset_wv
998: , LDWV);
999: this._dtrmm.Run("R", "U", "N", "N", JLEN, KNZ
1000: , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WV, 1+(1 + KZS) * LDWV + o_wv, LDWV);
1001: // *
1002: // * ==== Multiply by U11 ====
1003: // *
1004: this._dgemm.Run("N", "N", JLEN, I2, J2, ONE
1005: , H, JROW+(INCOL + 1) * LDH + o_h, LDH, U, offset_u, LDU, ONE, ref WV, offset_wv
1006: , LDWV);
1007: // *
1008: // * ==== Copy left of H to right of scratch ====
1009: // *
1010: this._dlacpy.Run("ALL", JLEN, J2, H, JROW+(INCOL + 1) * LDH + o_h, LDH, ref WV, 1+(1 + I2) * LDWV + o_wv
1011: , LDWV);
1012: // *
1013: // * ==== Multiply by U21 ====
1014: // *
1015: this._dtrmm.Run("R", "L", "N", "N", JLEN, I4 - I2
1016: , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WV, 1+(1 + I2) * LDWV + o_wv, LDWV);
1017: // *
1018: // * ==== Multiply by U22 ====
1019: // *
1020: this._dgemm.Run("N", "N", JLEN, I4 - I2, J4 - J2, ONE
1021: , H, JROW+(INCOL + 1 + J2) * LDH + o_h, LDH, U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, ONE, ref WV, 1+(1 + I2) * LDWV + o_wv
1022: , LDWV);
1023: // *
1024: // * ==== Copy it back ====
1025: // *
1026: this._dlacpy.Run("ALL", JLEN, KDU, WV, offset_wv, LDWV, ref H, JROW+(INCOL + 1) * LDH + o_h
1027: , LDH);
1028: }
1029: // *
1030: // * ==== Multiply Z (also vertical) ====
1031: // *
1032: if (WANTZ)
1033: {
1034: for (JROW = ILOZ; (NV >= 0) ? (JROW <= IHIZ) : (JROW >= IHIZ); JROW += NV)
1035: {
1036: JLEN = Math.Min(NV, IHIZ - JROW + 1);
1037: // *
1038: // * ==== Copy right of Z to left of scratch (first
1039: // * . KZS columns get multiplied by zero) ====
1040: // *
1041: this._dlacpy.Run("ALL", JLEN, KNZ, Z, JROW+(INCOL + 1 + J2) * LDZ + o_z, LDZ, ref WV, 1+(1 + KZS) * LDWV + o_wv
1042: , LDWV);
1043: // *
1044: // * ==== Multiply by U12 ====
1045: // *
1046: this._dlaset.Run("ALL", JLEN, KZS, ZERO, ZERO, ref WV, offset_wv
1047: , LDWV);
1048: this._dtrmm.Run("R", "U", "N", "N", JLEN, KNZ
1049: , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WV, 1+(1 + KZS) * LDWV + o_wv, LDWV);
1050: // *
1051: // * ==== Multiply by U11 ====
1052: // *
1053: this._dgemm.Run("N", "N", JLEN, I2, J2, ONE
1054: , Z, JROW+(INCOL + 1) * LDZ + o_z, LDZ, U, offset_u, LDU, ONE, ref WV, offset_wv
1055: , LDWV);
1056: // *
1057: // * ==== Copy left of Z to right of scratch ====
1058: // *
1059: this._dlacpy.Run("ALL", JLEN, J2, Z, JROW+(INCOL + 1) * LDZ + o_z, LDZ, ref WV, 1+(1 + I2) * LDWV + o_wv
1060: , LDWV);
1061: // *
1062: // * ==== Multiply by U21 ====
1063: // *
1064: this._dtrmm.Run("R", "L", "N", "N", JLEN, I4 - I2
1065: , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WV, 1+(1 + I2) * LDWV + o_wv, LDWV);
1066: // *
1067: // * ==== Multiply by U22 ====
1068: // *
1069: this._dgemm.Run("N", "N", JLEN, I4 - I2, J4 - J2, ONE
1070: , Z, JROW+(INCOL + 1 + J2) * LDZ + o_z, LDZ, U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, ONE, ref WV, 1+(1 + I2) * LDWV + o_wv
1071: , LDWV);
1072: // *
1073: // * ==== Copy the result back to Z ====
1074: // *
1075: this._dlacpy.Run("ALL", JLEN, KDU, WV, offset_wv, LDWV, ref Z, JROW+(INCOL + 1) * LDZ + o_z
1076: , LDZ);
1077: }
1078: }
1079: }
1080: }
1081: }
1082: // *
1083: // * ==== End of DLAQR5 ====
1084: // *
1085:
1086: #endregion
1087:
1088: }
1089: }
1090: }