LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ zlaqr5()

 subroutine zlaqr5 ( logical WANTT, logical WANTZ, integer KACC22, integer N, integer KTOP, integer KBOT, integer NSHFTS, complex*16, dimension( * ) S, complex*16, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16, dimension( ldv, * ) V, integer LDV, complex*16, dimension( ldu, * ) U, integer LDU, integer NV, complex*16, dimension( ldwv, * ) WV, integer LDWV, integer NH, complex*16, dimension( ldwh, * ) WH, integer LDWH )

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

Purpose:
```    ZLAQR5, called by ZLAQR0, 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: ZLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries. = 1: ZLAQR5 accumulates reflections and uses matrix-matrix multiply to update the far-from-diagonal matrix entries. = 2: ZLAQR5 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*16 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*16 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.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*16 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.GE.N.``` [out] V ` V is COMPLEX*16 array, dimension (LDV,NSHFTS/2)` [in] LDV ``` LDV is INTEGER LDV is the leading dimension of V as declared in the calling procedure. LDV.GE.3.``` [out] U ` U is COMPLEX*16 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.GE.3*NSHFTS-3.``` [in] NH ``` NH is INTEGER NH is the number of columns in array WH available for workspace. NH.GE.1.``` [out] WH ` WH is COMPLEX*16 array, dimension (LDWH,NH)` [in] LDWH ``` LDWH is INTEGER Leading dimension of WH just as declared in the calling procedure. LDWH.GE.3*NSHFTS-3.``` [in] NV ``` NV is INTEGER NV is the number of rows in WV agailable for workspace. NV.GE.1.``` [out] WV ` WV is COMPLEX*16 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.GE.NV.```
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 zlaqr5.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*16 h( ldh, * ), s( * ), u( ldu, * ), v( ldv, * ),
264  \$ wh( ldwh, * ), wv( ldwv, * ), z( ldz, * )
265 * ..
266 *
267 * ================================================================
268 * .. Parameters ..
269  COMPLEX*16 zero, one
270  parameter( zero = ( 0.0d0, 0.0d0 ),
271  \$ one = ( 1.0d0, 0.0d0 ) )
272  DOUBLE PRECISION rzero, rone
273  parameter( rzero = 0.0d0, rone = 1.0d0 )
274 * ..
275 * .. Local Scalars ..
276  COMPLEX*16 alpha, beta, cdum, refsum
277  DOUBLE PRECISION 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  DOUBLE PRECISION dlamch
287  EXTERNAL dlamch
288 * ..
289 * .. Intrinsic Functions ..
290 *
291  INTRINSIC abs, dble, dconjg, dimag, max, min, mod
292 * ..
293 * .. Local Arrays ..
294  COMPLEX*16 vt( 3 )
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dlabad, zgemm, zlacpy, zlaqr1, zlarfg, zlaset,
298  \$ ztrmm
299 * ..
300 * .. Statement Functions ..
301  DOUBLE PRECISION cabs1
302 * ..
303 * .. Statement Function definitions ..
304  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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 = dlamch( 'SAFE MINIMUM' )
327  safmax = rone / safmin
328  CALL dlabad( safmin, safmax )
329  ulp = dlamch( 'PRECISION' )
330  smlnum = safmin*( dble( 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 zlaset( '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 zlaqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
395  \$ s( 2*m ), v( 1, m ) )
396  alpha = v( 1, m )
397  CALL zlarfg( 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 zlarfg( 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 zlaqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
426  \$ s( 2*m ), vt )
427  alpha = vt( 1 )
428  CALL zlarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
429  refsum = dconjg( vt( 1 ) )*
430  \$ ( h( k+1, k )+dconjg( 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 zlaqr1( 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 zlarfg( 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 zlarfg( 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 = dconjg( v( 1, m ) )*
495  \$ ( h( k+1, j )+dconjg( v( 2, m ) )*
496  \$ h( k+2, j )+dconjg( 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 = dconjg( v( 1, m22 ) )*
506  \$ ( h( k+1, j )+dconjg( 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*dconjg( v( 2, m ) )
533  h( j, k+3 ) = h( j, k+3 ) -
534  \$ refsum*dconjg( 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*dconjg( v( 2, m ) )
550  u( j, kms+3 ) = u( j, kms+3 ) -
551  \$ refsum*dconjg( 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*dconjg( v( 2, m ) )
565  z( j, k+3 ) = z( j, k+3 ) -
566  \$ refsum*dconjg( 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*dconjg( 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*dconjg( 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*dconjg( 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*dconjg( v( 2, m ) )
671  h( k+4, k+3 ) = h( k+4, k+3 ) -
672  \$ refsum*dconjg( v( 3, m ) )
673  130 CONTINUE
674 *
675 * ==== End of near-the-diagonal bulge chase. ====
676 *
677  140 CONTINUE
678 *
679 * ==== Use U (if accumulated) to update far-from-diagonal
680 * . entries in H. If required, use U to update Z as
681 * . well. ====
682 *
683  IF( accum ) THEN
684  IF( wantt ) THEN
685  jtop = 1
686  jbot = n
687  ELSE
688  jtop = ktop
689  jbot = kbot
690  END IF
691  IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
692  \$ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) ) THEN
693 *
694 * ==== Updates not exploiting the 2-by-2 block
695 * . structure of U. K1 and NU keep track of
696 * . the location and size of U in the special
697 * . cases of introducing bulges and chasing
698 * . bulges off the bottom. In these special
699 * . cases and in case the number of shifts
700 * . is NS = 2, there is no 2-by-2 block
701 * . structure to exploit. ====
702 *
703  k1 = max( 1, ktop-incol )
704  nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
705 *
706 * ==== Horizontal Multiply ====
707 *
708  DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
709  jlen = min( nh, jbot-jcol+1 )
710  CALL zgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
711  \$ ldu, h( incol+k1, jcol ), ldh, zero, wh,
712  \$ ldwh )
713  CALL zlacpy( 'ALL', nu, jlen, wh, ldwh,
714  \$ h( incol+k1, jcol ), ldh )
715  150 CONTINUE
716 *
717 * ==== Vertical multiply ====
718 *
719  DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
720  jlen = min( nv, max( ktop, incol )-jrow )
721  CALL zgemm( 'N', 'N', jlen, nu, nu, one,
722  \$ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
723  \$ ldu, zero, wv, ldwv )
724  CALL zlacpy( 'ALL', jlen, nu, wv, ldwv,
725  \$ h( jrow, incol+k1 ), ldh )
726  160 CONTINUE
727 *
728 * ==== Z multiply (also vertical) ====
729 *
730  IF( wantz ) THEN
731  DO 170 jrow = iloz, ihiz, nv
732  jlen = min( nv, ihiz-jrow+1 )
733  CALL zgemm( 'N', 'N', jlen, nu, nu, one,
734  \$ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
735  \$ ldu, zero, wv, ldwv )
736  CALL zlacpy( 'ALL', jlen, nu, wv, ldwv,
737  \$ z( jrow, incol+k1 ), ldz )
738  170 CONTINUE
739  END IF
740  ELSE
741 *
742 * ==== Updates exploiting U's 2-by-2 block structure.
743 * . (I2, I4, J2, J4 are the last rows and columns
744 * . of the blocks.) ====
745 *
746  i2 = ( kdu+1 ) / 2
747  i4 = kdu
748  j2 = i4 - i2
749  j4 = kdu
750 *
751 * ==== KZS and KNZ deal with the band of zeros
752 * . along the diagonal of one of the triangular
753 * . blocks. ====
754 *
755  kzs = ( j4-j2 ) - ( ns+1 )
756  knz = ns + 1
757 *
758 * ==== Horizontal multiply ====
759 *
760  DO 180 jcol = min( ndcol, kbot ) + 1, jbot, nh
761  jlen = min( nh, jbot-jcol+1 )
762 *
763 * ==== Copy bottom of H to top+KZS of scratch ====
764 * (The first KZS rows get multiplied by zero.) ====
765 *
766  CALL zlacpy( 'ALL', knz, jlen, h( incol+1+j2, jcol ),
767  \$ ldh, wh( kzs+1, 1 ), ldwh )
768 *
769 * ==== Multiply by U21**H ====
770 *
771  CALL zlaset( 'ALL', kzs, jlen, zero, zero, wh, ldwh )
772  CALL ztrmm( 'L', 'U', 'C', 'N', knz, jlen, one,
773  \$ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
774  \$ ldwh )
775 *
776 * ==== Multiply top of H by U11**H ====
777 *
778  CALL zgemm( 'C', 'N', i2, jlen, j2, one, u, ldu,
779  \$ h( incol+1, jcol ), ldh, one, wh, ldwh )
780 *
781 * ==== Copy top of H to bottom of WH ====
782 *
783  CALL zlacpy( 'ALL', j2, jlen, h( incol+1, jcol ), ldh,
784  \$ wh( i2+1, 1 ), ldwh )
785 *
786 * ==== Multiply by U21**H ====
787 *
788  CALL ztrmm( 'L', 'L', 'C', 'N', j2, jlen, one,
789  \$ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
790 *
791 * ==== Multiply by U22 ====
792 *
793  CALL zgemm( 'C', 'N', i4-i2, jlen, j4-j2, one,
794  \$ u( j2+1, i2+1 ), ldu,
795  \$ h( incol+1+j2, jcol ), ldh, one,
796  \$ wh( i2+1, 1 ), ldwh )
797 *
798 * ==== Copy it back ====
799 *
800  CALL zlacpy( 'ALL', kdu, jlen, wh, ldwh,
801  \$ h( incol+1, jcol ), ldh )
802  180 CONTINUE
803 *
804 * ==== Vertical multiply ====
805 *
806  DO 190 jrow = jtop, max( incol, ktop ) - 1, nv
807  jlen = min( nv, max( incol, ktop )-jrow )
808 *
809 * ==== Copy right of H to scratch (the first KZS
810 * . columns get multiplied by zero) ====
811 *
812  CALL zlacpy( 'ALL', jlen, knz, h( jrow, incol+1+j2 ),
813  \$ ldh, wv( 1, 1+kzs ), ldwv )
814 *
815 * ==== Multiply by U21 ====
816 *
817  CALL zlaset( 'ALL', jlen, kzs, zero, zero, wv, ldwv )
818  CALL ztrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
819  \$ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
820  \$ ldwv )
821 *
822 * ==== Multiply by U11 ====
823 *
824  CALL zgemm( 'N', 'N', jlen, i2, j2, one,
825  \$ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
826  \$ ldwv )
827 *
828 * ==== Copy left of H to right of scratch ====
829 *
830  CALL zlacpy( 'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
831  \$ wv( 1, 1+i2 ), ldwv )
832 *
833 * ==== Multiply by U21 ====
834 *
835  CALL ztrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
836  \$ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
837 *
838 * ==== Multiply by U22 ====
839 *
840  CALL zgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
841  \$ h( jrow, incol+1+j2 ), ldh,
842  \$ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
843  \$ ldwv )
844 *
845 * ==== Copy it back ====
846 *
847  CALL zlacpy( 'ALL', jlen, kdu, wv, ldwv,
848  \$ h( jrow, incol+1 ), ldh )
849  190 CONTINUE
850 *
851 * ==== Multiply Z (also vertical) ====
852 *
853  IF( wantz ) THEN
854  DO 200 jrow = iloz, ihiz, nv
855  jlen = min( nv, ihiz-jrow+1 )
856 *
857 * ==== Copy right of Z to left of scratch (first
858 * . KZS columns get multiplied by zero) ====
859 *
860  CALL zlacpy( 'ALL', jlen, knz,
861  \$ z( jrow, incol+1+j2 ), ldz,
862  \$ wv( 1, 1+kzs ), ldwv )
863 *
864 * ==== Multiply by U12 ====
865 *
866  CALL zlaset( 'ALL', jlen, kzs, zero, zero, wv,
867  \$ ldwv )
868  CALL ztrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
869  \$ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
870  \$ ldwv )
871 *
872 * ==== Multiply by U11 ====
873 *
874  CALL zgemm( 'N', 'N', jlen, i2, j2, one,
875  \$ z( jrow, incol+1 ), ldz, u, ldu, one,
876  \$ wv, ldwv )
877 *
878 * ==== Copy left of Z to right of scratch ====
879 *
880  CALL zlacpy( 'ALL', jlen, j2, z( jrow, incol+1 ),
881  \$ ldz, wv( 1, 1+i2 ), ldwv )
882 *
883 * ==== Multiply by U21 ====
884 *
885  CALL ztrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
886  \$ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
887  \$ ldwv )
888 *
889 * ==== Multiply by U22 ====
890 *
891  CALL zgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
892  \$ z( jrow, incol+1+j2 ), ldz,
893  \$ u( j2+1, i2+1 ), ldu, one,
894  \$ wv( 1, 1+i2 ), ldwv )
895 *
896 * ==== Copy the result back to Z ====
897 *
898  CALL zlacpy( 'ALL', jlen, kdu, wv, ldwv,
899  \$ z( jrow, incol+1 ), ldz )
900  200 CONTINUE
901  END IF
902  END IF
903  END IF
904  210 CONTINUE
905 *
906 * ==== End of ZLAQR5 ====
907 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108