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

◆ zlaqr5()

subroutine zlaqr5 ( logical  wantt,
logical  wantz,
integer  kacc22,
integer  n,
integer  ktop,
integer  kbot,
integer  nshfts,
complex*16, dimension( * )  s,
complex*16, dimension( ldh, * )  h,
integer  ldh,
integer  iloz,
integer  ihiz,
complex*16, dimension( ldz, * )  z,
integer  ldz,
complex*16, dimension( ldv, * )  v,
integer  ldv,
complex*16, dimension( ldu, * )  u,
integer  ldu,
integer  nv,
complex*16, dimension( ldwv, * )  wv,
integer  ldwv,
integer  nh,
complex*16, dimension( ldwh, * )  wh,
integer  ldwh 
)

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

Download ZLAQR5 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
    ZLAQR5, called by ZLAQR0, performs a
    single small-bulge multi-shift QR sweep.
Parameters
[in]WANTT
          WANTT is LOGICAL
             WANTT = .true. if the triangular Schur factor
             is being computed.  WANTT is set to .false. otherwise.
[in]WANTZ
          WANTZ is LOGICAL
             WANTZ = .true. if the unitary 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: ZLAQR5 does not accumulate reflections and does not
             use matrix-matrix multiply to update far-from-diagonal
             matrix entries.
        = 1: ZLAQR5 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]S
          S is COMPLEX*16 array, dimension (NSHFTS)
             S contains the shifts of origin that define the multi-
             shift QR sweep.  On output S may be reordered.
[in,out]H
          H is COMPLEX*16 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 COMPLEX*16 array, dimension (LDZ,IHIZ)
             If WANTZ = .TRUE., then the QR Sweep unitary
             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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDWH,NH)
[in]LDWH
          LDWH is INTEGER
             Leading dimension of WH just as declared in the
             calling procedure.  LDWH >= 2*NSHFTS.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 254 of file zlaqr5.f.

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