LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 scalar
             WANTT = .true. if the triangular Schur factor
             is being computed.  WANTT is set to .false. otherwise.
[in]WANTZ
          WANTZ is logical scalar
             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 scalar
             N is the order of the Hessenberg matrix H upon which this
             subroutine operates.
[in]KTOP
          KTOP is integer scalar
[in]KBOT
          KBOT is integer scalar
             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 scalar
             NSHFTS gives the number of simultaneous shifts.  NSHFTS
             must be positive and even.
[in,out]S
          S is COMPLEX array of size (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 of size (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 scalar
             LDH is the leading dimension of H just as declared in the
             calling procedure.  LDH.GE.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 .LE. ILOZ .LE. IHIZ .LE. N
[in,out]Z
          Z is COMPLEX array of size (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 scalar
             LDA is the leading dimension of Z just as declared in
             the calling procedure. LDZ.GE.N.
[out]V
          V is COMPLEX array of size (LDV,NSHFTS/2)
[in]LDV
          LDV is integer scalar
             LDV is the leading dimension of V as declared in the
             calling procedure.  LDV.GE.3.
[out]U
          U is COMPLEX array of size
             (LDU,3*NSHFTS-3)
[in]LDU
          LDU is integer scalar
             LDU is the leading dimension of U just as declared in the
             in the calling subroutine.  LDU.GE.3*NSHFTS-3.
[in]NH
          NH is integer scalar
             NH is the number of columns in array WH available for
             workspace. NH.GE.1.
[out]WH
          WH is COMPLEX array of size (LDWH,NH)
[in]LDWH
          LDWH is integer scalar
             Leading dimension of WH just as declared in the
             calling procedure.  LDWH.GE.3*NSHFTS-3.
[in]NV
          NV is integer scalar
             NV is the number of rows in WV agailable for workspace.
             NV.GE.1.
[out]WV
          WV is COMPLEX array of size
             (LDWV,3*NSHFTS-3)
[in]LDWV
          LDWV is integer scalar
             LDWV is the leading dimension of WV as declared in the
             in the calling subroutine.  LDWV.GE.NV.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 253 of file claqr5.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: