LAPACK 3.11.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 $ slabad
305* ..
306* .. Statement Functions ..
307 REAL CABS1
308* ..
309* .. Statement Function definitions ..
310 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
311* ..
312* .. Executable Statements ..
313*
314* ==== If there are no shifts, then there is nothing to do. ====
315*
316 IF( nshfts.LT.2 )
317 $ RETURN
318*
319* ==== If the active block is empty or 1-by-1, then there
320* . is nothing to do. ====
321*
322 IF( ktop.GE.kbot )
323 $ RETURN
324*
325* ==== NSHFTS is supposed to be even, but if it is odd,
326* . then simply reduce it by one. ====
327*
328 ns = nshfts - mod( nshfts, 2 )
329*
330* ==== Machine constants for deflation ====
331*
332 safmin = slamch( 'SAFE MINIMUM' )
333 safmax = rone / safmin
334 CALL slabad( safmin, safmax )
335 ulp = slamch( 'PRECISION' )
336 smlnum = safmin*( real( n ) / ulp )
337*
338* ==== Use accumulated reflections to update far-from-diagonal
339* . entries ? ====
340*
341 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
342*
343* ==== clear trash ====
344*
345 IF( ktop+2.LE.kbot )
346 $ h( ktop+2, ktop ) = zero
347*
348* ==== NBMPS = number of 2-shift bulges in the chain ====
349*
350 nbmps = ns / 2
351*
352* ==== KDU = width of slab ====
353*
354 kdu = 4*nbmps
355*
356* ==== Create and chase chains of NBMPS bulges ====
357*
358 DO 180 incol = ktop - 2*nbmps + 1, kbot - 2, 2*nbmps
359*
360* JTOP = Index from which updates from the right start.
361*
362 IF( accum ) THEN
363 jtop = max( ktop, incol )
364 ELSE IF( wantt ) THEN
365 jtop = 1
366 ELSE
367 jtop = ktop
368 END IF
369*
370 ndcol = incol + kdu
371 IF( accum )
372 $ CALL claset( 'ALL', kdu, kdu, zero, one, u, ldu )
373*
374* ==== Near-the-diagonal bulge chase. The following loop
375* . performs the near-the-diagonal part of a small bulge
376* . multi-shift QR sweep. Each 4*NBMPS column diagonal
377* . chunk extends from column INCOL to column NDCOL
378* . (including both column INCOL and column NDCOL). The
379* . following loop chases a 2*NBMPS+1 column long chain of
380* . NBMPS bulges 2*NBMPS columns to the right. (INCOL
381* . may be less than KTOP and and NDCOL may be greater than
382* . KBOT indicating phantom columns from which to chase
383* . bulges before they are actually introduced or to which
384* . to chase bulges beyond column KBOT.) ====
385*
386 DO 145 krcol = incol, min( incol+2*nbmps-1, kbot-2 )
387*
388* ==== Bulges number MTOP to MBOT are active double implicit
389* . shift bulges. There may or may not also be small
390* . 2-by-2 bulge, if there is room. The inactive bulges
391* . (if any) must wait until the active bulges have moved
392* . down the diagonal to make room. The phantom matrix
393* . paradigm described above helps keep track. ====
394*
395 mtop = max( 1, ( ktop-krcol ) / 2+1 )
396 mbot = min( nbmps, ( kbot-krcol-1 ) / 2 )
397 m22 = mbot + 1
398 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+2*( m22-1 ) ).EQ.
399 $ ( kbot-2 )
400*
401* ==== Generate reflections to chase the chain right
402* . one column. (The minimum value of K is KTOP-1.) ====
403*
404 IF ( bmp22 ) THEN
405*
406* ==== Special case: 2-by-2 reflection at bottom treated
407* . separately ====
408*
409 k = krcol + 2*( m22-1 )
410 IF( k.EQ.ktop-1 ) THEN
411 CALL claqr1( 2, h( k+1, k+1 ), ldh, s( 2*m22-1 ),
412 $ s( 2*m22 ), v( 1, m22 ) )
413 beta = v( 1, m22 )
414 CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
415 ELSE
416 beta = h( k+1, k )
417 v( 2, m22 ) = h( k+2, k )
418 CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
419 h( k+1, k ) = beta
420 h( k+2, k ) = zero
421 END IF
422
423*
424* ==== Perform update from right within
425* . computational window. ====
426*
427 t1 = v( 1, m22 )
428 t2 = t1*conjg( v( 2, m22 ) )
429 DO 30 j = jtop, min( kbot, k+3 )
430 refsum = h( j, k+1 ) + v( 2, m22 )*h( j, k+2 )
431 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
432 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
433 30 CONTINUE
434*
435* ==== Perform update from left within
436* . computational window. ====
437*
438 IF( accum ) THEN
439 jbot = min( ndcol, kbot )
440 ELSE IF( wantt ) THEN
441 jbot = n
442 ELSE
443 jbot = kbot
444 END IF
445 t1 = conjg( v( 1, m22 ) )
446 t2 = t1*v( 2, m22 )
447 DO 40 j = k+1, jbot
448 refsum = h( k+1, j ) +
449 $ conjg( v( 2, m22 ) )*h( k+2, j )
450 h( k+1, j ) = h( k+1, j ) - refsum*t1
451 h( k+2, j ) = h( k+2, j ) - refsum*t2
452 40 CONTINUE
453*
454* ==== The following convergence test requires that
455* . the tradition small-compared-to-nearby-diagonals
456* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
457* . criteria both be satisfied. The latter improves
458* . accuracy in some examples. Falling back on an
459* . alternate convergence criterion when TST1 or TST2
460* . is zero (as done here) is traditional but probably
461* . unnecessary. ====
462*
463 IF( k.GE.ktop) THEN
464 IF( h( k+1, k ).NE.zero ) THEN
465 tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
466 IF( tst1.EQ.rzero ) THEN
467 IF( k.GE.ktop+1 )
468 $ tst1 = tst1 + cabs1( h( k, k-1 ) )
469 IF( k.GE.ktop+2 )
470 $ tst1 = tst1 + cabs1( h( k, k-2 ) )
471 IF( k.GE.ktop+3 )
472 $ tst1 = tst1 + cabs1( h( k, k-3 ) )
473 IF( k.LE.kbot-2 )
474 $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
475 IF( k.LE.kbot-3 )
476 $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
477 IF( k.LE.kbot-4 )
478 $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
479 END IF
480 IF( cabs1( h( k+1, k ) )
481 $ .LE.max( smlnum, ulp*tst1 ) ) THEN
482 h12 = max( cabs1( h( k+1, k ) ),
483 $ cabs1( h( k, k+1 ) ) )
484 h21 = min( cabs1( h( k+1, k ) ),
485 $ cabs1( h( k, k+1 ) ) )
486 h11 = max( cabs1( h( k+1, k+1 ) ),
487 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
488 h22 = min( cabs1( h( k+1, k+1 ) ),
489 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
490 scl = h11 + h12
491 tst2 = h22*( h11 / scl )
492*
493 IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
494 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
495 END IF
496 END IF
497 END IF
498*
499* ==== Accumulate orthogonal transformations. ====
500*
501 IF( accum ) THEN
502 kms = k - incol
503 DO 50 j = max( 1, ktop-incol ), kdu
504 refsum = v( 1, m22 )*( u( j, kms+1 )+
505 $ v( 2, m22 )*u( j, kms+2 ) )
506 u( j, kms+1 ) = u( j, kms+1 ) - refsum
507 u( j, kms+2 ) = u( j, kms+2 ) -
508 $ refsum*conjg( v( 2, m22 ) )
509 50 CONTINUE
510 ELSE IF( wantz ) THEN
511 DO 60 j = iloz, ihiz
512 refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
513 $ z( j, k+2 ) )
514 z( j, k+1 ) = z( j, k+1 ) - refsum
515 z( j, k+2 ) = z( j, k+2 ) -
516 $ refsum*conjg( v( 2, m22 ) )
517 60 CONTINUE
518 END IF
519 END IF
520*
521* ==== Normal case: Chain of 3-by-3 reflections ====
522*
523 DO 80 m = mbot, mtop, -1
524 k = krcol + 2*( m-1 )
525 IF( k.EQ.ktop-1 ) THEN
526 CALL claqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
527 $ s( 2*m ), v( 1, m ) )
528 alpha = v( 1, m )
529 CALL clarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
530 ELSE
531*
532* ==== Perform delayed transformation of row below
533* . Mth bulge. Exploit fact that first two elements
534* . of row are actually zero. ====
535*
536 refsum = v( 1, m )*v( 3, m )*h( k+3, k+2 )
537 h( k+3, k ) = -refsum
538 h( k+3, k+1 ) = -refsum*conjg( v( 2, m ) )
539 h( k+3, k+2 ) = h( k+3, k+2 ) -
540 $ refsum*conjg( v( 3, m ) )
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 refsum = conjg( vt( 1 ) )*
576 $ ( h( k+1, k )+conjg( vt( 2 ) )*
577 $ h( k+2, k ) )
578*
579 IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
580 $ cabs1( refsum*vt( 3 ) ).GT.ulp*
581 $ ( cabs1( h( k, k ) )+cabs1( h( k+1,
582 $ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) ) THEN
583*
584* ==== Starting a new bulge here would
585* . create non-negligible fill. Use
586* . the old one with trepidation. ====
587*
588 h( k+1, k ) = beta
589 h( k+2, k ) = zero
590 h( k+3, k ) = zero
591 ELSE
592*
593* ==== Starting a new bulge here would
594* . create only negligible fill.
595* . Replace the old reflector with
596* . the new one. ====
597*
598 h( k+1, k ) = h( k+1, k ) - refsum
599 h( k+2, k ) = zero
600 h( k+3, k ) = zero
601 v( 1, m ) = vt( 1 )
602 v( 2, m ) = vt( 2 )
603 v( 3, m ) = vt( 3 )
604 END IF
605 END IF
606 END IF
607*
608* ==== Apply reflection from the right and
609* . the first column of update from the left.
610* . These updates are required for the vigilant
611* . deflation check. We still delay most of the
612* . updates from the left for efficiency. ====
613*
614 t1 = v( 1, m )
615 t2 = t1*conjg( v( 2, m ) )
616 t3 = t1*conjg( v( 3, m ) )
617 DO 70 j = jtop, min( kbot, k+3 )
618 refsum = h( j, k+1 ) + v( 2, m )*h( j, k+2 )
619 $ + v( 3, m )*h( j, k+3 )
620 h( j, k+1 ) = h( j, k+1 ) - refsum*t1
621 h( j, k+2 ) = h( j, k+2 ) - refsum*t2
622 h( j, k+3 ) = h( j, k+3 ) - refsum*t3
623 70 CONTINUE
624*
625* ==== Perform update from left for subsequent
626* . column. ====
627*
628 t1 = conjg( v( 1, m ) )
629 t2 = t1*v( 2, m )
630 t3 = t1*v( 3, m )
631 refsum = h( k+1, k+1 ) + conjg( v( 2, m ) )*h( k+2, k+1 )
632 $ + conjg( v( 3, m ) )*h( k+3, k+1 )
633 h( k+1, k+1 ) = h( k+1, k+1 ) - refsum*t1
634 h( k+2, k+1 ) = h( k+2, k+1 ) - refsum*t2
635 h( k+3, k+1 ) = h( k+3, k+1 ) - refsum*t3
636*
637* ==== The following convergence test requires that
638* . the tradition small-compared-to-nearby-diagonals
639* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
640* . criteria both be satisfied. The latter improves
641* . accuracy in some examples. Falling back on an
642* . alternate convergence criterion when TST1 or TST2
643* . is zero (as done here) is traditional but probably
644* . unnecessary. ====
645*
646 IF( k.LT.ktop)
647 $ cycle
648 IF( h( k+1, k ).NE.zero ) THEN
649 tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
650 IF( tst1.EQ.rzero ) THEN
651 IF( k.GE.ktop+1 )
652 $ tst1 = tst1 + cabs1( h( k, k-1 ) )
653 IF( k.GE.ktop+2 )
654 $ tst1 = tst1 + cabs1( h( k, k-2 ) )
655 IF( k.GE.ktop+3 )
656 $ tst1 = tst1 + cabs1( h( k, k-3 ) )
657 IF( k.LE.kbot-2 )
658 $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
659 IF( k.LE.kbot-3 )
660 $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
661 IF( k.LE.kbot-4 )
662 $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
663 END IF
664 IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
665 $ THEN
666 h12 = max( cabs1( h( k+1, k ) ),
667 $ cabs1( h( k, k+1 ) ) )
668 h21 = min( cabs1( h( k+1, k ) ),
669 $ cabs1( h( k, k+1 ) ) )
670 h11 = max( cabs1( h( k+1, k+1 ) ),
671 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
672 h22 = min( cabs1( h( k+1, k+1 ) ),
673 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
674 scl = h11 + h12
675 tst2 = h22*( h11 / scl )
676*
677 IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
678 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
679 END IF
680 END IF
681 80 CONTINUE
682*
683* ==== Multiply H by reflections from the left ====
684*
685 IF( accum ) THEN
686 jbot = min( ndcol, kbot )
687 ELSE IF( wantt ) THEN
688 jbot = n
689 ELSE
690 jbot = kbot
691 END IF
692*
693 DO 100 m = mbot, mtop, -1
694 k = krcol + 2*( m-1 )
695 t1 = conjg( v( 1, m ) )
696 t2 = t1*v( 2, m )
697 t3 = t1*v( 3, m )
698 DO 90 j = max( ktop, krcol + 2*m ), jbot
699 refsum = h( k+1, j ) + conjg( v( 2, m ) )*
700 $ h( k+2, j ) + conjg( v( 3, m ) )*h( k+3, j )
701 h( k+1, j ) = h( k+1, j ) - refsum*t1
702 h( k+2, j ) = h( k+2, j ) - refsum*t2
703 h( k+3, j ) = h( k+3, j ) - refsum*t3
704 90 CONTINUE
705 100 CONTINUE
706*
707* ==== Accumulate orthogonal transformations. ====
708*
709 IF( accum ) THEN
710*
711* ==== Accumulate U. (If needed, update Z later
712* . with an efficient matrix-matrix
713* . multiply.) ====
714*
715 DO 120 m = mbot, mtop, -1
716 k = krcol + 2*( m-1 )
717 kms = k - incol
718 i2 = max( 1, ktop-incol )
719 i2 = max( i2, kms-(krcol-incol)+1 )
720 i4 = min( kdu, krcol + 2*( mbot-1 ) - incol + 5 )
721 t1 = v( 1, m )
722 t2 = t1*conjg( v( 2, m ) )
723 t3 = t1*conjg( v( 3, m ) )
724 DO 110 j = i2, i4
725 refsum = u( j, kms+1 ) + v( 2, m )*u( j, kms+2 )
726 $ + v( 3, m )*u( j, kms+3 )
727 u( j, kms+1 ) = u( j, kms+1 ) - refsum*t1
728 u( j, kms+2 ) = u( j, kms+2 ) - refsum*t2
729 u( j, kms+3 ) = u( j, kms+3 ) - refsum*t3
730 110 CONTINUE
731 120 CONTINUE
732 ELSE IF( wantz ) THEN
733*
734* ==== U is not accumulated, so update Z
735* . now by multiplying by reflections
736* . from the right. ====
737*
738 DO 140 m = mbot, mtop, -1
739 k = krcol + 2*( m-1 )
740 t1 = v( 1, m )
741 t2 = t1*conjg( v( 2, m ) )
742 t3 = t1*conjg( v( 3, m ) )
743 DO 130 j = iloz, ihiz
744 refsum = z( j, k+1 ) + v( 2, m )*z( j, k+2 )
745 $ + v( 3, m )*z( j, k+3 )
746 z( j, k+1 ) = z( j, k+1 ) - refsum*t1
747 z( j, k+2 ) = z( j, k+2 ) - refsum*t2
748 z( j, k+3 ) = z( j, k+3 ) - refsum*t3
749 130 CONTINUE
750 140 CONTINUE
751 END IF
752*
753* ==== End of near-the-diagonal bulge chase. ====
754*
755 145 CONTINUE
756*
757* ==== Use U (if accumulated) to update far-from-diagonal
758* . entries in H. If required, use U to update Z as
759* . well. ====
760*
761 IF( accum ) THEN
762 IF( wantt ) THEN
763 jtop = 1
764 jbot = n
765 ELSE
766 jtop = ktop
767 jbot = kbot
768 END IF
769 k1 = max( 1, ktop-incol )
770 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
771*
772* ==== Horizontal Multiply ====
773*
774 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
775 jlen = min( nh, jbot-jcol+1 )
776 CALL cgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
777 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
778 $ ldwh )
779 CALL clacpy( 'ALL', nu, jlen, wh, ldwh,
780 $ h( incol+k1, jcol ), ldh )
781 150 CONTINUE
782*
783* ==== Vertical multiply ====
784*
785 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
786 jlen = min( nv, max( ktop, incol )-jrow )
787 CALL cgemm( 'N', 'N', jlen, nu, nu, one,
788 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
789 $ ldu, zero, wv, ldwv )
790 CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
791 $ h( jrow, incol+k1 ), ldh )
792 160 CONTINUE
793*
794* ==== Z multiply (also vertical) ====
795*
796 IF( wantz ) THEN
797 DO 170 jrow = iloz, ihiz, nv
798 jlen = min( nv, ihiz-jrow+1 )
799 CALL cgemm( 'N', 'N', jlen, nu, nu, one,
800 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
801 $ ldu, zero, wv, ldwv )
802 CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
803 $ z( jrow, incol+k1 ), ldz )
804 170 CONTINUE
805 END IF
806 END IF
807 180 CONTINUE
808*
809* ==== End of CLAQR5 ====
810*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:177
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
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
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
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: