LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ dlaqr5()

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

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

Purpose:
```    DLAQR5, called by DLAQR0, performs a
single small-bulge multi-shift QR sweep.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL WANTT = .true. if the quasi-triangular Schur factor is being computed. WANTT is set to .false. otherwise.``` [in] WANTZ ``` WANTZ is LOGICAL WANTZ = .true. if the orthogonal 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: DLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries. = 1: DLAQR5 accumulates reflections and uses matrix-matrix multiply to update the far-from-diagonal matrix entries. = 2: DLAQR5 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] SR ` SR is DOUBLE PRECISION array, dimension (NSHFTS)` [in,out] SI ``` SI is DOUBLE PRECISION array, dimension (NSHFTS) SR contains the real parts and SI contains the imaginary parts of the NSHFTS shifts of origin that define the multi-shift QR sweep. On output SR and SI may be reordered.``` [in,out] H ``` H is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDZ,IHIZ) If WANTZ = .TRUE., then the QR Sweep orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 259 of file dlaqr5.f.

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