Annotation of rpl/lapack/lapack/dlaqr5.f, revision 1.13
1.13 ! bertrand 1: *> \brief \b DLAQR5 performs a single small-bulge multi-shift QR sweep.
1.10 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAQR5 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr5.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr5.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr5.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
22: * SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U,
23: * LDU, NV, WV, LDWV, NH, WH, LDWH )
24: *
25: * .. Scalar Arguments ..
26: * INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
27: * $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
28: * LOGICAL WANTT, WANTZ
29: * ..
30: * .. Array Arguments ..
31: * DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ),
32: * $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ),
33: * $ Z( LDZ, * )
34: * ..
35: *
36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> DLAQR5, called by DLAQR0, performs a
43: *> single small-bulge multi-shift QR sweep.
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] WANTT
50: *> \verbatim
51: *> WANTT is logical scalar
52: *> WANTT = .true. if the quasi-triangular Schur factor
53: *> is being computed. WANTT is set to .false. otherwise.
54: *> \endverbatim
55: *>
56: *> \param[in] WANTZ
57: *> \verbatim
58: *> WANTZ is logical scalar
59: *> WANTZ = .true. if the orthogonal Schur factor is being
60: *> computed. WANTZ is set to .false. otherwise.
61: *> \endverbatim
62: *>
63: *> \param[in] KACC22
64: *> \verbatim
65: *> KACC22 is integer with value 0, 1, or 2.
66: *> Specifies the computation mode of far-from-diagonal
67: *> orthogonal updates.
68: *> = 0: DLAQR5 does not accumulate reflections and does not
69: *> use matrix-matrix multiply to update far-from-diagonal
70: *> matrix entries.
71: *> = 1: DLAQR5 accumulates reflections and uses matrix-matrix
72: *> multiply to update the far-from-diagonal matrix entries.
73: *> = 2: DLAQR5 accumulates reflections, uses matrix-matrix
74: *> multiply to update the far-from-diagonal matrix entries,
75: *> and takes advantage of 2-by-2 block structure during
76: *> matrix multiplies.
77: *> \endverbatim
78: *>
79: *> \param[in] N
80: *> \verbatim
81: *> N is integer scalar
82: *> N is the order of the Hessenberg matrix H upon which this
83: *> subroutine operates.
84: *> \endverbatim
85: *>
86: *> \param[in] KTOP
87: *> \verbatim
88: *> KTOP is integer scalar
89: *> \endverbatim
90: *>
91: *> \param[in] KBOT
92: *> \verbatim
93: *> KBOT is integer scalar
94: *> These are the first and last rows and columns of an
95: *> isolated diagonal block upon which the QR sweep is to be
96: *> applied. It is assumed without a check that
97: *> either KTOP = 1 or H(KTOP,KTOP-1) = 0
98: *> and
99: *> either KBOT = N or H(KBOT+1,KBOT) = 0.
100: *> \endverbatim
101: *>
102: *> \param[in] NSHFTS
103: *> \verbatim
104: *> NSHFTS is integer scalar
105: *> NSHFTS gives the number of simultaneous shifts. NSHFTS
106: *> must be positive and even.
107: *> \endverbatim
108: *>
109: *> \param[in,out] SR
110: *> \verbatim
111: *> SR is DOUBLE PRECISION array of size (NSHFTS)
112: *> \endverbatim
113: *>
114: *> \param[in,out] SI
115: *> \verbatim
116: *> SI is DOUBLE PRECISION array of size (NSHFTS)
117: *> SR contains the real parts and SI contains the imaginary
118: *> parts of the NSHFTS shifts of origin that define the
119: *> multi-shift QR sweep. On output SR and SI may be
120: *> reordered.
121: *> \endverbatim
122: *>
123: *> \param[in,out] H
124: *> \verbatim
125: *> H is DOUBLE PRECISION array of size (LDH,N)
126: *> On input H contains a Hessenberg matrix. On output a
127: *> multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
128: *> to the isolated diagonal block in rows and columns KTOP
129: *> through KBOT.
130: *> \endverbatim
131: *>
132: *> \param[in] LDH
133: *> \verbatim
134: *> LDH is integer scalar
135: *> LDH is the leading dimension of H just as declared in the
136: *> calling procedure. LDH.GE.MAX(1,N).
137: *> \endverbatim
138: *>
139: *> \param[in] ILOZ
140: *> \verbatim
141: *> ILOZ is INTEGER
142: *> \endverbatim
143: *>
144: *> \param[in] IHIZ
145: *> \verbatim
146: *> IHIZ is INTEGER
147: *> Specify the rows of Z to which transformations must be
148: *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
149: *> \endverbatim
150: *>
151: *> \param[in,out] Z
152: *> \verbatim
153: *> Z is DOUBLE PRECISION array of size (LDZ,IHI)
154: *> If WANTZ = .TRUE., then the QR Sweep orthogonal
155: *> similarity transformation is accumulated into
156: *> Z(ILOZ:IHIZ,ILO:IHI) from the right.
157: *> If WANTZ = .FALSE., then Z is unreferenced.
158: *> \endverbatim
159: *>
160: *> \param[in] LDZ
161: *> \verbatim
162: *> LDZ is integer scalar
163: *> LDA is the leading dimension of Z just as declared in
164: *> the calling procedure. LDZ.GE.N.
165: *> \endverbatim
166: *>
167: *> \param[out] V
168: *> \verbatim
169: *> V is DOUBLE PRECISION array of size (LDV,NSHFTS/2)
170: *> \endverbatim
171: *>
172: *> \param[in] LDV
173: *> \verbatim
174: *> LDV is integer scalar
175: *> LDV is the leading dimension of V as declared in the
176: *> calling procedure. LDV.GE.3.
177: *> \endverbatim
178: *>
179: *> \param[out] U
180: *> \verbatim
181: *> U is DOUBLE PRECISION array of size
182: *> (LDU,3*NSHFTS-3)
183: *> \endverbatim
184: *>
185: *> \param[in] LDU
186: *> \verbatim
187: *> LDU is integer scalar
188: *> LDU is the leading dimension of U just as declared in the
189: *> in the calling subroutine. LDU.GE.3*NSHFTS-3.
190: *> \endverbatim
191: *>
192: *> \param[in] NH
193: *> \verbatim
194: *> NH is integer scalar
195: *> NH is the number of columns in array WH available for
196: *> workspace. NH.GE.1.
197: *> \endverbatim
198: *>
199: *> \param[out] WH
200: *> \verbatim
201: *> WH is DOUBLE PRECISION array of size (LDWH,NH)
202: *> \endverbatim
203: *>
204: *> \param[in] LDWH
205: *> \verbatim
206: *> LDWH is integer scalar
207: *> Leading dimension of WH just as declared in the
208: *> calling procedure. LDWH.GE.3*NSHFTS-3.
209: *> \endverbatim
210: *>
211: *> \param[in] NV
212: *> \verbatim
213: *> NV is integer scalar
214: *> NV is the number of rows in WV agailable for workspace.
215: *> NV.GE.1.
216: *> \endverbatim
217: *>
218: *> \param[out] WV
219: *> \verbatim
220: *> WV is DOUBLE PRECISION array of size
221: *> (LDWV,3*NSHFTS-3)
222: *> \endverbatim
223: *>
224: *> \param[in] LDWV
225: *> \verbatim
226: *> LDWV is integer scalar
227: *> LDWV is the leading dimension of WV as declared in the
228: *> in the calling subroutine. LDWV.GE.NV.
229: *> \endverbatim
230: *
231: * Authors:
232: * ========
233: *
234: *> \author Univ. of Tennessee
235: *> \author Univ. of California Berkeley
236: *> \author Univ. of Colorado Denver
237: *> \author NAG Ltd.
238: *
1.13 ! bertrand 239: *> \date September 2012
1.10 bertrand 240: *
241: *> \ingroup doubleOTHERauxiliary
242: *
243: *> \par Contributors:
244: * ==================
245: *>
246: *> Karen Braman and Ralph Byers, Department of Mathematics,
247: *> University of Kansas, USA
248: *
249: *> \par References:
250: * ================
251: *>
252: *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
253: *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
254: *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
255: *> 929--947, 2002.
256: *>
257: * =====================================================================
1.1 bertrand 258: SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
259: $ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U,
260: $ LDU, NV, WV, LDWV, NH, WH, LDWH )
261: *
1.13 ! bertrand 262: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 263: * -- LAPACK is a software package provided by Univ. of Tennessee, --
264: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13 ! bertrand 265: * September 2012
1.1 bertrand 266: *
267: * .. Scalar Arguments ..
268: INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
269: $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
270: LOGICAL WANTT, WANTZ
271: * ..
272: * .. Array Arguments ..
273: DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ),
274: $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ),
275: $ Z( LDZ, * )
276: * ..
277: *
1.10 bertrand 278: * ================================================================
1.1 bertrand 279: * .. Parameters ..
280: DOUBLE PRECISION ZERO, ONE
281: PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
282: * ..
283: * .. Local Scalars ..
284: DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM,
285: $ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, TST1, TST2,
286: $ ULP
287: INTEGER I, I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
288: $ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
289: $ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
290: $ NS, NU
291: LOGICAL ACCUM, BLK22, BMP22
292: * ..
293: * .. External Functions ..
294: DOUBLE PRECISION DLAMCH
295: EXTERNAL DLAMCH
296: * ..
297: * .. Intrinsic Functions ..
298: *
299: INTRINSIC ABS, DBLE, MAX, MIN, MOD
300: * ..
301: * .. Local Arrays ..
302: DOUBLE PRECISION VT( 3 )
303: * ..
304: * .. External Subroutines ..
305: EXTERNAL DGEMM, DLABAD, DLACPY, DLAQR1, DLARFG, DLASET,
306: $ DTRMM
307: * ..
308: * .. Executable Statements ..
309: *
310: * ==== If there are no shifts, then there is nothing to do. ====
311: *
312: IF( NSHFTS.LT.2 )
313: $ RETURN
314: *
315: * ==== If the active block is empty or 1-by-1, then there
316: * . is nothing to do. ====
317: *
318: IF( KTOP.GE.KBOT )
319: $ RETURN
320: *
321: * ==== Shuffle shifts into pairs of real shifts and pairs
322: * . of complex conjugate shifts assuming complex
323: * . conjugate shifts are already adjacent to one
324: * . another. ====
325: *
326: DO 10 I = 1, NSHFTS - 2, 2
327: IF( SI( I ).NE.-SI( I+1 ) ) THEN
328: *
329: SWAP = SR( I )
330: SR( I ) = SR( I+1 )
331: SR( I+1 ) = SR( I+2 )
332: SR( I+2 ) = SWAP
333: *
334: SWAP = SI( I )
335: SI( I ) = SI( I+1 )
336: SI( I+1 ) = SI( I+2 )
337: SI( I+2 ) = SWAP
338: END IF
339: 10 CONTINUE
340: *
341: * ==== NSHFTS is supposed to be even, but if it is odd,
342: * . then simply reduce it by one. The shuffle above
343: * . ensures that the dropped shift is real and that
344: * . the remaining shifts are paired. ====
345: *
346: NS = NSHFTS - MOD( NSHFTS, 2 )
347: *
348: * ==== Machine constants for deflation ====
349: *
350: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
351: SAFMAX = ONE / SAFMIN
352: CALL DLABAD( SAFMIN, SAFMAX )
353: ULP = DLAMCH( 'PRECISION' )
354: SMLNUM = SAFMIN*( DBLE( N ) / ULP )
355: *
356: * ==== Use accumulated reflections to update far-from-diagonal
357: * . entries ? ====
358: *
359: ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
360: *
361: * ==== If so, exploit the 2-by-2 block structure? ====
362: *
363: BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
364: *
365: * ==== clear trash ====
366: *
367: IF( KTOP+2.LE.KBOT )
368: $ H( KTOP+2, KTOP ) = ZERO
369: *
370: * ==== NBMPS = number of 2-shift bulges in the chain ====
371: *
372: NBMPS = NS / 2
373: *
374: * ==== KDU = width of slab ====
375: *
376: KDU = 6*NBMPS - 3
377: *
378: * ==== Create and chase chains of NBMPS bulges ====
379: *
380: DO 220 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
381: NDCOL = INCOL + KDU
382: IF( ACCUM )
383: $ CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )
384: *
385: * ==== Near-the-diagonal bulge chase. The following loop
386: * . performs the near-the-diagonal part of a small bulge
387: * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
388: * . chunk extends from column INCOL to column NDCOL
389: * . (including both column INCOL and column NDCOL). The
390: * . following loop chases a 3*NBMPS column long chain of
391: * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
392: * . may be less than KTOP and and NDCOL may be greater than
393: * . KBOT indicating phantom columns from which to chase
394: * . bulges before they are actually introduced or to which
395: * . to chase bulges beyond column KBOT.) ====
396: *
397: DO 150 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
398: *
399: * ==== Bulges number MTOP to MBOT are active double implicit
400: * . shift bulges. There may or may not also be small
401: * . 2-by-2 bulge, if there is room. The inactive bulges
402: * . (if any) must wait until the active bulges have moved
403: * . down the diagonal to make room. The phantom matrix
404: * . paradigm described above helps keep track. ====
405: *
406: MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
407: MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
408: M22 = MBOT + 1
409: BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
410: $ ( KBOT-2 )
411: *
412: * ==== Generate reflections to chase the chain right
413: * . one column. (The minimum value of K is KTOP-1.) ====
414: *
415: DO 20 M = MTOP, MBOT
416: K = KRCOL + 3*( M-1 )
417: IF( K.EQ.KTOP-1 ) THEN
418: CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ),
419: $ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ),
420: $ V( 1, M ) )
421: ALPHA = V( 1, M )
422: CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
423: ELSE
424: BETA = H( K+1, K )
425: V( 2, M ) = H( K+2, K )
426: V( 3, M ) = H( K+3, K )
427: CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )
428: *
429: * ==== A Bulge may collapse because of vigilant
430: * . deflation or destructive underflow. In the
431: * . underflow case, try the two-small-subdiagonals
432: * . trick to try to reinflate the bulge. ====
433: *
434: IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE.
435: $ ZERO .OR. H( K+3, K+2 ).EQ.ZERO ) THEN
436: *
437: * ==== Typical case: not collapsed (yet). ====
438: *
439: H( K+1, K ) = BETA
440: H( K+2, K ) = ZERO
441: H( K+3, K ) = ZERO
442: ELSE
443: *
444: * ==== Atypical case: collapsed. Attempt to
445: * . reintroduce ignoring H(K+1,K) and H(K+2,K).
446: * . If the fill resulting from the new
447: * . reflector is too large, then abandon it.
448: * . Otherwise, use the new one. ====
449: *
450: CALL DLAQR1( 3, H( K+1, K+1 ), LDH, SR( 2*M-1 ),
451: $ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ),
452: $ VT )
453: ALPHA = VT( 1 )
454: CALL DLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
455: REFSUM = VT( 1 )*( H( K+1, K )+VT( 2 )*
456: $ H( K+2, K ) )
457: *
458: IF( ABS( H( K+2, K )-REFSUM*VT( 2 ) )+
459: $ ABS( REFSUM*VT( 3 ) ).GT.ULP*
460: $ ( ABS( H( K, K ) )+ABS( H( K+1,
461: $ K+1 ) )+ABS( H( K+2, K+2 ) ) ) ) THEN
462: *
463: * ==== Starting a new bulge here would
464: * . create non-negligible fill. Use
465: * . the old one with trepidation. ====
466: *
467: H( K+1, K ) = BETA
468: H( K+2, K ) = ZERO
469: H( K+3, K ) = ZERO
470: ELSE
471: *
472: * ==== Stating a new bulge here would
473: * . create only negligible fill.
474: * . Replace the old reflector with
475: * . the new one. ====
476: *
477: H( K+1, K ) = H( K+1, K ) - REFSUM
478: H( K+2, K ) = ZERO
479: H( K+3, K ) = ZERO
480: V( 1, M ) = VT( 1 )
481: V( 2, M ) = VT( 2 )
482: V( 3, M ) = VT( 3 )
483: END IF
484: END IF
485: END IF
486: 20 CONTINUE
487: *
488: * ==== Generate a 2-by-2 reflection, if needed. ====
489: *
490: K = KRCOL + 3*( M22-1 )
491: IF( BMP22 ) THEN
492: IF( K.EQ.KTOP-1 ) THEN
493: CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ),
494: $ SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ),
495: $ V( 1, M22 ) )
496: BETA = V( 1, M22 )
497: CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
498: ELSE
499: BETA = H( K+1, K )
500: V( 2, M22 ) = H( K+2, K )
501: CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
502: H( K+1, K ) = BETA
503: H( K+2, K ) = ZERO
504: END IF
505: END IF
506: *
507: * ==== Multiply H by reflections from the left ====
508: *
509: IF( ACCUM ) THEN
510: JBOT = MIN( NDCOL, KBOT )
511: ELSE IF( WANTT ) THEN
512: JBOT = N
513: ELSE
514: JBOT = KBOT
515: END IF
516: DO 40 J = MAX( KTOP, KRCOL ), JBOT
517: MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
518: DO 30 M = MTOP, MEND
519: K = KRCOL + 3*( M-1 )
520: REFSUM = V( 1, M )*( H( K+1, J )+V( 2, M )*
521: $ H( K+2, J )+V( 3, M )*H( K+3, J ) )
522: H( K+1, J ) = H( K+1, J ) - REFSUM
523: H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
524: H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
525: 30 CONTINUE
526: 40 CONTINUE
527: IF( BMP22 ) THEN
528: K = KRCOL + 3*( M22-1 )
529: DO 50 J = MAX( K+1, KTOP ), JBOT
530: REFSUM = V( 1, M22 )*( H( K+1, J )+V( 2, M22 )*
531: $ H( K+2, J ) )
532: H( K+1, J ) = H( K+1, J ) - REFSUM
533: H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
534: 50 CONTINUE
535: END IF
536: *
537: * ==== Multiply H by reflections from the right.
538: * . Delay filling in the last row until the
539: * . vigilant deflation check is complete. ====
540: *
541: IF( ACCUM ) THEN
542: JTOP = MAX( KTOP, INCOL )
543: ELSE IF( WANTT ) THEN
544: JTOP = 1
545: ELSE
546: JTOP = KTOP
547: END IF
548: DO 90 M = MTOP, MBOT
549: IF( V( 1, M ).NE.ZERO ) THEN
550: K = KRCOL + 3*( M-1 )
551: DO 60 J = JTOP, MIN( KBOT, K+3 )
552: REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
553: $ H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
554: H( J, K+1 ) = H( J, K+1 ) - REFSUM
555: H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M )
556: H( J, K+3 ) = H( J, K+3 ) - REFSUM*V( 3, M )
557: 60 CONTINUE
558: *
559: IF( ACCUM ) THEN
560: *
561: * ==== Accumulate U. (If necessary, update Z later
562: * . with with an efficient matrix-matrix
563: * . multiply.) ====
564: *
565: KMS = K - INCOL
566: DO 70 J = MAX( 1, KTOP-INCOL ), KDU
567: REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
568: $ U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
569: U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
570: U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M )
571: U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*V( 3, M )
572: 70 CONTINUE
573: ELSE IF( WANTZ ) THEN
574: *
575: * ==== U is not accumulated, so update Z
576: * . now by multiplying by reflections
577: * . from the right. ====
578: *
579: DO 80 J = ILOZ, IHIZ
580: REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
581: $ Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
582: Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
583: Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M )
584: Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M )
585: 80 CONTINUE
586: END IF
587: END IF
588: 90 CONTINUE
589: *
590: * ==== Special case: 2-by-2 reflection (if needed) ====
591: *
592: K = KRCOL + 3*( M22-1 )
1.7 bertrand 593: IF( BMP22 ) THEN
594: IF ( V( 1, M22 ).NE.ZERO ) THEN
595: DO 100 J = JTOP, MIN( KBOT, K+3 )
596: REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
597: $ H( J, K+2 ) )
598: H( J, K+1 ) = H( J, K+1 ) - REFSUM
599: H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 )
600: 100 CONTINUE
601: *
602: IF( ACCUM ) THEN
603: KMS = K - INCOL
604: DO 110 J = MAX( 1, KTOP-INCOL ), KDU
605: REFSUM = V( 1, M22 )*( U( J, KMS+1 )+
606: $ V( 2, M22 )*U( J, KMS+2 ) )
607: U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
608: U( J, KMS+2 ) = U( J, KMS+2 ) -
609: $ REFSUM*V( 2, M22 )
1.1 bertrand 610: 110 CONTINUE
1.7 bertrand 611: ELSE IF( WANTZ ) THEN
612: DO 120 J = ILOZ, IHIZ
613: REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
614: $ Z( J, K+2 ) )
615: Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
616: Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 )
617: 120 CONTINUE
618: END IF
1.1 bertrand 619: END IF
620: END IF
621: *
622: * ==== Vigilant deflation check ====
623: *
624: MSTART = MTOP
625: IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
626: $ MSTART = MSTART + 1
627: MEND = MBOT
628: IF( BMP22 )
629: $ MEND = MEND + 1
630: IF( KRCOL.EQ.KBOT-2 )
631: $ MEND = MEND + 1
632: DO 130 M = MSTART, MEND
633: K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
634: *
635: * ==== The following convergence test requires that
636: * . the tradition small-compared-to-nearby-diagonals
637: * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
638: * . criteria both be satisfied. The latter improves
639: * . accuracy in some examples. Falling back on an
640: * . alternate convergence criterion when TST1 or TST2
641: * . is zero (as done here) is traditional but probably
642: * . unnecessary. ====
643: *
644: IF( H( K+1, K ).NE.ZERO ) THEN
645: TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) )
646: IF( TST1.EQ.ZERO ) THEN
647: IF( K.GE.KTOP+1 )
648: $ TST1 = TST1 + ABS( H( K, K-1 ) )
649: IF( K.GE.KTOP+2 )
650: $ TST1 = TST1 + ABS( H( K, K-2 ) )
651: IF( K.GE.KTOP+3 )
652: $ TST1 = TST1 + ABS( H( K, K-3 ) )
653: IF( K.LE.KBOT-2 )
654: $ TST1 = TST1 + ABS( H( K+2, K+1 ) )
655: IF( K.LE.KBOT-3 )
656: $ TST1 = TST1 + ABS( H( K+3, K+1 ) )
657: IF( K.LE.KBOT-4 )
658: $ TST1 = TST1 + ABS( H( K+4, K+1 ) )
659: END IF
660: IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
661: $ THEN
662: H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )
663: H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )
664: H11 = MAX( ABS( H( K+1, K+1 ) ),
665: $ ABS( H( K, K )-H( K+1, K+1 ) ) )
666: H22 = MIN( ABS( H( K+1, K+1 ) ),
667: $ ABS( H( K, K )-H( K+1, K+1 ) ) )
668: SCL = H11 + H12
669: TST2 = H22*( H11 / SCL )
670: *
671: IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE.
672: $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
673: END IF
674: END IF
675: 130 CONTINUE
676: *
677: * ==== Fill in the last row of each bulge. ====
678: *
679: MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
680: DO 140 M = MTOP, MEND
681: K = KRCOL + 3*( M-1 )
682: REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
683: H( K+4, K+1 ) = -REFSUM
684: H( K+4, K+2 ) = -REFSUM*V( 2, M )
685: H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M )
686: 140 CONTINUE
687: *
688: * ==== End of near-the-diagonal bulge chase. ====
689: *
690: 150 CONTINUE
691: *
692: * ==== Use U (if accumulated) to update far-from-diagonal
693: * . entries in H. If required, use U to update Z as
694: * . well. ====
695: *
696: IF( ACCUM ) THEN
697: IF( WANTT ) THEN
698: JTOP = 1
699: JBOT = N
700: ELSE
701: JTOP = KTOP
702: JBOT = KBOT
703: END IF
704: IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
705: $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
706: *
707: * ==== Updates not exploiting the 2-by-2 block
708: * . structure of U. K1 and NU keep track of
709: * . the location and size of U in the special
710: * . cases of introducing bulges and chasing
711: * . bulges off the bottom. In these special
712: * . cases and in case the number of shifts
713: * . is NS = 2, there is no 2-by-2 block
714: * . structure to exploit. ====
715: *
716: K1 = MAX( 1, KTOP-INCOL )
717: NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
718: *
719: * ==== Horizontal Multiply ====
720: *
721: DO 160 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
722: JLEN = MIN( NH, JBOT-JCOL+1 )
723: CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
724: $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
725: $ LDWH )
726: CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH,
727: $ H( INCOL+K1, JCOL ), LDH )
728: 160 CONTINUE
729: *
730: * ==== Vertical multiply ====
731: *
732: DO 170 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
733: JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
734: CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE,
735: $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
736: $ LDU, ZERO, WV, LDWV )
737: CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV,
738: $ H( JROW, INCOL+K1 ), LDH )
739: 170 CONTINUE
740: *
741: * ==== Z multiply (also vertical) ====
742: *
743: IF( WANTZ ) THEN
744: DO 180 JROW = ILOZ, IHIZ, NV
745: JLEN = MIN( NV, IHIZ-JROW+1 )
746: CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE,
747: $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
748: $ LDU, ZERO, WV, LDWV )
749: CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV,
750: $ Z( JROW, INCOL+K1 ), LDZ )
751: 180 CONTINUE
752: END IF
753: ELSE
754: *
755: * ==== Updates exploiting U's 2-by-2 block structure.
756: * . (I2, I4, J2, J4 are the last rows and columns
757: * . of the blocks.) ====
758: *
759: I2 = ( KDU+1 ) / 2
760: I4 = KDU
761: J2 = I4 - I2
762: J4 = KDU
763: *
764: * ==== KZS and KNZ deal with the band of zeros
765: * . along the diagonal of one of the triangular
766: * . blocks. ====
767: *
768: KZS = ( J4-J2 ) - ( NS+1 )
769: KNZ = NS + 1
770: *
771: * ==== Horizontal multiply ====
772: *
773: DO 190 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
774: JLEN = MIN( NH, JBOT-JCOL+1 )
775: *
776: * ==== Copy bottom of H to top+KZS of scratch ====
777: * (The first KZS rows get multiplied by zero.) ====
778: *
779: CALL DLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
780: $ LDH, WH( KZS+1, 1 ), LDWH )
781: *
1.9 bertrand 782: * ==== Multiply by U21**T ====
1.1 bertrand 783: *
784: CALL DLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
785: CALL DTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
786: $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
787: $ LDWH )
788: *
1.9 bertrand 789: * ==== Multiply top of H by U11**T ====
1.1 bertrand 790: *
791: CALL DGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
792: $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
793: *
794: * ==== Copy top of H to bottom of WH ====
795: *
796: CALL DLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
797: $ WH( I2+1, 1 ), LDWH )
798: *
1.9 bertrand 799: * ==== Multiply by U21**T ====
1.1 bertrand 800: *
801: CALL DTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
802: $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
803: *
804: * ==== Multiply by U22 ====
805: *
806: CALL DGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
807: $ U( J2+1, I2+1 ), LDU,
808: $ H( INCOL+1+J2, JCOL ), LDH, ONE,
809: $ WH( I2+1, 1 ), LDWH )
810: *
811: * ==== Copy it back ====
812: *
813: CALL DLACPY( 'ALL', KDU, JLEN, WH, LDWH,
814: $ H( INCOL+1, JCOL ), LDH )
815: 190 CONTINUE
816: *
817: * ==== Vertical multiply ====
818: *
819: DO 200 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
820: JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
821: *
822: * ==== Copy right of H to scratch (the first KZS
823: * . columns get multiplied by zero) ====
824: *
825: CALL DLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
826: $ LDH, WV( 1, 1+KZS ), LDWV )
827: *
828: * ==== Multiply by U21 ====
829: *
830: CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
831: CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
832: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
833: $ LDWV )
834: *
835: * ==== Multiply by U11 ====
836: *
837: CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE,
838: $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
839: $ LDWV )
840: *
841: * ==== Copy left of H to right of scratch ====
842: *
843: CALL DLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
844: $ WV( 1, 1+I2 ), LDWV )
845: *
846: * ==== Multiply by U21 ====
847: *
848: CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
849: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
850: *
851: * ==== Multiply by U22 ====
852: *
853: CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
854: $ H( JROW, INCOL+1+J2 ), LDH,
855: $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
856: $ LDWV )
857: *
858: * ==== Copy it back ====
859: *
860: CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV,
861: $ H( JROW, INCOL+1 ), LDH )
862: 200 CONTINUE
863: *
864: * ==== Multiply Z (also vertical) ====
865: *
866: IF( WANTZ ) THEN
867: DO 210 JROW = ILOZ, IHIZ, NV
868: JLEN = MIN( NV, IHIZ-JROW+1 )
869: *
870: * ==== Copy right of Z to left of scratch (first
871: * . KZS columns get multiplied by zero) ====
872: *
873: CALL DLACPY( 'ALL', JLEN, KNZ,
874: $ Z( JROW, INCOL+1+J2 ), LDZ,
875: $ WV( 1, 1+KZS ), LDWV )
876: *
877: * ==== Multiply by U12 ====
878: *
879: CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
880: $ LDWV )
881: CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
882: $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
883: $ LDWV )
884: *
885: * ==== Multiply by U11 ====
886: *
887: CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE,
888: $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
889: $ WV, LDWV )
890: *
891: * ==== Copy left of Z to right of scratch ====
892: *
893: CALL DLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
894: $ LDZ, WV( 1, 1+I2 ), LDWV )
895: *
896: * ==== Multiply by U21 ====
897: *
898: CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
899: $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
900: $ LDWV )
901: *
902: * ==== Multiply by U22 ====
903: *
904: CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
905: $ Z( JROW, INCOL+1+J2 ), LDZ,
906: $ U( J2+1, I2+1 ), LDU, ONE,
907: $ WV( 1, 1+I2 ), LDWV )
908: *
909: * ==== Copy the result back to Z ====
910: *
911: CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV,
912: $ Z( JROW, INCOL+1 ), LDZ )
913: 210 CONTINUE
914: END IF
915: END IF
916: END IF
917: 220 CONTINUE
918: *
919: * ==== End of DLAQR5 ====
920: *
921: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaqr5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack