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

◆ claqr5()

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

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

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

Purpose:
    CLAQR5 called by CLAQR0 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: CLAQR5 does not accumulate reflections and does not
             use matrix-matrix multiply to update far-from-diagonal
             matrix entries.
        = 1: CLAQR5 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 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 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 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 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 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 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 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 claqr5.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 H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
270 $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
271* ..
272*
273* ================================================================
274* .. Parameters ..
275 COMPLEX ZERO, ONE
276 parameter( zero = ( 0.0e0, 0.0e0 ),
277 $ one = ( 1.0e0, 0.0e0 ) )
278 REAL RZERO, RONE
279 parameter( rzero = 0.0e0, rone = 1.0e0 )
280* ..
281* .. Local Scalars ..
282 COMPLEX ALPHA, BETA, CDUM, REFSUM, T1, T2, T3
283 REAL 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 REAL SLAMCH
293 EXTERNAL slamch
294* ..
295* .. Intrinsic Functions ..
296*
297 INTRINSIC abs, aimag, conjg, max, min, mod, real
298* ..
299* .. Local Arrays ..
300 COMPLEX VT( 3 )
301* ..
302* .. External Subroutines ..
303 EXTERNAL cgemm, clacpy, claqr1, clarfg, claset, ctrmm
304* ..
305* .. Statement Functions ..
306 REAL CABS1
307* ..
308* .. Statement Function definitions ..
309 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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 = slamch( 'SAFE MINIMUM' )
332 safmax = rone / safmin
333 ulp = slamch( 'PRECISION' )
334 smlnum = safmin*( real( 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 claset( '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 claqr1( 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 clarfg( 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 clarfg( 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*conjg( 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 = conjg( v( 1, m22 ) )
444 t2 = t1*v( 2, m22 )
445 DO 40 j = k+1, jbot
446 refsum = h( k+1, j ) +
447 $ conjg( 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*conjg( 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*conjg( 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 claqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
525 $ s( 2*m ), v( 1, m ) )
526 alpha = v( 1, m )
527 CALL clarfg( 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*conjg( v( 2, m ) )
536 t3 = t1*conjg( 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 clarfg( 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 claqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
572 $ s( 2*m ), vt )
573 alpha = vt( 1 )
574 CALL clarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
575 t1 = conjg( vt( 1 ) )
576 t2 = t1*vt( 2 )
577 t3 = t1*vt( 3 )
578 refsum = h( k+1, k )+conjg( 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*conjg( v( 2, m ) )
617 t3 = t1*conjg( 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 = conjg( v( 1, m ) )
630 t2 = t1*v( 2, m )
631 t3 = t1*v( 3, m )
632 refsum = h( k+1, k+1 ) + conjg( v( 2, m ) )*h( k+2, k+1 )
633 $ + conjg( v( 3, m ) )*h( k+3, k+1 )
634 h( k+1, k+1 ) = h( k+1, k+1 ) - refsum*t1
635 h( k+2, k+1 ) = h( k+2, k+1 ) - refsum*t2
636 h( k+3, k+1 ) = h( k+3, k+1 ) - refsum*t3
637*
638* ==== The following convergence test requires that
639* . the tradition small-compared-to-nearby-diagonals
640* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
641* . criteria both be satisfied. The latter improves
642* . accuracy in some examples. Falling back on an
643* . alternate convergence criterion when TST1 or TST2
644* . is zero (as done here) is traditional but probably
645* . unnecessary. ====
646*
647 IF( k.LT.ktop)
648 $ cycle
649 IF( h( k+1, k ).NE.zero ) THEN
650 tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
651 IF( tst1.EQ.rzero ) THEN
652 IF( k.GE.ktop+1 )
653 $ tst1 = tst1 + cabs1( h( k, k-1 ) )
654 IF( k.GE.ktop+2 )
655 $ tst1 = tst1 + cabs1( h( k, k-2 ) )
656 IF( k.GE.ktop+3 )
657 $ tst1 = tst1 + cabs1( h( k, k-3 ) )
658 IF( k.LE.kbot-2 )
659 $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
660 IF( k.LE.kbot-3 )
661 $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
662 IF( k.LE.kbot-4 )
663 $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
664 END IF
665 IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
666 $ THEN
667 h12 = max( cabs1( h( k+1, k ) ),
668 $ cabs1( h( k, k+1 ) ) )
669 h21 = min( cabs1( h( k+1, k ) ),
670 $ cabs1( h( k, k+1 ) ) )
671 h11 = max( cabs1( h( k+1, k+1 ) ),
672 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
673 h22 = min( cabs1( h( k+1, k+1 ) ),
674 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
675 scl = h11 + h12
676 tst2 = h22*( h11 / scl )
677*
678 IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
679 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
680 END IF
681 END IF
682 80 CONTINUE
683*
684* ==== Multiply H by reflections from the left ====
685*
686 IF( accum ) THEN
687 jbot = min( ndcol, kbot )
688 ELSE IF( wantt ) THEN
689 jbot = n
690 ELSE
691 jbot = kbot
692 END IF
693*
694 DO 100 m = mbot, mtop, -1
695 k = krcol + 2*( m-1 )
696 t1 = conjg( v( 1, m ) )
697 t2 = t1*v( 2, m )
698 t3 = t1*v( 3, m )
699 DO 90 j = max( ktop, krcol + 2*m ), jbot
700 refsum = h( k+1, j ) + conjg( v( 2, m ) )*
701 $ h( k+2, j ) + conjg( v( 3, m ) )*h( k+3, j )
702 h( k+1, j ) = h( k+1, j ) - refsum*t1
703 h( k+2, j ) = h( k+2, j ) - refsum*t2
704 h( k+3, j ) = h( k+3, j ) - refsum*t3
705 90 CONTINUE
706 100 CONTINUE
707*
708* ==== Accumulate orthogonal transformations. ====
709*
710 IF( accum ) THEN
711*
712* ==== Accumulate U. (If needed, update Z later
713* . with an efficient matrix-matrix
714* . multiply.) ====
715*
716 DO 120 m = mbot, mtop, -1
717 k = krcol + 2*( m-1 )
718 kms = k - incol
719 i2 = max( 1, ktop-incol )
720 i2 = max( i2, kms-(krcol-incol)+1 )
721 i4 = min( kdu, krcol + 2*( mbot-1 ) - incol + 5 )
722 t1 = v( 1, m )
723 t2 = t1*conjg( v( 2, m ) )
724 t3 = t1*conjg( v( 3, m ) )
725 DO 110 j = i2, i4
726 refsum = u( j, kms+1 ) + v( 2, m )*u( j, kms+2 )
727 $ + v( 3, m )*u( j, kms+3 )
728 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
729 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
730 u( j, kms+3 ) = u( j, kms+3 ) - refsum*t3
731 110 CONTINUE
732 120 CONTINUE
733 ELSE IF( wantz ) THEN
734*
735* ==== U is not accumulated, so update Z
736* . now by multiplying by reflections
737* . from the right. ====
738*
739 DO 140 m = mbot, mtop, -1
740 k = krcol + 2*( m-1 )
741 t1 = v( 1, m )
742 t2 = t1*conjg( v( 2, m ) )
743 t3 = t1*conjg( v( 3, m ) )
744 DO 130 j = iloz, ihiz
745 refsum = z( j, k+1 ) + v( 2, m )*z( j, k+2 )
746 $ + v( 3, m )*z( j, k+3 )
747 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
748 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
749 z( j, k+3 ) = z( j, k+3 ) - refsum*t3
750 130 CONTINUE
751 140 CONTINUE
752 END IF
753*
754* ==== End of near-the-diagonal bulge chase. ====
755*
756 145 CONTINUE
757*
758* ==== Use U (if accumulated) to update far-from-diagonal
759* . entries in H. If required, use U to update Z as
760* . well. ====
761*
762 IF( accum ) THEN
763 IF( wantt ) THEN
764 jtop = 1
765 jbot = n
766 ELSE
767 jtop = ktop
768 jbot = kbot
769 END IF
770 k1 = max( 1, ktop-incol )
771 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
772*
773* ==== Horizontal Multiply ====
774*
775 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
776 jlen = min( nh, jbot-jcol+1 )
777 CALL cgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
778 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
779 $ ldwh )
780 CALL clacpy( 'ALL', nu, jlen, wh, ldwh,
781 $ h( incol+k1, jcol ), ldh )
782 150 CONTINUE
783*
784* ==== Vertical multiply ====
785*
786 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
787 jlen = min( nv, max( ktop, incol )-jrow )
788 CALL cgemm( 'N', 'N', jlen, nu, nu, one,
789 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
790 $ ldu, zero, wv, ldwv )
791 CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
792 $ h( jrow, incol+k1 ), ldh )
793 160 CONTINUE
794*
795* ==== Z multiply (also vertical) ====
796*
797 IF( wantz ) THEN
798 DO 170 jrow = iloz, ihiz, nv
799 jlen = min( nv, ihiz-jrow+1 )
800 CALL cgemm( 'N', 'N', jlen, nu, nu, one,
801 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
802 $ ldu, zero, wv, ldwv )
803 CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
804 $ z( jrow, incol+k1 ), ldz )
805 170 CONTINUE
806 END IF
807 END IF
808 180 CONTINUE
809*
810* ==== End of CLAQR5 ====
811*
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
claqr1
subroutine claqr1(n, h, ldh, s1, s2, v)
CLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
Definition claqr1.f:107
clarfg
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
ctrmm
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: