LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ 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.

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.```
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
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,
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  DO 30 j = jtop, min( kbot, k+3 )
428  refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
429  \$ h( j, k+2 ) )
430  h( j, k+1 ) = h( j, k+1 ) - refsum
431  h( j, k+2 ) = h( j, k+2 ) -
432  \$ refsum*conjg( v( 2, m22 ) )
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  DO 40 j = k+1, jbot
446  refsum = conjg( v( 1, m22 ) )*
447  \$ ( h( k+1, j )+conjg( v( 2, m22 ) )*
448  \$ h( k+2, j ) )
449  h( k+1, j ) = h( k+1, j ) - refsum
450  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
451  40 CONTINUE
452 *
453 * ==== The following convergence test requires that
454 * . the tradition small-compared-to-nearby-diagonals
455 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
456 * . criteria both be satisfied. The latter improves
457 * . accuracy in some examples. Falling back on an
458 * . alternate convergence criterion when TST1 or TST2
459 * . is zero (as done here) is traditional but probably
460 * . unnecessary. ====
461 *
462  IF( k.GE.ktop) THEN
463  IF( h( k+1, k ).NE.zero ) THEN
464  tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
465  IF( tst1.EQ.rzero ) THEN
466  IF( k.GE.ktop+1 )
467  \$ tst1 = tst1 + cabs1( h( k, k-1 ) )
468  IF( k.GE.ktop+2 )
469  \$ tst1 = tst1 + cabs1( h( k, k-2 ) )
470  IF( k.GE.ktop+3 )
471  \$ tst1 = tst1 + cabs1( h( k, k-3 ) )
472  IF( k.LE.kbot-2 )
473  \$ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
474  IF( k.LE.kbot-3 )
475  \$ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
476  IF( k.LE.kbot-4 )
477  \$ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
478  END IF
479  IF( cabs1( h( k+1, k ) )
480  \$ .LE.max( smlnum, ulp*tst1 ) ) THEN
481  h12 = max( cabs1( h( k+1, k ) ),
482  \$ cabs1( h( k, k+1 ) ) )
483  h21 = min( cabs1( h( k+1, k ) ),
484  \$ cabs1( h( k, k+1 ) ) )
485  h11 = max( cabs1( h( k+1, k+1 ) ),
486  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
487  h22 = min( cabs1( h( k+1, k+1 ) ),
488  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
489  scl = h11 + h12
490  tst2 = h22*( h11 / scl )
491 *
492  IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
493  \$ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
494  END IF
495  END IF
496  END IF
497 *
498 * ==== Accumulate orthogonal transformations. ====
499 *
500  IF( accum ) THEN
501  kms = k - incol
502  DO 50 j = max( 1, ktop-incol ), kdu
503  refsum = v( 1, m22 )*( u( j, kms+1 )+
504  \$ v( 2, m22 )*u( j, kms+2 ) )
505  u( j, kms+1 ) = u( j, kms+1 ) - refsum
506  u( j, kms+2 ) = u( j, kms+2 ) -
507  \$ refsum*conjg( v( 2, m22 ) )
508  50 CONTINUE
509  ELSE IF( wantz ) THEN
510  DO 60 j = iloz, ihiz
511  refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
512  \$ z( j, k+2 ) )
513  z( j, k+1 ) = z( j, k+1 ) - refsum
514  z( j, k+2 ) = z( j, k+2 ) -
515  \$ refsum*conjg( v( 2, m22 ) )
516  60 CONTINUE
517  END IF
518  END IF
519 *
520 * ==== Normal case: Chain of 3-by-3 reflections ====
521 *
522  DO 80 m = mbot, mtop, -1
523  k = krcol + 2*( m-1 )
524  IF( k.EQ.ktop-1 ) THEN
525  CALL claqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
526  \$ s( 2*m ), v( 1, m ) )
527  alpha = v( 1, m )
528  CALL clarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
529  ELSE
530 *
531 * ==== Perform delayed transformation of row below
532 * . Mth bulge. Exploit fact that first two elements
533 * . of row are actually zero. ====
534 *
535  refsum = v( 1, m )*v( 3, m )*h( k+3, k+2 )
536  h( k+3, k ) = -refsum
537  h( k+3, k+1 ) = -refsum*conjg( v( 2, m ) )
538  h( k+3, k+2 ) = h( k+3, k+2 ) -
539  \$ refsum*conjg( v( 3, m ) )
540 *
541 * ==== Calculate reflection to move
542 * . Mth bulge one step. ====
543 *
544  beta = h( k+1, k )
545  v( 2, m ) = h( k+2, k )
546  v( 3, m ) = h( k+3, k )
547  CALL clarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
548 *
549 * ==== A Bulge may collapse because of vigilant
550 * . deflation or destructive underflow. In the
551 * . underflow case, try the two-small-subdiagonals
552 * . trick to try to reinflate the bulge. ====
553 *
554  IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
555  \$ zero .OR. h( k+3, k+2 ).EQ.zero ) THEN
556 *
557 * ==== Typical case: not collapsed (yet). ====
558 *
559  h( k+1, k ) = beta
560  h( k+2, k ) = zero
561  h( k+3, k ) = zero
562  ELSE
563 *
564 * ==== Atypical case: collapsed. Attempt to
565 * . reintroduce ignoring H(K+1,K) and H(K+2,K).
566 * . If the fill resulting from the new
567 * . reflector is too large, then abandon it.
568 * . Otherwise, use the new one. ====
569 *
570  CALL claqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
571  \$ s( 2*m ), vt )
572  alpha = vt( 1 )
573  CALL clarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
574  refsum = conjg( vt( 1 ) )*
575  \$ ( h( k+1, k )+conjg( vt( 2 ) )*
576  \$ h( k+2, k ) )
577 *
578  IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
579  \$ cabs1( refsum*vt( 3 ) ).GT.ulp*
580  \$ ( cabs1( h( k, k ) )+cabs1( h( k+1,
581  \$ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) ) THEN
582 *
583 * ==== Starting a new bulge here would
584 * . create non-negligible fill. Use
585 * . the old one with trepidation. ====
586 *
587  h( k+1, k ) = beta
588  h( k+2, k ) = zero
589  h( k+3, k ) = zero
590  ELSE
591 *
592 * ==== Starting a new bulge here would
593 * . create only negligible fill.
594 * . Replace the old reflector with
595 * . the new one. ====
596 *
597  h( k+1, k ) = h( k+1, k ) - refsum
598  h( k+2, k ) = zero
599  h( k+3, k ) = zero
600  v( 1, m ) = vt( 1 )
601  v( 2, m ) = vt( 2 )
602  v( 3, m ) = vt( 3 )
603  END IF
604  END IF
605  END IF
606 *
607 * ==== Apply reflection from the right and
608 * . the first column of update from the left.
609 * . These updates are required for the vigilant
610 * . deflation check. We still delay most of the
611 * . updates from the left for efficiency. ====
612 *
613  DO 70 j = jtop, min( kbot, k+3 )
614  refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
615  \$ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
616  h( j, k+1 ) = h( j, k+1 ) - refsum
617  h( j, k+2 ) = h( j, k+2 ) -
618  \$ refsum*conjg( v( 2, m ) )
619  h( j, k+3 ) = h( j, k+3 ) -
620  \$ refsum*conjg( v( 3, m ) )
621  70 CONTINUE
622 *
623 * ==== Perform update from left for subsequent
624 * . column. ====
625 *
626  refsum = conjg( v( 1, m ) )*( h( k+1, k+1 )
627  \$ +conjg( v( 2, m ) )*h( k+2, k+1 )
628  \$ +conjg( v( 3, m ) )*h( k+3, k+1 ) )
629  h( k+1, k+1 ) = h( k+1, k+1 ) - refsum
630  h( k+2, k+1 ) = h( k+2, k+1 ) - refsum*v( 2, m )
631  h( k+3, k+1 ) = h( k+3, k+1 ) - refsum*v( 3, m )
632 *
633 * ==== The following convergence test requires that
634 * . the tradition small-compared-to-nearby-diagonals
635 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
636 * . criteria both be satisfied. The latter improves
637 * . accuracy in some examples. Falling back on an
638 * . alternate convergence criterion when TST1 or TST2
639 * . is zero (as done here) is traditional but probably
640 * . unnecessary. ====
641 *
642  IF( k.LT.ktop)
643  \$ cycle
644  IF( h( k+1, k ).NE.zero ) THEN
645  tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
646  IF( tst1.EQ.rzero ) THEN
647  IF( k.GE.ktop+1 )
648  \$ tst1 = tst1 + cabs1( h( k, k-1 ) )
649  IF( k.GE.ktop+2 )
650  \$ tst1 = tst1 + cabs1( h( k, k-2 ) )
651  IF( k.GE.ktop+3 )
652  \$ tst1 = tst1 + cabs1( h( k, k-3 ) )
653  IF( k.LE.kbot-2 )
654  \$ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
655  IF( k.LE.kbot-3 )
656  \$ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
657  IF( k.LE.kbot-4 )
658  \$ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
659  END IF
660  IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
661  \$ THEN
662  h12 = max( cabs1( h( k+1, k ) ),
663  \$ cabs1( h( k, k+1 ) ) )
664  h21 = min( cabs1( h( k+1, k ) ),
665  \$ cabs1( h( k, k+1 ) ) )
666  h11 = max( cabs1( h( k+1, k+1 ) ),
667  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
668  h22 = min( cabs1( h( k+1, k+1 ) ),
669  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
670  scl = h11 + h12
671  tst2 = h22*( h11 / scl )
672 *
673  IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
674  \$ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
675  END IF
676  END IF
677  80 CONTINUE
678 *
679 * ==== Multiply H by reflections from the left ====
680 *
681  IF( accum ) THEN
682  jbot = min( ndcol, kbot )
683  ELSE IF( wantt ) THEN
684  jbot = n
685  ELSE
686  jbot = kbot
687  END IF
688 *
689  DO 100 m = mbot, mtop, -1
690  k = krcol + 2*( m-1 )
691  DO 90 j = max( ktop, krcol + 2*m ), jbot
692  refsum = conjg( v( 1, m ) )*
693  \$ ( h( k+1, j )+conjg( v( 2, m ) )*
694  \$ h( k+2, j )+conjg( v( 3, m ) )*h( k+3, j ) )
695  h( k+1, j ) = h( k+1, j ) - refsum
696  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
697  h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
698  90 CONTINUE
699  100 CONTINUE
700 *
701 * ==== Accumulate orthogonal transformations. ====
702 *
703  IF( accum ) THEN
704 *
705 * ==== Accumulate U. (If needed, update Z later
706 * . with an efficient matrix-matrix
707 * . multiply.) ====
708 *
709  DO 120 m = mbot, mtop, -1
710  k = krcol + 2*( m-1 )
711  kms = k - incol
712  i2 = max( 1, ktop-incol )
713  i2 = max( i2, kms-(krcol-incol)+1 )
714  i4 = min( kdu, krcol + 2*( mbot-1 ) - incol + 5 )
715  DO 110 j = i2, i4
716  refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
717  \$ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
718  u( j, kms+1 ) = u( j, kms+1 ) - refsum
719  u( j, kms+2 ) = u( j, kms+2 ) -
720  \$ refsum*conjg( v( 2, m ) )
721  u( j, kms+3 ) = u( j, kms+3 ) -
722  \$ refsum*conjg( v( 3, m ) )
723  110 CONTINUE
724  120 CONTINUE
725  ELSE IF( wantz ) THEN
726 *
727 * ==== U is not accumulated, so update Z
728 * . now by multiplying by reflections
729 * . from the right. ====
730 *
731  DO 140 m = mbot, mtop, -1
732  k = krcol + 2*( m-1 )
733  DO 130 j = iloz, ihiz
734  refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
735  \$ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
736  z( j, k+1 ) = z( j, k+1 ) - refsum
737  z( j, k+2 ) = z( j, k+2 ) -
738  \$ refsum*conjg( v( 2, m ) )
739  z( j, k+3 ) = z( j, k+3 ) -
740  \$ refsum*conjg( v( 3, m ) )
741  130 CONTINUE
742  140 CONTINUE
743  END IF
744 *
745 * ==== End of near-the-diagonal bulge chase. ====
746 *
747  145 CONTINUE
748 *
749 * ==== Use U (if accumulated) to update far-from-diagonal
750 * . entries in H. If required, use U to update Z as
751 * . well. ====
752 *
753  IF( accum ) THEN
754  IF( wantt ) THEN
755  jtop = 1
756  jbot = n
757  ELSE
758  jtop = ktop
759  jbot = kbot
760  END IF
761  k1 = max( 1, ktop-incol )
762  nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
763 *
764 * ==== Horizontal Multiply ====
765 *
766  DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
767  jlen = min( nh, jbot-jcol+1 )
768  CALL cgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
769  \$ ldu, h( incol+k1, jcol ), ldh, zero, wh,
770  \$ ldwh )
771  CALL clacpy( 'ALL', nu, jlen, wh, ldwh,
772  \$ h( incol+k1, jcol ), ldh )
773  150 CONTINUE
774 *
775 * ==== Vertical multiply ====
776 *
777  DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
778  jlen = min( nv, max( ktop, incol )-jrow )
779  CALL cgemm( 'N', 'N', jlen, nu, nu, one,
780  \$ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
781  \$ ldu, zero, wv, ldwv )
782  CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
783  \$ h( jrow, incol+k1 ), ldh )
784  160 CONTINUE
785 *
786 * ==== Z multiply (also vertical) ====
787 *
788  IF( wantz ) THEN
789  DO 170 jrow = iloz, ihiz, nv
790  jlen = min( nv, ihiz-jrow+1 )
791  CALL cgemm( 'N', 'N', jlen, nu, nu, one,
792  \$ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
793  \$ ldu, zero, wv, ldwv )
794  CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
795  \$ z( jrow, incol+k1 ), ldz )
796  170 CONTINUE
797  END IF
798  END IF
799  180 CONTINUE
800 *
801 * ==== End of CLAQR5 ====
802 *
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: