LAPACK  3.9.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: CLAQR5 accumulates reflections, uses matrix-matrix multiply to update the far-from-diagonal matrix entries, and takes advantage of 2-by-2 block structure during matrix multiplies.``` [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,3*NSHFTS-3)` [in] LDU ``` LDU is INTEGER LDU is the leading dimension of U just as declared in the in the calling subroutine. LDU >= 3*NSHFTS-3.``` [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,3*NSHFTS-3)` [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 >= 3*NSHFTS-3.```
Date
June 2016
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
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.

Definition at line 251 of file claqr5.f.

251 *
252 * -- LAPACK auxiliary routine (version 3.7.1) --
253 * -- LAPACK is a software package provided by Univ. of Tennessee, --
254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255 * June 2016
256 *
257 * .. Scalar Arguments ..
258  INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
259  \$ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
260  LOGICAL WANTT, WANTZ
261 * ..
262 * .. Array Arguments ..
263  COMPLEX H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
264  \$ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
265 * ..
266 *
267 * ================================================================
268 * .. Parameters ..
269  COMPLEX ZERO, ONE
270  parameter( zero = ( 0.0e0, 0.0e0 ),
271  \$ one = ( 1.0e0, 0.0e0 ) )
272  REAL RZERO, RONE
273  parameter( rzero = 0.0e0, rone = 1.0e0 )
274 * ..
275 * .. Local Scalars ..
276  COMPLEX ALPHA, BETA, CDUM, REFSUM
277  REAL H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
278  \$ SMLNUM, TST1, TST2, ULP
279  INTEGER I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
280  \$ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
281  \$ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
282  \$ NS, NU
283  LOGICAL ACCUM, BLK22, BMP22
284 * ..
285 * .. External Functions ..
286  REAL SLAMCH
287  EXTERNAL slamch
288 * ..
289 * .. Intrinsic Functions ..
290 *
291  INTRINSIC abs, aimag, conjg, max, min, mod, real
292 * ..
293 * .. Local Arrays ..
294  COMPLEX VT( 3 )
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL cgemm, clacpy, claqr1, clarfg, claset, ctrmm,
299 * ..
300 * .. Statement Functions ..
301  REAL CABS1
302 * ..
303 * .. Statement Function definitions ..
304  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
305 * ..
306 * .. Executable Statements ..
307 *
308 * ==== If there are no shifts, then there is nothing to do. ====
309 *
310  IF( nshfts.LT.2 )
311  \$ RETURN
312 *
313 * ==== If the active block is empty or 1-by-1, then there
314 * . is nothing to do. ====
315 *
316  IF( ktop.GE.kbot )
317  \$ RETURN
318 *
319 * ==== NSHFTS is supposed to be even, but if it is odd,
320 * . then simply reduce it by one. ====
321 *
322  ns = nshfts - mod( nshfts, 2 )
323 *
324 * ==== Machine constants for deflation ====
325 *
326  safmin = slamch( 'SAFE MINIMUM' )
327  safmax = rone / safmin
328  CALL slabad( safmin, safmax )
329  ulp = slamch( 'PRECISION' )
330  smlnum = safmin*( real( n ) / ulp )
331 *
332 * ==== Use accumulated reflections to update far-from-diagonal
333 * . entries ? ====
334 *
335  accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
336 *
337 * ==== If so, exploit the 2-by-2 block structure? ====
338 *
339  blk22 = ( ns.GT.2 ) .AND. ( 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 = 6*nbmps - 3
353 *
354 * ==== Create and chase chains of NBMPS bulges ====
355 *
356  DO 210 incol = 3*( 1-nbmps ) + ktop - 1, kbot - 2, 3*nbmps - 2
357  ndcol = incol + kdu
358  IF( accum )
359  \$ CALL claset( 'ALL', kdu, kdu, zero, one, u, ldu )
360 *
361 * ==== Near-the-diagonal bulge chase. The following loop
362 * . performs the near-the-diagonal part of a small bulge
363 * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
364 * . chunk extends from column INCOL to column NDCOL
365 * . (including both column INCOL and column NDCOL). The
366 * . following loop chases a 3*NBMPS column long chain of
367 * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
368 * . may be less than KTOP and and NDCOL may be greater than
369 * . KBOT indicating phantom columns from which to chase
370 * . bulges before they are actually introduced or to which
371 * . to chase bulges beyond column KBOT.) ====
372 *
373  DO 140 krcol = incol, min( incol+3*nbmps-3, kbot-2 )
374 *
375 * ==== Bulges number MTOP to MBOT are active double implicit
376 * . shift bulges. There may or may not also be small
377 * . 2-by-2 bulge, if there is room. The inactive bulges
378 * . (if any) must wait until the active bulges have moved
379 * . down the diagonal to make room. The phantom matrix
380 * . paradigm described above helps keep track. ====
381 *
382  mtop = max( 1, ( ( ktop-1 )-krcol+2 ) / 3+1 )
383  mbot = min( nbmps, ( kbot-krcol ) / 3 )
384  m22 = mbot + 1
385  bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+3*( m22-1 ) ).EQ.
386  \$ ( kbot-2 )
387 *
388 * ==== Generate reflections to chase the chain right
389 * . one column. (The minimum value of K is KTOP-1.) ====
390 *
391  DO 10 m = mtop, mbot
392  k = krcol + 3*( m-1 )
393  IF( k.EQ.ktop-1 ) THEN
394  CALL claqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
395  \$ s( 2*m ), v( 1, m ) )
396  alpha = v( 1, m )
397  CALL clarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
398  ELSE
399  beta = h( k+1, k )
400  v( 2, m ) = h( k+2, k )
401  v( 3, m ) = h( k+3, k )
402  CALL clarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
403 *
404 * ==== A Bulge may collapse because of vigilant
405 * . deflation or destructive underflow. In the
406 * . underflow case, try the two-small-subdiagonals
407 * . trick to try to reinflate the bulge. ====
408 *
409  IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
410  \$ zero .OR. h( k+3, k+2 ).EQ.zero ) THEN
411 *
412 * ==== Typical case: not collapsed (yet). ====
413 *
414  h( k+1, k ) = beta
415  h( k+2, k ) = zero
416  h( k+3, k ) = zero
417  ELSE
418 *
419 * ==== Atypical case: collapsed. Attempt to
420 * . reintroduce ignoring H(K+1,K) and H(K+2,K).
421 * . If the fill resulting from the new
422 * . reflector is too large, then abandon it.
423 * . Otherwise, use the new one. ====
424 *
425  CALL claqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
426  \$ s( 2*m ), vt )
427  alpha = vt( 1 )
428  CALL clarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
429  refsum = conjg( vt( 1 ) )*
430  \$ ( h( k+1, k )+conjg( vt( 2 ) )*
431  \$ h( k+2, k ) )
432 *
433  IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
434  \$ cabs1( refsum*vt( 3 ) ).GT.ulp*
435  \$ ( cabs1( h( k, k ) )+cabs1( h( k+1,
436  \$ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) ) THEN
437 *
438 * ==== Starting a new bulge here would
439 * . create non-negligible fill. Use
440 * . the old one with trepidation. ====
441 *
442  h( k+1, k ) = beta
443  h( k+2, k ) = zero
444  h( k+3, k ) = zero
445  ELSE
446 *
447 * ==== Stating a new bulge here would
448 * . create only negligible fill.
449 * . Replace the old reflector with
450 * . the new one. ====
451 *
452  h( k+1, k ) = h( k+1, k ) - refsum
453  h( k+2, k ) = zero
454  h( k+3, k ) = zero
455  v( 1, m ) = vt( 1 )
456  v( 2, m ) = vt( 2 )
457  v( 3, m ) = vt( 3 )
458  END IF
459  END IF
460  END IF
461  10 CONTINUE
462 *
463 * ==== Generate a 2-by-2 reflection, if needed. ====
464 *
465  k = krcol + 3*( m22-1 )
466  IF( bmp22 ) THEN
467  IF( k.EQ.ktop-1 ) THEN
468  CALL claqr1( 2, h( k+1, k+1 ), ldh, s( 2*m22-1 ),
469  \$ s( 2*m22 ), v( 1, m22 ) )
470  beta = v( 1, m22 )
471  CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
472  ELSE
473  beta = h( k+1, k )
474  v( 2, m22 ) = h( k+2, k )
475  CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
476  h( k+1, k ) = beta
477  h( k+2, k ) = zero
478  END IF
479  END IF
480 *
481 * ==== Multiply H by reflections from the left ====
482 *
483  IF( accum ) THEN
484  jbot = min( ndcol, kbot )
485  ELSE IF( wantt ) THEN
486  jbot = n
487  ELSE
488  jbot = kbot
489  END IF
490  DO 30 j = max( ktop, krcol ), jbot
491  mend = min( mbot, ( j-krcol+2 ) / 3 )
492  DO 20 m = mtop, mend
493  k = krcol + 3*( m-1 )
494  refsum = conjg( v( 1, m ) )*
495  \$ ( h( k+1, j )+conjg( v( 2, m ) )*h( k+2, j )+
496  \$ conjg( v( 3, m ) )*h( k+3, j ) )
497  h( k+1, j ) = h( k+1, j ) - refsum
498  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
499  h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
500  20 CONTINUE
501  30 CONTINUE
502  IF( bmp22 ) THEN
503  k = krcol + 3*( m22-1 )
504  DO 40 j = max( k+1, ktop ), jbot
505  refsum = conjg( v( 1, m22 ) )*
506  \$ ( h( k+1, j )+conjg( v( 2, m22 ) )*
507  \$ h( k+2, j ) )
508  h( k+1, j ) = h( k+1, j ) - refsum
509  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
510  40 CONTINUE
511  END IF
512 *
513 * ==== Multiply H by reflections from the right.
514 * . Delay filling in the last row until the
515 * . vigilant deflation check is complete. ====
516 *
517  IF( accum ) THEN
518  jtop = max( ktop, incol )
519  ELSE IF( wantt ) THEN
520  jtop = 1
521  ELSE
522  jtop = ktop
523  END IF
524  DO 80 m = mtop, mbot
525  IF( v( 1, m ).NE.zero ) THEN
526  k = krcol + 3*( m-1 )
527  DO 50 j = jtop, min( kbot, k+3 )
528  refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
529  \$ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
530  h( j, k+1 ) = h( j, k+1 ) - refsum
531  h( j, k+2 ) = h( j, k+2 ) -
532  \$ refsum*conjg( v( 2, m ) )
533  h( j, k+3 ) = h( j, k+3 ) -
534  \$ refsum*conjg( v( 3, m ) )
535  50 CONTINUE
536 *
537  IF( accum ) THEN
538 *
539 * ==== Accumulate U. (If necessary, update Z later
540 * . with with an efficient matrix-matrix
541 * . multiply.) ====
542 *
543  kms = k - incol
544  DO 60 j = max( 1, ktop-incol ), kdu
545  refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
546  \$ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
547  u( j, kms+1 ) = u( j, kms+1 ) - refsum
548  u( j, kms+2 ) = u( j, kms+2 ) -
549  \$ refsum*conjg( v( 2, m ) )
550  u( j, kms+3 ) = u( j, kms+3 ) -
551  \$ refsum*conjg( v( 3, m ) )
552  60 CONTINUE
553  ELSE IF( wantz ) THEN
554 *
555 * ==== U is not accumulated, so update Z
556 * . now by multiplying by reflections
557 * . from the right. ====
558 *
559  DO 70 j = iloz, ihiz
560  refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
561  \$ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
562  z( j, k+1 ) = z( j, k+1 ) - refsum
563  z( j, k+2 ) = z( j, k+2 ) -
564  \$ refsum*conjg( v( 2, m ) )
565  z( j, k+3 ) = z( j, k+3 ) -
566  \$ refsum*conjg( v( 3, m ) )
567  70 CONTINUE
568  END IF
569  END IF
570  80 CONTINUE
571 *
572 * ==== Special case: 2-by-2 reflection (if needed) ====
573 *
574  k = krcol + 3*( m22-1 )
575  IF( bmp22 ) THEN
576  IF ( v( 1, m22 ).NE.zero ) THEN
577  DO 90 j = jtop, min( kbot, k+3 )
578  refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
579  \$ h( j, k+2 ) )
580  h( j, k+1 ) = h( j, k+1 ) - refsum
581  h( j, k+2 ) = h( j, k+2 ) -
582  \$ refsum*conjg( v( 2, m22 ) )
583  90 CONTINUE
584 *
585  IF( accum ) THEN
586  kms = k - incol
587  DO 100 j = max( 1, ktop-incol ), kdu
588  refsum = v( 1, m22 )*( u( j, kms+1 )+
589  \$ v( 2, m22 )*u( j, kms+2 ) )
590  u( j, kms+1 ) = u( j, kms+1 ) - refsum
591  u( j, kms+2 ) = u( j, kms+2 ) -
592  \$ refsum*conjg( v( 2, m22 ) )
593  100 CONTINUE
594  ELSE IF( wantz ) THEN
595  DO 110 j = iloz, ihiz
596  refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
597  \$ z( j, k+2 ) )
598  z( j, k+1 ) = z( j, k+1 ) - refsum
599  z( j, k+2 ) = z( j, k+2 ) -
600  \$ refsum*conjg( v( 2, m22 ) )
601  110 CONTINUE
602  END IF
603  END IF
604  END IF
605 *
606 * ==== Vigilant deflation check ====
607 *
608  mstart = mtop
609  IF( krcol+3*( mstart-1 ).LT.ktop )
610  \$ mstart = mstart + 1
611  mend = mbot
612  IF( bmp22 )
613  \$ mend = mend + 1
614  IF( krcol.EQ.kbot-2 )
615  \$ mend = mend + 1
616  DO 120 m = mstart, mend
617  k = min( kbot-1, krcol+3*( m-1 ) )
618 *
619 * ==== The following convergence test requires that
620 * . the tradition small-compared-to-nearby-diagonals
621 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
622 * . criteria both be satisfied. The latter improves
623 * . accuracy in some examples. Falling back on an
624 * . alternate convergence criterion when TST1 or TST2
625 * . is zero (as done here) is traditional but probably
626 * . unnecessary. ====
627 *
628  IF( h( k+1, k ).NE.zero ) THEN
629  tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
630  IF( tst1.EQ.rzero ) THEN
631  IF( k.GE.ktop+1 )
632  \$ tst1 = tst1 + cabs1( h( k, k-1 ) )
633  IF( k.GE.ktop+2 )
634  \$ tst1 = tst1 + cabs1( h( k, k-2 ) )
635  IF( k.GE.ktop+3 )
636  \$ tst1 = tst1 + cabs1( h( k, k-3 ) )
637  IF( k.LE.kbot-2 )
638  \$ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
639  IF( k.LE.kbot-3 )
640  \$ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
641  IF( k.LE.kbot-4 )
642  \$ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
643  END IF
644  IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
645  \$ THEN
646  h12 = max( cabs1( h( k+1, k ) ),
647  \$ cabs1( h( k, k+1 ) ) )
648  h21 = min( cabs1( h( k+1, k ) ),
649  \$ cabs1( h( k, k+1 ) ) )
650  h11 = max( cabs1( h( k+1, k+1 ) ),
651  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
652  h22 = min( cabs1( h( k+1, k+1 ) ),
653  \$ cabs1( h( k, k )-h( k+1, k+1 ) ) )
654  scl = h11 + h12
655  tst2 = h22*( h11 / scl )
656 *
657  IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
658  \$ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
659  END IF
660  END IF
661  120 CONTINUE
662 *
663 * ==== Fill in the last row of each bulge. ====
664 *
665  mend = min( nbmps, ( kbot-krcol-1 ) / 3 )
666  DO 130 m = mtop, mend
667  k = krcol + 3*( m-1 )
668  refsum = v( 1, m )*v( 3, m )*h( k+4, k+3 )
669  h( k+4, k+1 ) = -refsum
670  h( k+4, k+2 ) = -refsum*conjg( v( 2, m ) )
671  h( k+4, k+3 ) = h( k+4, k+3 ) - refsum*conjg( v( 3, m ) )
672  130 CONTINUE
673 *
674 * ==== End of near-the-diagonal bulge chase. ====
675 *
676  140 CONTINUE
677 *
678 * ==== Use U (if accumulated) to update far-from-diagonal
679 * . entries in H. If required, use U to update Z as
680 * . well. ====
681 *
682  IF( accum ) THEN
683  IF( wantt ) THEN
684  jtop = 1
685  jbot = n
686  ELSE
687  jtop = ktop
688  jbot = kbot
689  END IF
690  IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
691  \$ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) ) THEN
692 *
693 * ==== Updates not exploiting the 2-by-2 block
694 * . structure of U. K1 and NU keep track of
695 * . the location and size of U in the special
696 * . cases of introducing bulges and chasing
697 * . bulges off the bottom. In these special
698 * . cases and in case the number of shifts
699 * . is NS = 2, there is no 2-by-2 block
700 * . structure to exploit. ====
701 *
702  k1 = max( 1, ktop-incol )
703  nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
704 *
705 * ==== Horizontal Multiply ====
706 *
707  DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
708  jlen = min( nh, jbot-jcol+1 )
709  CALL cgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
710  \$ ldu, h( incol+k1, jcol ), ldh, zero, wh,
711  \$ ldwh )
712  CALL clacpy( 'ALL', nu, jlen, wh, ldwh,
713  \$ h( incol+k1, jcol ), ldh )
714  150 CONTINUE
715 *
716 * ==== Vertical multiply ====
717 *
718  DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
719  jlen = min( nv, max( ktop, incol )-jrow )
720  CALL cgemm( 'N', 'N', jlen, nu, nu, one,
721  \$ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
722  \$ ldu, zero, wv, ldwv )
723  CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
724  \$ h( jrow, incol+k1 ), ldh )
725  160 CONTINUE
726 *
727 * ==== Z multiply (also vertical) ====
728 *
729  IF( wantz ) THEN
730  DO 170 jrow = iloz, ihiz, nv
731  jlen = min( nv, ihiz-jrow+1 )
732  CALL cgemm( 'N', 'N', jlen, nu, nu, one,
733  \$ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
734  \$ ldu, zero, wv, ldwv )
735  CALL clacpy( 'ALL', jlen, nu, wv, ldwv,
736  \$ z( jrow, incol+k1 ), ldz )
737  170 CONTINUE
738  END IF
739  ELSE
740 *
741 * ==== Updates exploiting U's 2-by-2 block structure.
742 * . (I2, I4, J2, J4 are the last rows and columns
743 * . of the blocks.) ====
744 *
745  i2 = ( kdu+1 ) / 2
746  i4 = kdu
747  j2 = i4 - i2
748  j4 = kdu
749 *
750 * ==== KZS and KNZ deal with the band of zeros
751 * . along the diagonal of one of the triangular
752 * . blocks. ====
753 *
754  kzs = ( j4-j2 ) - ( ns+1 )
755  knz = ns + 1
756 *
757 * ==== Horizontal multiply ====
758 *
759  DO 180 jcol = min( ndcol, kbot ) + 1, jbot, nh
760  jlen = min( nh, jbot-jcol+1 )
761 *
762 * ==== Copy bottom of H to top+KZS of scratch ====
763 * (The first KZS rows get multiplied by zero.) ====
764 *
765  CALL clacpy( 'ALL', knz, jlen, h( incol+1+j2, jcol ),
766  \$ ldh, wh( kzs+1, 1 ), ldwh )
767 *
768 * ==== Multiply by U21**H ====
769 *
770  CALL claset( 'ALL', kzs, jlen, zero, zero, wh, ldwh )
771  CALL ctrmm( 'L', 'U', 'C', 'N', knz, jlen, one,
772  \$ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
773  \$ ldwh )
774 *
775 * ==== Multiply top of H by U11**H ====
776 *
777  CALL cgemm( 'C', 'N', i2, jlen, j2, one, u, ldu,
778  \$ h( incol+1, jcol ), ldh, one, wh, ldwh )
779 *
780 * ==== Copy top of H to bottom of WH ====
781 *
782  CALL clacpy( 'ALL', j2, jlen, h( incol+1, jcol ), ldh,
783  \$ wh( i2+1, 1 ), ldwh )
784 *
785 * ==== Multiply by U21**H ====
786 *
787  CALL ctrmm( 'L', 'L', 'C', 'N', j2, jlen, one,
788  \$ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
789 *
790 * ==== Multiply by U22 ====
791 *
792  CALL cgemm( 'C', 'N', i4-i2, jlen, j4-j2, one,
793  \$ u( j2+1, i2+1 ), ldu,
794  \$ h( incol+1+j2, jcol ), ldh, one,
795  \$ wh( i2+1, 1 ), ldwh )
796 *
797 * ==== Copy it back ====
798 *
799  CALL clacpy( 'ALL', kdu, jlen, wh, ldwh,
800  \$ h( incol+1, jcol ), ldh )
801  180 CONTINUE
802 *
803 * ==== Vertical multiply ====
804 *
805  DO 190 jrow = jtop, max( incol, ktop ) - 1, nv
806  jlen = min( nv, max( incol, ktop )-jrow )
807 *
808 * ==== Copy right of H to scratch (the first KZS
809 * . columns get multiplied by zero) ====
810 *
811  CALL clacpy( 'ALL', jlen, knz, h( jrow, incol+1+j2 ),
812  \$ ldh, wv( 1, 1+kzs ), ldwv )
813 *
814 * ==== Multiply by U21 ====
815 *
816  CALL claset( 'ALL', jlen, kzs, zero, zero, wv, ldwv )
817  CALL ctrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
818  \$ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
819  \$ ldwv )
820 *
821 * ==== Multiply by U11 ====
822 *
823  CALL cgemm( 'N', 'N', jlen, i2, j2, one,
824  \$ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
825  \$ ldwv )
826 *
827 * ==== Copy left of H to right of scratch ====
828 *
829  CALL clacpy( 'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
830  \$ wv( 1, 1+i2 ), ldwv )
831 *
832 * ==== Multiply by U21 ====
833 *
834  CALL ctrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
835  \$ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
836 *
837 * ==== Multiply by U22 ====
838 *
839  CALL cgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
840  \$ h( jrow, incol+1+j2 ), ldh,
841  \$ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
842  \$ ldwv )
843 *
844 * ==== Copy it back ====
845 *
846  CALL clacpy( 'ALL', jlen, kdu, wv, ldwv,
847  \$ h( jrow, incol+1 ), ldh )
848  190 CONTINUE
849 *
850 * ==== Multiply Z (also vertical) ====
851 *
852  IF( wantz ) THEN
853  DO 200 jrow = iloz, ihiz, nv
854  jlen = min( nv, ihiz-jrow+1 )
855 *
856 * ==== Copy right of Z to left of scratch (first
857 * . KZS columns get multiplied by zero) ====
858 *
859  CALL clacpy( 'ALL', jlen, knz,
860  \$ z( jrow, incol+1+j2 ), ldz,
861  \$ wv( 1, 1+kzs ), ldwv )
862 *
863 * ==== Multiply by U12 ====
864 *
865  CALL claset( 'ALL', jlen, kzs, zero, zero, wv,
866  \$ ldwv )
867  CALL ctrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
868  \$ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
869  \$ ldwv )
870 *
871 * ==== Multiply by U11 ====
872 *
873  CALL cgemm( 'N', 'N', jlen, i2, j2, one,
874  \$ z( jrow, incol+1 ), ldz, u, ldu, one,
875  \$ wv, ldwv )
876 *
877 * ==== Copy left of Z to right of scratch ====
878 *
879  CALL clacpy( 'ALL', jlen, j2, z( jrow, incol+1 ),
880  \$ ldz, wv( 1, 1+i2 ), ldwv )
881 *
882 * ==== Multiply by U21 ====
883 *
884  CALL ctrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
885  \$ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
886  \$ ldwv )
887 *
888 * ==== Multiply by U22 ====
889 *
890  CALL cgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
891  \$ z( jrow, incol+1+j2 ), ldz,
892  \$ u( j2+1, i2+1 ), ldu, one,
893  \$ wv( 1, 1+i2 ), ldwv )
894 *
895 * ==== Copy the result back to Z ====
896 *
897  CALL clacpy( 'ALL', jlen, kdu, wv, ldwv,
898  \$ z( jrow, incol+1 ), ldz )
899  200 CONTINUE
900  END IF
901  END IF
902  END IF
903  210 CONTINUE
904 *
905 * ==== End of CLAQR5 ====
906 *
Here is the call graph for this function:
Here is the caller graph for this function:
clarfg
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
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:109
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:105
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:108
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
ctrmm
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179