LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

◆ dlaqr5()

 subroutine dlaqr5 ( logical wantt, logical wantz, integer kacc22, integer n, integer ktop, integer kbot, integer nshfts, double precision, dimension( * ) sr, double precision, dimension( * ) si, double precision, dimension( ldh, * ) h, integer ldh, integer iloz, integer ihiz, double precision, dimension( ldz, * ) z, integer ldz, double precision, dimension( ldv, * ) v, integer ldv, double precision, dimension( ldu, * ) u, integer ldu, integer nv, double precision, dimension( ldwv, * ) wv, integer ldwv, integer nh, double precision, dimension( ldwh, * ) wh, integer ldwh )

DLAQR5 performs a single small-bulge multi-shift QR sweep.

Purpose:
```    DLAQR5, called by DLAQR0, performs a
single small-bulge multi-shift QR sweep.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL WANTT = .true. if the quasi-triangular Schur factor is being computed. WANTT is set to .false. otherwise.``` [in] WANTZ ``` WANTZ is LOGICAL WANTZ = .true. if the orthogonal Schur factor is being computed. WANTZ is set to .false. otherwise.``` [in] KACC22 ``` KACC22 is INTEGER with value 0, 1, or 2. Specifies the computation mode of far-from-diagonal orthogonal updates. = 0: DLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries. = 1: DLAQR5 accumulates reflections and uses matrix-matrix multiply to update the far-from-diagonal matrix entries. = 2: Same as KACC22 = 1. This option used to enable exploiting the 2-by-2 structure during matrix multiplications, but this is no longer supported.``` [in] N ``` N is INTEGER N is the order of the Hessenberg matrix H upon which this subroutine operates.``` [in] KTOP ` KTOP is INTEGER` [in] KBOT ``` KBOT is INTEGER These are the first and last rows and columns of an isolated diagonal block upon which the QR sweep is to be applied. It is assumed without a check that either KTOP = 1 or H(KTOP,KTOP-1) = 0 and either KBOT = N or H(KBOT+1,KBOT) = 0.``` [in] NSHFTS ``` NSHFTS is INTEGER NSHFTS gives the number of simultaneous shifts. NSHFTS must be positive and even.``` [in,out] SR ` SR is DOUBLE PRECISION array, dimension (NSHFTS)` [in,out] SI ``` SI is DOUBLE PRECISION array, dimension (NSHFTS) SR contains the real parts and SI contains the imaginary parts of the NSHFTS shifts of origin that define the multi-shift QR sweep. On output SR and SI may be reordered.``` [in,out] H ``` H is DOUBLE PRECISION array, dimension (LDH,N) On input H contains a Hessenberg matrix. On output a multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied to the isolated diagonal block in rows and columns KTOP through KBOT.``` [in] LDH ``` LDH is INTEGER LDH is the leading dimension of H just as declared in the calling procedure. LDH >= MAX(1,N).``` [in] ILOZ ` ILOZ is INTEGER` [in] IHIZ ``` IHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N``` [in,out] Z ``` Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ) If WANTZ = .TRUE., then the QR Sweep orthogonal similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ = .FALSE., then Z is unreferenced.``` [in] LDZ ``` LDZ is INTEGER LDA is the leading dimension of Z just as declared in the calling procedure. LDZ >= N.``` [out] V ` V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2)` [in] LDV ``` LDV is INTEGER LDV is the leading dimension of V as declared in the calling procedure. LDV >= 3.``` [out] U ` U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS)` [in] LDU ``` LDU is INTEGER LDU is the leading dimension of U just as declared in the in the calling subroutine. LDU >= 2*NSHFTS.``` [in] NV ``` NV is INTEGER NV is the number of rows in WV agailable for workspace. NV >= 1.``` [out] WV ` WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS)` [in] LDWV ``` LDWV is INTEGER LDWV is the leading dimension of WV as declared in the in the calling subroutine. LDWV >= NV.``` [in] NH ``` NH is INTEGER NH is the number of columns in array WH available for workspace. NH >= 1.``` [out] WH ` WH is DOUBLE PRECISION array, dimension (LDWH,NH)` [in] LDWH ``` LDWH is INTEGER Leading dimension of WH just as declared in the calling procedure. LDWH >= 2*NSHFTS.```
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Lars Karlsson, Daniel Kressner, and Bruno Lang

Thijs Steel, Department of Computer science, KU Leuven, Belgium

References:
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume 23, pages 929–947, 2002.

Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed chains of bulges in multishift QR algorithms. ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014).

Definition at line 262 of file dlaqr5.f.

265 IMPLICIT NONE
266*
267* -- LAPACK auxiliary routine --
268* -- LAPACK is a software package provided by Univ. of Tennessee, --
269* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
270*
271* .. Scalar Arguments ..
272 INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
273 \$ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
274 LOGICAL WANTT, WANTZ
275* ..
276* .. Array Arguments ..
277 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ),
278 \$ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ),
279 \$ Z( LDZ, * )
280* ..
281*
282* ================================================================
283* .. Parameters ..
284 DOUBLE PRECISION ZERO, ONE
285 parameter( zero = 0.0d0, one = 1.0d0 )
286* ..
287* .. Local Scalars ..
288 DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM,
289 \$ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, T1, T2,
290 \$ T3, TST1, TST2, ULP
291 INTEGER I, I2, I4, INCOL, J, JBOT, JCOL, JLEN,
292 \$ JROW, JTOP, K, K1, KDU, KMS, KRCOL,
293 \$ M, M22, MBOT, MTOP, NBMPS, NDCOL,
294 \$ NS, NU
295 LOGICAL ACCUM, BMP22
296* ..
297* .. External Functions ..
298 DOUBLE PRECISION DLAMCH
299 EXTERNAL dlamch
300* ..
301* .. Intrinsic Functions ..
302*
303 INTRINSIC abs, dble, max, min, mod
304* ..
305* .. Local Arrays ..
306 DOUBLE PRECISION VT( 3 )
307* ..
308* .. External Subroutines ..
309 EXTERNAL dgemm, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm
310* ..
311* .. Executable Statements ..
312*
313* ==== If there are no shifts, then there is nothing to do. ====
314*
315 IF( nshfts.LT.2 )
316 \$ RETURN
317*
318* ==== If the active block is empty or 1-by-1, then there
319* . is nothing to do. ====
320*
321 IF( ktop.GE.kbot )
322 \$ RETURN
323*
324* ==== Shuffle shifts into pairs of real shifts and pairs
325* . of complex conjugate shifts assuming complex
327* . another. ====
328*
329 DO 10 i = 1, nshfts - 2, 2
330 IF( si( i ).NE.-si( i+1 ) ) THEN
331*
332 swap = sr( i )
333 sr( i ) = sr( i+1 )
334 sr( i+1 ) = sr( i+2 )
335 sr( i+2 ) = swap
336*
337 swap = si( i )
338 si( i ) = si( i+1 )
339 si( i+1 ) = si( i+2 )
340 si( i+2 ) = swap
341 END IF
342 10 CONTINUE
343*
344* ==== NSHFTS is supposed to be even, but if it is odd,
345* . then simply reduce it by one. The shuffle above
346* . ensures that the dropped shift is real and that
347* . the remaining shifts are paired. ====
348*
349 ns = nshfts - mod( nshfts, 2 )
350*
351* ==== Machine constants for deflation ====
352*
353 safmin = dlamch( 'SAFE MINIMUM' )
354 safmax = one / safmin
355 ulp = dlamch( 'PRECISION' )
356 smlnum = safmin*( dble( n ) / ulp )
357*
358* ==== Use accumulated reflections to update far-from-diagonal
359* . entries ? ====
360*
361 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
362*
363* ==== clear trash ====
364*
365 IF( ktop+2.LE.kbot )
366 \$ h( ktop+2, ktop ) = zero
367*
368* ==== NBMPS = number of 2-shift bulges in the chain ====
369*
370 nbmps = ns / 2
371*
372* ==== KDU = width of slab ====
373*
374 kdu = 4*nbmps
375*
376* ==== Create and chase chains of NBMPS bulges ====
377*
378 DO 180 incol = ktop - 2*nbmps + 1, kbot - 2, 2*nbmps
379*
380* JTOP = Index from which updates from the right start.
381*
382 IF( accum ) THEN
383 jtop = max( ktop, incol )
384 ELSE IF( wantt ) THEN
385 jtop = 1
386 ELSE
387 jtop = ktop
388 END IF
389*
390 ndcol = incol + kdu
391 IF( accum )
392 \$ CALL dlaset( 'ALL', kdu, kdu, zero, one, u, ldu )
393*
394* ==== Near-the-diagonal bulge chase. The following loop
395* . performs the near-the-diagonal part of a small bulge
396* . multi-shift QR sweep. Each 4*NBMPS column diagonal
397* . chunk extends from column INCOL to column NDCOL
398* . (including both column INCOL and column NDCOL). The
399* . following loop chases a 2*NBMPS+1 column long chain of
400* . NBMPS bulges 2*NBMPS columns to the right. (INCOL
401* . may be less than KTOP and and NDCOL may be greater than
402* . KBOT indicating phantom columns from which to chase
403* . bulges before they are actually introduced or to which
404* . to chase bulges beyond column KBOT.) ====
405*
406 DO 145 krcol = incol, min( incol+2*nbmps-1, kbot-2 )
407*
408* ==== Bulges number MTOP to MBOT are active double implicit
409* . shift bulges. There may or may not also be small
410* . 2-by-2 bulge, if there is room. The inactive bulges
411* . (if any) must wait until the active bulges have moved
412* . down the diagonal to make room. The phantom matrix
413* . paradigm described above helps keep track. ====
414*
415 mtop = max( 1, ( ktop-krcol ) / 2+1 )
416 mbot = min( nbmps, ( kbot-krcol-1 ) / 2 )
417 m22 = mbot + 1
418 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+2*( m22-1 ) ).EQ.
419 \$ ( kbot-2 )
420*
421* ==== Generate reflections to chase the chain right
422* . one column. (The minimum value of K is KTOP-1.) ====
423*
424 IF ( bmp22 ) THEN
425*
426* ==== Special case: 2-by-2 reflection at bottom treated
427* . separately ====
428*
429 k = krcol + 2*( m22-1 )
430 IF( k.EQ.ktop-1 ) THEN
431 CALL dlaqr1( 2, h( k+1, k+1 ), ldh, sr( 2*m22-1 ),
432 \$ si( 2*m22-1 ), sr( 2*m22 ), si( 2*m22 ),
433 \$ v( 1, m22 ) )
434 beta = v( 1, m22 )
435 CALL dlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
436 ELSE
437 beta = h( k+1, k )
438 v( 2, m22 ) = h( k+2, k )
439 CALL dlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
440 h( k+1, k ) = beta
441 h( k+2, k ) = zero
442 END IF
443
444*
445* ==== Perform update from right within
446* . computational window. ====
447*
448 t1 = v( 1, m22 )
449 t2 = t1*v( 2, m22 )
450 DO 30 j = jtop, min( kbot, k+3 )
451 refsum = h( j, k+1 ) + v( 2, m22 )*h( j, k+2 )
452 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
453 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
454 30 CONTINUE
455*
456* ==== Perform update from left within
457* . computational window. ====
458*
459 IF( accum ) THEN
460 jbot = min( ndcol, kbot )
461 ELSE IF( wantt ) THEN
462 jbot = n
463 ELSE
464 jbot = kbot
465 END IF
466 t1 = v( 1, m22 )
467 t2 = t1*v( 2, m22 )
468 DO 40 j = k+1, jbot
469 refsum = h( k+1, j ) + v( 2, m22 )*h( k+2, j )
470 h( k+1, j ) = h( k+1, j ) - refsum*t1
471 h( k+2, j ) = h( k+2, j ) - refsum*t2
472 40 CONTINUE
473*
474* ==== The following convergence test requires that
476* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
477* . criteria both be satisfied. The latter improves
478* . accuracy in some examples. Falling back on an
479* . alternate convergence criterion when TST1 or TST2
480* . is zero (as done here) is traditional but probably
481* . unnecessary. ====
482*
483 IF( k.GE.ktop ) THEN
484 IF( h( k+1, k ).NE.zero ) THEN
485 tst1 = abs( h( k, k ) ) + abs( h( k+1, k+1 ) )
486 IF( tst1.EQ.zero ) THEN
487 IF( k.GE.ktop+1 )
488 \$ tst1 = tst1 + abs( h( k, k-1 ) )
489 IF( k.GE.ktop+2 )
490 \$ tst1 = tst1 + abs( h( k, k-2 ) )
491 IF( k.GE.ktop+3 )
492 \$ tst1 = tst1 + abs( h( k, k-3 ) )
493 IF( k.LE.kbot-2 )
494 \$ tst1 = tst1 + abs( h( k+2, k+1 ) )
495 IF( k.LE.kbot-3 )
496 \$ tst1 = tst1 + abs( h( k+3, k+1 ) )
497 IF( k.LE.kbot-4 )
498 \$ tst1 = tst1 + abs( h( k+4, k+1 ) )
499 END IF
500 IF( abs( h( k+1, k ) )
501 \$ .LE.max( smlnum, ulp*tst1 ) ) THEN
502 h12 = max( abs( h( k+1, k ) ),
503 \$ abs( h( k, k+1 ) ) )
504 h21 = min( abs( h( k+1, k ) ),
505 \$ abs( h( k, k+1 ) ) )
506 h11 = max( abs( h( k+1, k+1 ) ),
507 \$ abs( h( k, k )-h( k+1, k+1 ) ) )
508 h22 = min( abs( h( k+1, k+1 ) ),
509 \$ abs( h( k, k )-h( k+1, k+1 ) ) )
510 scl = h11 + h12
511 tst2 = h22*( h11 / scl )
512*
513 IF( tst2.EQ.zero .OR. h21*( h12 / scl ).LE.
514 \$ max( smlnum, ulp*tst2 ) ) THEN
515 h( k+1, k ) = zero
516 END IF
517 END IF
518 END IF
519 END IF
520*
521* ==== Accumulate orthogonal transformations. ====
522*
523 IF( accum ) THEN
524 kms = k - incol
525 t1 = v( 1, m22 )
526 t2 = t1*v( 2, m22 )
527 DO 50 j = max( 1, ktop-incol ), kdu
528 refsum = u( j, kms+1 ) + v( 2, m22 )*u( j, kms+2 )
529 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
530 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
531 50 CONTINUE
532 ELSE IF( wantz ) THEN
533 t1 = v( 1, m22 )
534 t2 = t1*v( 2, m22 )
535 DO 60 j = iloz, ihiz
536 refsum = z( j, k+1 )+v( 2, m22 )*z( j, k+2 )
537 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
538 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
539 60 CONTINUE
540 END IF
541 END IF
542*
543* ==== Normal case: Chain of 3-by-3 reflections ====
544*
545 DO 80 m = mbot, mtop, -1
546 k = krcol + 2*( m-1 )
547 IF( k.EQ.ktop-1 ) THEN
548 CALL dlaqr1( 3, h( ktop, ktop ), ldh, sr( 2*m-1 ),
549 \$ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
550 \$ v( 1, m ) )
551 alpha = v( 1, m )
552 CALL dlarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
553 ELSE
554*
555* ==== Perform delayed transformation of row below
556* . Mth bulge. Exploit fact that first two elements
557* . of row are actually zero. ====
558*
559 t1 = v( 1, m )
560 t2 = t1*v( 2, m )
561 t3 = t1*v( 3, m )
562 refsum = v( 3, m )*h( k+3, k+2 )
563 h( k+3, k ) = -refsum*t1
564 h( k+3, k+1 ) = -refsum*t2
565 h( k+3, k+2 ) = h( k+3, k+2 ) - refsum*t3
566*
567* ==== Calculate reflection to move
568* . Mth bulge one step. ====
569*
570 beta = h( k+1, k )
571 v( 2, m ) = h( k+2, k )
572 v( 3, m ) = h( k+3, k )
573 CALL dlarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
574*
575* ==== A Bulge may collapse because of vigilant
576* . deflation or destructive underflow. In the
577* . underflow case, try the two-small-subdiagonals
578* . trick to try to reinflate the bulge. ====
579*
580 IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
581 \$ zero .OR. h( k+3, k+2 ).EQ.zero ) THEN
582*
583* ==== Typical case: not collapsed (yet). ====
584*
585 h( k+1, k ) = beta
586 h( k+2, k ) = zero
587 h( k+3, k ) = zero
588 ELSE
589*
590* ==== Atypical case: collapsed. Attempt to
591* . reintroduce ignoring H(K+1,K) and H(K+2,K).
592* . If the fill resulting from the new
593* . reflector is too large, then abandon it.
594* . Otherwise, use the new one. ====
595*
596 CALL dlaqr1( 3, h( k+1, k+1 ), ldh, sr( 2*m-1 ),
597 \$ si( 2*m-1 ), sr( 2*m ), si( 2*m ),
598 \$ vt )
599 alpha = vt( 1 )
600 CALL dlarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
601 t1 = vt( 1 )
602 t2 = t1*vt( 2 )
603 t3 = t1*vt( 3 )
604 refsum = h( k+1, k ) + vt( 2 )*h( k+2, k )
605*
606 IF( abs( h( k+2, k )-refsum*t2 )+
607 \$ abs( refsum*t3 ).GT.ulp*
608 \$ ( abs( h( k, k ) )+abs( h( k+1,
609 \$ k+1 ) )+abs( h( k+2, k+2 ) ) ) ) THEN
610*
611* ==== Starting a new bulge here would
612* . create non-negligible fill. Use
613* . the old one with trepidation. ====
614*
615 h( k+1, k ) = beta
616 h( k+2, k ) = zero
617 h( k+3, k ) = zero
618 ELSE
619*
620* ==== Starting a new bulge here would
621* . create only negligible fill.
622* . Replace the old reflector with
623* . the new one. ====
624*
625 h( k+1, k ) = h( k+1, k ) - refsum*t1
626 h( k+2, k ) = zero
627 h( k+3, k ) = zero
628 v( 1, m ) = vt( 1 )
629 v( 2, m ) = vt( 2 )
630 v( 3, m ) = vt( 3 )
631 END IF
632 END IF
633 END IF
634*
635* ==== Apply reflection from the right and
636* . the first column of update from the left.
637* . These updates are required for the vigilant
638* . deflation check. We still delay most of the
639* . updates from the left for efficiency. ====
640*
641 t1 = v( 1, m )
642 t2 = t1*v( 2, m )
643 t3 = t1*v( 3, m )
644 DO 70 j = jtop, min( kbot, k+3 )
645 refsum = h( j, k+1 ) + v( 2, m )*h( j, k+2 )
646 \$ + v( 3, m )*h( j, k+3 )
647 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
648 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
649 h( j, k+3 ) = h( j, k+3 ) - refsum*t3
650 70 CONTINUE
651*
652* ==== Perform update from left for subsequent
653* . column. ====
654*
655 refsum = h( k+1, k+1 ) + v( 2, m )*h( k+2, k+1 )
656 \$ + v( 3, m )*h( k+3, k+1 )
657 h( k+1, k+1 ) = h( k+1, k+1 ) - refsum*t1
658 h( k+2, k+1 ) = h( k+2, k+1 ) - refsum*t2
659 h( k+3, k+1 ) = h( k+3, k+1 ) - refsum*t3
660*
661* ==== The following convergence test requires that
663* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
664* . criteria both be satisfied. The latter improves
665* . accuracy in some examples. Falling back on an
666* . alternate convergence criterion when TST1 or TST2
667* . is zero (as done here) is traditional but probably
668* . unnecessary. ====
669*
670 IF( k.LT.ktop)
671 \$ cycle
672 IF( h( k+1, k ).NE.zero ) THEN
673 tst1 = abs( h( k, k ) ) + abs( h( k+1, k+1 ) )
674 IF( tst1.EQ.zero ) THEN
675 IF( k.GE.ktop+1 )
676 \$ tst1 = tst1 + abs( h( k, k-1 ) )
677 IF( k.GE.ktop+2 )
678 \$ tst1 = tst1 + abs( h( k, k-2 ) )
679 IF( k.GE.ktop+3 )
680 \$ tst1 = tst1 + abs( h( k, k-3 ) )
681 IF( k.LE.kbot-2 )
682 \$ tst1 = tst1 + abs( h( k+2, k+1 ) )
683 IF( k.LE.kbot-3 )
684 \$ tst1 = tst1 + abs( h( k+3, k+1 ) )
685 IF( k.LE.kbot-4 )
686 \$ tst1 = tst1 + abs( h( k+4, k+1 ) )
687 END IF
688 IF( abs( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
689 \$ THEN
690 h12 = max( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
691 h21 = min( abs( h( k+1, k ) ), abs( h( k, k+1 ) ) )
692 h11 = max( abs( h( k+1, k+1 ) ),
693 \$ abs( h( k, k )-h( k+1, k+1 ) ) )
694 h22 = min( abs( h( k+1, k+1 ) ),
695 \$ abs( h( k, k )-h( k+1, k+1 ) ) )
696 scl = h11 + h12
697 tst2 = h22*( h11 / scl )
698*
699 IF( tst2.EQ.zero .OR. h21*( h12 / scl ).LE.
700 \$ max( smlnum, ulp*tst2 ) ) THEN
701 h( k+1, k ) = zero
702 END IF
703 END IF
704 END IF
705 80 CONTINUE
706*
707* ==== Multiply H by reflections from the left ====
708*
709 IF( accum ) THEN
710 jbot = min( ndcol, kbot )
711 ELSE IF( wantt ) THEN
712 jbot = n
713 ELSE
714 jbot = kbot
715 END IF
716*
717 DO 100 m = mbot, mtop, -1
718 k = krcol + 2*( m-1 )
719 t1 = v( 1, m )
720 t2 = t1*v( 2, m )
721 t3 = t1*v( 3, m )
722 DO 90 j = max( ktop, krcol + 2*m ), jbot
723 refsum = h( k+1, j ) + v( 2, m )*h( k+2, j )
724 \$ + v( 3, m )*h( k+3, j )
725 h( k+1, j ) = h( k+1, j ) - refsum*t1
726 h( k+2, j ) = h( k+2, j ) - refsum*t2
727 h( k+3, j ) = h( k+3, j ) - refsum*t3
728 90 CONTINUE
729 100 CONTINUE
730*
731* ==== Accumulate orthogonal transformations. ====
732*
733 IF( accum ) THEN
734*
735* ==== Accumulate U. (If needed, update Z later
736* . with an efficient matrix-matrix
737* . multiply.) ====
738*
739 DO 120 m = mbot, mtop, -1
740 k = krcol + 2*( m-1 )
741 kms = k - incol
742 i2 = max( 1, ktop-incol )
743 i2 = max( i2, kms-(krcol-incol)+1 )
744 i4 = min( kdu, krcol + 2*( mbot-1 ) - incol + 5 )
745 t1 = v( 1, m )
746 t2 = t1*v( 2, m )
747 t3 = t1*v( 3, m )
748 DO 110 j = i2, i4
749 refsum = u( j, kms+1 ) + v( 2, m )*u( j, kms+2 )
750 \$ + v( 3, m )*u( j, kms+3 )
751 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
752 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
753 u( j, kms+3 ) = u( j, kms+3 ) - refsum*t3
754 110 CONTINUE
755 120 CONTINUE
756 ELSE IF( wantz ) THEN
757*
758* ==== U is not accumulated, so update Z
759* . now by multiplying by reflections
760* . from the right. ====
761*
762 DO 140 m = mbot, mtop, -1
763 k = krcol + 2*( m-1 )
764 t1 = v( 1, m )
765 t2 = t1*v( 2, m )
766 t3 = t1*v( 3, m )
767 DO 130 j = iloz, ihiz
768 refsum = z( j, k+1 ) + v( 2, m )*z( j, k+2 )
769 \$ + v( 3, m )*z( j, k+3 )
770 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
771 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
772 z( j, k+3 ) = z( j, k+3 ) - refsum*t3
773 130 CONTINUE
774 140 CONTINUE
775 END IF
776*
777* ==== End of near-the-diagonal bulge chase. ====
778*
779 145 CONTINUE
780*
781* ==== Use U (if accumulated) to update far-from-diagonal
782* . entries in H. If required, use U to update Z as
783* . well. ====
784*
785 IF( accum ) THEN
786 IF( wantt ) THEN
787 jtop = 1
788 jbot = n
789 ELSE
790 jtop = ktop
791 jbot = kbot
792 END IF
793 k1 = max( 1, ktop-incol )
794 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
795*
796* ==== Horizontal Multiply ====
797*
798 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
799 jlen = min( nh, jbot-jcol+1 )
800 CALL dgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
801 \$ ldu, h( incol+k1, jcol ), ldh, zero, wh,
802 \$ ldwh )
803 CALL dlacpy( 'ALL', nu, jlen, wh, ldwh,
804 \$ h( incol+k1, jcol ), ldh )
805 150 CONTINUE
806*
807* ==== Vertical multiply ====
808*
809 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
810 jlen = min( nv, max( ktop, incol )-jrow )
811 CALL dgemm( 'N', 'N', jlen, nu, nu, one,
812 \$ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
813 \$ ldu, zero, wv, ldwv )
814 CALL dlacpy( 'ALL', jlen, nu, wv, ldwv,
815 \$ h( jrow, incol+k1 ), ldh )
816 160 CONTINUE
817*
818* ==== Z multiply (also vertical) ====
819*
820 IF( wantz ) THEN
821 DO 170 jrow = iloz, ihiz, nv
822 jlen = min( nv, ihiz-jrow+1 )
823 CALL dgemm( 'N', 'N', jlen, nu, nu, one,
824 \$ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
825 \$ ldu, zero, wv, ldwv )
826 CALL dlacpy( 'ALL', jlen, nu, wv, ldwv,
827 \$ z( jrow, incol+k1 ), ldz )
828 170 CONTINUE
829 END IF
830 END IF
831 180 CONTINUE
832*
833* ==== End of DLAQR5 ====
834*
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
subroutine dlaqr1(n, h, ldh, sr1, si1, sr2, si2, v)
DLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
Definition dlaqr1.f:121
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: