--- rpl/lapack/lapack/dlaqr5.f 2010/12/21 13:53:31 1.8
+++ rpl/lapack/lapack/dlaqr5.f 2023/08/07 08:38:56 1.23
@@ -1,11 +1,272 @@
+*> \brief \b DLAQR5 performs a single small-bulge multi-shift QR sweep.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLAQR5 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
+* SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U,
+* LDU, NV, WV, LDWV, NH, WH, LDWH )
+*
+* .. Scalar Arguments ..
+* INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
+* $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
+* LOGICAL WANTT, WANTZ
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), U( LDU, * ),
+* $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ),
+* $ Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DLAQR5, called by DLAQR0, performs a
+*> single small-bulge multi-shift QR sweep.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] WANTT
+*> \verbatim
+*> WANTT is LOGICAL
+*> WANTT = .true. if the quasi-triangular Schur factor
+*> is being computed. WANTT is set to .false. otherwise.
+*> \endverbatim
+*>
+*> \param[in] WANTZ
+*> \verbatim
+*> WANTZ is LOGICAL
+*> WANTZ = .true. if the orthogonal Schur factor is being
+*> computed. WANTZ is set to .false. otherwise.
+*> \endverbatim
+*>
+*> \param[in] KACC22
+*> \verbatim
+*> 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: Same as KACC22 = 1. This option used to enable exploiting
+*> the 2-by-2 structure during matrix multiplications, but
+*> this is no longer supported.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> N is the order of the Hessenberg matrix H upon which this
+*> subroutine operates.
+*> \endverbatim
+*>
+*> \param[in] KTOP
+*> \verbatim
+*> KTOP is INTEGER
+*> \endverbatim
+*>
+*> \param[in] KBOT
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in] NSHFTS
+*> \verbatim
+*> NSHFTS is INTEGER
+*> NSHFTS gives the number of simultaneous shifts. NSHFTS
+*> must be positive and even.
+*> \endverbatim
+*>
+*> \param[in,out] SR
+*> \verbatim
+*> SR is DOUBLE PRECISION array, dimension (NSHFTS)
+*> \endverbatim
+*>
+*> \param[in,out] SI
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in,out] H
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in] LDH
+*> \verbatim
+*> LDH is INTEGER
+*> LDH is the leading dimension of H just as declared in the
+*> calling procedure. LDH >= MAX(1,N).
+*> \endverbatim
+*>
+*> \param[in] ILOZ
+*> \verbatim
+*> ILOZ is INTEGER
+*> \endverbatim
+*>
+*> \param[in] IHIZ
+*> \verbatim
+*> IHIZ is INTEGER
+*> Specify the rows of Z to which transformations must be
+*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N
+*> \endverbatim
+*>
+*> \param[in,out] Z
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> LDA is the leading dimension of Z just as declared in
+*> the calling procedure. LDZ >= N.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*> V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2)
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*> LDV is INTEGER
+*> LDV is the leading dimension of V as declared in the
+*> calling procedure. LDV >= 3.
+*> \endverbatim
+*>
+*> \param[out] U
+*> \verbatim
+*> U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS)
+*> \endverbatim
+*>
+*> \param[in] LDU
+*> \verbatim
+*> LDU is INTEGER
+*> LDU is the leading dimension of U just as declared in the
+*> in the calling subroutine. LDU >= 2*NSHFTS.
+*> \endverbatim
+*>
+*> \param[in] NV
+*> \verbatim
+*> NV is INTEGER
+*> NV is the number of rows in WV agailable for workspace.
+*> NV >= 1.
+*> \endverbatim
+*>
+*> \param[out] WV
+*> \verbatim
+*> WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS)
+*> \endverbatim
+*>
+*> \param[in] LDWV
+*> \verbatim
+*> LDWV is INTEGER
+*> LDWV is the leading dimension of WV as declared in the
+*> in the calling subroutine. LDWV >= NV.
+*> \endverbatim
+*
+*> \param[in] NH
+*> \verbatim
+*> NH is INTEGER
+*> NH is the number of columns in array WH available for
+*> workspace. NH >= 1.
+*> \endverbatim
+*>
+*> \param[out] WH
+*> \verbatim
+*> WH is DOUBLE PRECISION array, dimension (LDWH,NH)
+*> \endverbatim
+*>
+*> \param[in] LDWH
+*> \verbatim
+*> LDWH is INTEGER
+*> Leading dimension of WH just as declared in the
+*> calling procedure. LDWH >= 2*NSHFTS.
+*> \endverbatim
+*>
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*> \par Contributors:
+* ==================
+*>
+*> Karen Braman and Ralph Byers, Department of Mathematics,
+*> University of Kansas, USA
+*>
+*> Lars Karlsson, Daniel Kressner, and Bruno Lang
+*>
+*> Thijs Steel, Department of Computer science,
+*> KU Leuven, Belgium
+*
+*> \par References:
+* ================
+*>
+*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
+*> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
+*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
+*> 929--947, 2002.
+*>
+*> Lars Karlsson, Daniel Kressner, and Bruno Lang, Optimally packed
+*> chains of bulges in multishift QR algorithms.
+*> ACM Trans. Math. Softw. 40, 2, Article 12 (February 2014).
+*>
+* =====================================================================
SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS,
$ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U,
$ LDU, NV, WV, LDWV, NH, WH, LDWH )
+ IMPLICIT NONE
*
-* -- LAPACK auxiliary routine (version 3.3.0) --
+* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2010
*
* .. Scalar Arguments ..
INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
@@ -18,140 +279,20 @@
$ Z( LDZ, * )
* ..
*
-* This auxiliary subroutine called by DLAQR0 performs a
-* single small-bulge multi-shift QR sweep.
-*
-* WANTT (input) logical scalar
-* WANTT = .true. if the quasi-triangular Schur factor
-* is being computed. WANTT is set to .false. otherwise.
-*
-* WANTZ (input) logical scalar
-* WANTZ = .true. if the orthogonal Schur factor is being
-* computed. WANTZ is set to .false. otherwise.
-*
-* KACC22 (input) 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.
-*
-* N (input) integer scalar
-* N is the order of the Hessenberg matrix H upon which this
-* subroutine operates.
-*
-* KTOP (input) integer scalar
-* KBOT (input) 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.
-*
-* NSHFTS (input) integer scalar
-* NSHFTS gives the number of simultaneous shifts. NSHFTS
-* must be positive and even.
-*
-* SR (input/output) DOUBLE PRECISION array of size (NSHFTS)
-* SI (input/output) DOUBLE PRECISION array of size (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.
-*
-* H (input/output) DOUBLE PRECISION 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.
-*
-* LDH (input) integer scalar
-* LDH is the leading dimension of H just as declared in the
-* calling procedure. LDH.GE.MAX(1,N).
-*
-* ILOZ (input) INTEGER
-* IHIZ (input) INTEGER
-* Specify the rows of Z to which transformations must be
-* applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
-*
-* Z (input/output) DOUBLE PRECISION array of size (LDZ,IHI)
-* If WANTZ = .TRUE., then the QR Sweep orthogonal
-* similarity transformation is accumulated into
-* Z(ILOZ:IHIZ,ILO:IHI) from the right.
-* If WANTZ = .FALSE., then Z is unreferenced.
-*
-* LDZ (input) integer scalar
-* LDA is the leading dimension of Z just as declared in
-* the calling procedure. LDZ.GE.N.
-*
-* V (workspace) DOUBLE PRECISION array of size (LDV,NSHFTS/2)
-*
-* LDV (input) integer scalar
-* LDV is the leading dimension of V as declared in the
-* calling procedure. LDV.GE.3.
-*
-* U (workspace) DOUBLE PRECISION array of size
-* (LDU,3*NSHFTS-3)
-*
-* LDU (input) integer scalar
-* LDU is the leading dimension of U just as declared in the
-* in the calling subroutine. LDU.GE.3*NSHFTS-3.
-*
-* NH (input) integer scalar
-* NH is the number of columns in array WH available for
-* workspace. NH.GE.1.
-*
-* WH (workspace) DOUBLE PRECISION array of size (LDWH,NH)
-*
-* LDWH (input) integer scalar
-* Leading dimension of WH just as declared in the
-* calling procedure. LDWH.GE.3*NSHFTS-3.
-*
-* NV (input) integer scalar
-* NV is the number of rows in WV agailable for workspace.
-* NV.GE.1.
-*
-* WV (workspace) DOUBLE PRECISION array of size
-* (LDWV,3*NSHFTS-3)
-*
-* LDWV (input) integer scalar
-* LDWV is the leading dimension of WV as declared in the
-* in the calling subroutine. LDWV.GE.NV.
-*
-* ================================================================
-* Based on contributions by
-* Karen Braman and Ralph Byers, Department of Mathematics,
-* University of Kansas, USA
-*
-* ================================================================
-* Reference:
-*
-* 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.
-*
-* ================================================================
+* ================================================================
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
* ..
* .. Local Scalars ..
DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM,
- $ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, TST1, TST2,
- $ ULP
- INTEGER I, I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
- $ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
- $ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
+ $ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, T1, T2,
+ $ T3, TST1, TST2, ULP
+ INTEGER I, I2, I4, INCOL, J, JBOT, JCOL, JLEN,
+ $ JROW, JTOP, K, K1, KDU, KMS, KRCOL,
+ $ M, M22, MBOT, MTOP, NBMPS, NDCOL,
$ NS, NU
- LOGICAL ACCUM, BLK22, BMP22
+ LOGICAL ACCUM, BMP22
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
@@ -221,10 +362,6 @@
*
ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
*
-* ==== If so, exploit the 2-by-2 block structure? ====
-*
- BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
-*
* ==== clear trash ====
*
IF( KTOP+2.LE.KBOT )
@@ -236,28 +373,39 @@
*
* ==== KDU = width of slab ====
*
- KDU = 6*NBMPS - 3
+ KDU = 4*NBMPS
*
* ==== Create and chase chains of NBMPS bulges ====
*
- DO 220 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
+ DO 180 INCOL = KTOP - 2*NBMPS + 1, KBOT - 2, 2*NBMPS
+*
+* JTOP = Index from which updates from the right start.
+*
+ IF( ACCUM ) THEN
+ JTOP = MAX( KTOP, INCOL )
+ ELSE IF( WANTT ) THEN
+ JTOP = 1
+ ELSE
+ JTOP = KTOP
+ END IF
+*
NDCOL = INCOL + KDU
IF( ACCUM )
$ CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )
*
* ==== Near-the-diagonal bulge chase. The following loop
* . performs the near-the-diagonal part of a small bulge
-* . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
+* . multi-shift QR sweep. Each 4*NBMPS column diagonal
* . chunk extends from column INCOL to column NDCOL
* . (including both column INCOL and column NDCOL). The
-* . following loop chases a 3*NBMPS column long chain of
-* . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
+* . following loop chases a 2*NBMPS+1 column long chain of
+* . NBMPS bulges 2*NBMPS columns to the right. (INCOL
* . may be less than KTOP and and NDCOL may be greater than
* . KBOT indicating phantom columns from which to chase
* . bulges before they are actually introduced or to which
* . to chase bulges beyond column KBOT.) ====
*
- DO 150 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
+ DO 145 KRCOL = INCOL, MIN( INCOL+2*NBMPS-1, KBOT-2 )
*
* ==== Bulges number MTOP to MBOT are active double implicit
* . shift bulges. There may or may not also be small
@@ -266,17 +414,138 @@
* . down the diagonal to make room. The phantom matrix
* . paradigm described above helps keep track. ====
*
- MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
- MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
+ MTOP = MAX( 1, ( KTOP-KRCOL ) / 2+1 )
+ MBOT = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 2 )
M22 = MBOT + 1
- BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
+ BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+2*( M22-1 ) ).EQ.
$ ( KBOT-2 )
*
* ==== Generate reflections to chase the chain right
* . one column. (The minimum value of K is KTOP-1.) ====
*
- DO 20 M = MTOP, MBOT
- K = KRCOL + 3*( M-1 )
+ IF ( BMP22 ) THEN
+*
+* ==== Special case: 2-by-2 reflection at bottom treated
+* . separately ====
+*
+ K = KRCOL + 2*( M22-1 )
+ IF( K.EQ.KTOP-1 ) THEN
+ CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ),
+ $ SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ),
+ $ V( 1, M22 ) )
+ BETA = V( 1, M22 )
+ CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
+ ELSE
+ BETA = H( K+1, K )
+ V( 2, M22 ) = H( K+2, K )
+ CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
+ H( K+1, K ) = BETA
+ H( K+2, K ) = ZERO
+ END IF
+
+*
+* ==== Perform update from right within
+* . computational window. ====
+*
+ T1 = V( 1, M22 )
+ T2 = T1*V( 2, M22 )
+ DO 30 J = JTOP, MIN( KBOT, K+3 )
+ REFSUM = H( J, K+1 ) + V( 2, M22 )*H( J, K+2 )
+ H( J, K+1 ) = H( J, K+1 ) - REFSUM*T1
+ H( J, K+2 ) = H( J, K+2 ) - REFSUM*T2
+ 30 CONTINUE
+*
+* ==== Perform update from left within
+* . computational window. ====
+*
+ IF( ACCUM ) THEN
+ JBOT = MIN( NDCOL, KBOT )
+ ELSE IF( WANTT ) THEN
+ JBOT = N
+ ELSE
+ JBOT = KBOT
+ END IF
+ T1 = V( 1, M22 )
+ T2 = T1*V( 2, M22 )
+ DO 40 J = K+1, JBOT
+ REFSUM = H( K+1, J ) + V( 2, M22 )*H( K+2, J )
+ H( K+1, J ) = H( K+1, J ) - REFSUM*T1
+ H( K+2, J ) = H( K+2, J ) - REFSUM*T2
+ 40 CONTINUE
+*
+* ==== The following convergence test requires that
+* . the tradition small-compared-to-nearby-diagonals
+* . criterion and the Ahues & Tisseur (LAWN 122, 1997)
+* . criteria both be satisfied. The latter improves
+* . accuracy in some examples. Falling back on an
+* . alternate convergence criterion when TST1 or TST2
+* . is zero (as done here) is traditional but probably
+* . unnecessary. ====
+*
+ IF( K.GE.KTOP ) THEN
+ IF( H( K+1, K ).NE.ZERO ) THEN
+ TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) )
+ IF( TST1.EQ.ZERO ) THEN
+ IF( K.GE.KTOP+1 )
+ $ TST1 = TST1 + ABS( H( K, K-1 ) )
+ IF( K.GE.KTOP+2 )
+ $ TST1 = TST1 + ABS( H( K, K-2 ) )
+ IF( K.GE.KTOP+3 )
+ $ TST1 = TST1 + ABS( H( K, K-3 ) )
+ IF( K.LE.KBOT-2 )
+ $ TST1 = TST1 + ABS( H( K+2, K+1 ) )
+ IF( K.LE.KBOT-3 )
+ $ TST1 = TST1 + ABS( H( K+3, K+1 ) )
+ IF( K.LE.KBOT-4 )
+ $ TST1 = TST1 + ABS( H( K+4, K+1 ) )
+ END IF
+ IF( ABS( H( K+1, K ) )
+ $ .LE.MAX( SMLNUM, ULP*TST1 ) ) THEN
+ H12 = MAX( ABS( H( K+1, K ) ),
+ $ ABS( H( K, K+1 ) ) )
+ H21 = MIN( ABS( H( K+1, K ) ),
+ $ ABS( H( K, K+1 ) ) )
+ H11 = MAX( ABS( H( K+1, K+1 ) ),
+ $ ABS( H( K, K )-H( K+1, K+1 ) ) )
+ H22 = MIN( ABS( H( K+1, K+1 ) ),
+ $ ABS( H( K, K )-H( K+1, K+1 ) ) )
+ SCL = H11 + H12
+ TST2 = H22*( H11 / SCL )
+*
+ IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE.
+ $ MAX( SMLNUM, ULP*TST2 ) ) THEN
+ H( K+1, K ) = ZERO
+ END IF
+ END IF
+ END IF
+ END IF
+*
+* ==== Accumulate orthogonal transformations. ====
+*
+ IF( ACCUM ) THEN
+ KMS = K - INCOL
+ T1 = V( 1, M22 )
+ T2 = T1*V( 2, M22 )
+ DO 50 J = MAX( 1, KTOP-INCOL ), KDU
+ REFSUM = U( J, KMS+1 ) + V( 2, M22 )*U( J, KMS+2 )
+ U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM*T1
+ U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*T2
+ 50 CONTINUE
+ ELSE IF( WANTZ ) THEN
+ T1 = V( 1, M22 )
+ T2 = T1*V( 2, M22 )
+ DO 60 J = ILOZ, IHIZ
+ REFSUM = Z( J, K+1 )+V( 2, M22 )*Z( J, K+2 )
+ Z( J, K+1 ) = Z( J, K+1 ) - REFSUM*T1
+ Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*T2
+ 60 CONTINUE
+ END IF
+ END IF
+*
+* ==== Normal case: Chain of 3-by-3 reflections ====
+*
+ DO 80 M = MBOT, MTOP, -1
+ K = KRCOL + 2*( M-1 )
IF( K.EQ.KTOP-1 ) THEN
CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ),
$ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ),
@@ -284,7 +553,20 @@
ALPHA = V( 1, M )
CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
ELSE
- BETA = H( K+1, K )
+*
+* ==== Perform delayed transformation of row below
+* . Mth bulge. Exploit fact that first two elements
+* . of row are actually zero. ====
+*
+ REFSUM = V( 1, M )*V( 3, M )*H( K+3, K+2 )
+ H( K+3, K ) = -REFSUM
+ H( K+3, K+1 ) = -REFSUM*V( 2, M )
+ H( K+3, K+2 ) = H( K+3, K+2 ) - REFSUM*V( 3, M )
+*
+* ==== Calculate reflection to move
+* . Mth bulge one step. ====
+*
+ BETA = H( K+1, K )
V( 2, M ) = H( K+2, K )
V( 3, M ) = H( K+3, K )
CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )
@@ -332,7 +614,7 @@
H( K+3, K ) = ZERO
ELSE
*
-* ==== Stating a new bulge here would
+* ==== Starting a new bulge here would
* . create only negligible fill.
* . Replace the old reflector with
* . the new one. ====
@@ -346,154 +628,32 @@
END IF
END IF
END IF
- 20 CONTINUE
-*
-* ==== Generate a 2-by-2 reflection, if needed. ====
-*
- K = KRCOL + 3*( M22-1 )
- IF( BMP22 ) THEN
- IF( K.EQ.KTOP-1 ) THEN
- CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ),
- $ SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ),
- $ V( 1, M22 ) )
- BETA = V( 1, M22 )
- CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
- ELSE
- BETA = H( K+1, K )
- V( 2, M22 ) = H( K+2, K )
- CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
- H( K+1, K ) = BETA
- H( K+2, K ) = ZERO
- END IF
- END IF
-*
-* ==== Multiply H by reflections from the left ====
-*
- IF( ACCUM ) THEN
- JBOT = MIN( NDCOL, KBOT )
- ELSE IF( WANTT ) THEN
- JBOT = N
- ELSE
- JBOT = KBOT
- END IF
- DO 40 J = MAX( KTOP, KRCOL ), JBOT
- MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
- DO 30 M = MTOP, MEND
- K = KRCOL + 3*( M-1 )
- REFSUM = V( 1, M )*( H( K+1, J )+V( 2, M )*
- $ H( K+2, J )+V( 3, M )*H( K+3, J ) )
- H( K+1, J ) = H( K+1, J ) - REFSUM
- H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
- H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
- 30 CONTINUE
- 40 CONTINUE
- IF( BMP22 ) THEN
- K = KRCOL + 3*( M22-1 )
- DO 50 J = MAX( K+1, KTOP ), JBOT
- REFSUM = V( 1, M22 )*( H( K+1, J )+V( 2, M22 )*
- $ H( K+2, J ) )
- H( K+1, J ) = H( K+1, J ) - REFSUM
- H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
- 50 CONTINUE
- END IF
-*
-* ==== Multiply H by reflections from the right.
-* . Delay filling in the last row until the
-* . vigilant deflation check is complete. ====
-*
- IF( ACCUM ) THEN
- JTOP = MAX( KTOP, INCOL )
- ELSE IF( WANTT ) THEN
- JTOP = 1
- ELSE
- JTOP = KTOP
- END IF
- DO 90 M = MTOP, MBOT
- IF( V( 1, M ).NE.ZERO ) THEN
- K = KRCOL + 3*( M-1 )
- DO 60 J = JTOP, MIN( KBOT, K+3 )
- REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
- $ H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
- H( J, K+1 ) = H( J, K+1 ) - REFSUM
- H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M )
- H( J, K+3 ) = H( J, K+3 ) - REFSUM*V( 3, M )
- 60 CONTINUE
-*
- IF( ACCUM ) THEN
-*
-* ==== Accumulate U. (If necessary, update Z later
-* . with with an efficient matrix-matrix
-* . multiply.) ====
-*
- KMS = K - INCOL
- DO 70 J = MAX( 1, KTOP-INCOL ), KDU
- REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
- $ U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
- U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
- U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M )
- U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*V( 3, M )
- 70 CONTINUE
- ELSE IF( WANTZ ) THEN
-*
-* ==== U is not accumulated, so update Z
-* . now by multiplying by reflections
-* . from the right. ====
-*
- DO 80 J = ILOZ, IHIZ
- REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
- $ Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
- Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
- Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M )
- Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M )
- 80 CONTINUE
- END IF
- END IF
- 90 CONTINUE
-*
-* ==== Special case: 2-by-2 reflection (if needed) ====
*
- K = KRCOL + 3*( M22-1 )
- IF( BMP22 ) THEN
- IF ( V( 1, M22 ).NE.ZERO ) THEN
- DO 100 J = JTOP, MIN( KBOT, K+3 )
- REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
- $ H( J, K+2 ) )
- H( J, K+1 ) = H( J, K+1 ) - REFSUM
- H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 )
- 100 CONTINUE
-*
- IF( ACCUM ) THEN
- KMS = K - INCOL
- DO 110 J = MAX( 1, KTOP-INCOL ), KDU
- REFSUM = V( 1, M22 )*( U( J, KMS+1 )+
- $ V( 2, M22 )*U( J, KMS+2 ) )
- U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
- U( J, KMS+2 ) = U( J, KMS+2 ) -
- $ REFSUM*V( 2, M22 )
- 110 CONTINUE
- ELSE IF( WANTZ ) THEN
- DO 120 J = ILOZ, IHIZ
- REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
- $ Z( J, K+2 ) )
- Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
- Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 )
- 120 CONTINUE
- END IF
- END IF
- END IF
-*
-* ==== Vigilant deflation check ====
-*
- MSTART = MTOP
- IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
- $ MSTART = MSTART + 1
- MEND = MBOT
- IF( BMP22 )
- $ MEND = MEND + 1
- IF( KRCOL.EQ.KBOT-2 )
- $ MEND = MEND + 1
- DO 130 M = MSTART, MEND
- K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
+* ==== Apply reflection from the right and
+* . the first column of update from the left.
+* . These updates are required for the vigilant
+* . deflation check. We still delay most of the
+* . updates from the left for efficiency. ====
+*
+ T1 = V( 1, M )
+ T2 = T1*V( 2, M )
+ T3 = T1*V( 3, M )
+ DO 70 J = JTOP, MIN( KBOT, K+3 )
+ REFSUM = H( J, K+1 ) + V( 2, M )*H( J, K+2 )
+ $ + V( 3, M )*H( J, K+3 )
+ H( J, K+1 ) = H( J, K+1 ) - REFSUM*T1
+ H( J, K+2 ) = H( J, K+2 ) - REFSUM*T2
+ H( J, K+3 ) = H( J, K+3 ) - REFSUM*T3
+ 70 CONTINUE
+*
+* ==== Perform update from left for subsequent
+* . column. ====
+*
+ REFSUM = H( K+1, K+1 ) + V( 2, M )*H( K+2, K+1 )
+ $ + V( 3, M )*H( K+3, K+1 )
+ H( K+1, K+1 ) = H( K+1, K+1 ) - REFSUM*T1
+ H( K+2, K+1 ) = H( K+2, K+1 ) - REFSUM*T2
+ H( K+3, K+1 ) = H( K+3, K+1 ) - REFSUM*T3
*
* ==== The following convergence test requires that
* . the tradition small-compared-to-nearby-diagonals
@@ -504,6 +664,8 @@
* . is zero (as done here) is traditional but probably
* . unnecessary. ====
*
+ IF( K.LT.KTOP)
+ $ CYCLE
IF( H( K+1, K ).NE.ZERO ) THEN
TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) )
IF( TST1.EQ.ZERO ) THEN
@@ -532,25 +694,86 @@
TST2 = H22*( H11 / SCL )
*
IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE.
- $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
+ $ MAX( SMLNUM, ULP*TST2 ) ) THEN
+ H( K+1, K ) = ZERO
+ END IF
END IF
END IF
- 130 CONTINUE
+ 80 CONTINUE
*
-* ==== Fill in the last row of each bulge. ====
+* ==== Multiply H by reflections from the left ====
*
- MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
- DO 140 M = MTOP, MEND
- K = KRCOL + 3*( M-1 )
- REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
- H( K+4, K+1 ) = -REFSUM
- H( K+4, K+2 ) = -REFSUM*V( 2, M )
- H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M )
- 140 CONTINUE
+ IF( ACCUM ) THEN
+ JBOT = MIN( NDCOL, KBOT )
+ ELSE IF( WANTT ) THEN
+ JBOT = N
+ ELSE
+ JBOT = KBOT
+ END IF
+*
+ DO 100 M = MBOT, MTOP, -1
+ K = KRCOL + 2*( M-1 )
+ T1 = V( 1, M )
+ T2 = T1*V( 2, M )
+ T3 = T1*V( 3, M )
+ DO 90 J = MAX( KTOP, KRCOL + 2*M ), JBOT
+ REFSUM = H( K+1, J ) + V( 2, M )*H( K+2, J )
+ $ + V( 3, M )*H( K+3, J )
+ H( K+1, J ) = H( K+1, J ) - REFSUM*T1
+ H( K+2, J ) = H( K+2, J ) - REFSUM*T2
+ H( K+3, J ) = H( K+3, J ) - REFSUM*T3
+ 90 CONTINUE
+ 100 CONTINUE
+*
+* ==== Accumulate orthogonal transformations. ====
+*
+ IF( ACCUM ) THEN
+*
+* ==== Accumulate U. (If needed, update Z later
+* . with an efficient matrix-matrix
+* . multiply.) ====
+*
+ DO 120 M = MBOT, MTOP, -1
+ K = KRCOL + 2*( M-1 )
+ KMS = K - INCOL
+ I2 = MAX( 1, KTOP-INCOL )
+ I2 = MAX( I2, KMS-(KRCOL-INCOL)+1 )
+ I4 = MIN( KDU, KRCOL + 2*( MBOT-1 ) - INCOL + 5 )
+ T1 = V( 1, M )
+ T2 = T1*V( 2, M )
+ T3 = T1*V( 3, M )
+ DO 110 J = I2, I4
+ REFSUM = U( J, KMS+1 ) + V( 2, M )*U( J, KMS+2 )
+ $ + V( 3, M )*U( J, KMS+3 )
+ U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM*T1
+ U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*T2
+ U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*T3
+ 110 CONTINUE
+ 120 CONTINUE
+ ELSE IF( WANTZ ) THEN
+*
+* ==== U is not accumulated, so update Z
+* . now by multiplying by reflections
+* . from the right. ====
+*
+ DO 140 M = MBOT, MTOP, -1
+ K = KRCOL + 2*( M-1 )
+ T1 = V( 1, M )
+ T2 = T1*V( 2, M )
+ T3 = T1*V( 3, M )
+ DO 130 J = ILOZ, IHIZ
+ REFSUM = Z( J, K+1 ) + V( 2, M )*Z( J, K+2 )
+ $ + V( 3, M )*Z( J, K+3 )
+ Z( J, K+1 ) = Z( J, K+1 ) - REFSUM*T1
+ Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*T2
+ Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*T3
+ 130 CONTINUE
+ 140 CONTINUE
+ END IF
*
* ==== End of near-the-diagonal bulge chase. ====
*
- 150 CONTINUE
+ 145 CONTINUE
*
* ==== Use U (if accumulated) to update far-from-diagonal
* . entries in H. If required, use U to update Z as
@@ -564,220 +787,45 @@
JTOP = KTOP
JBOT = KBOT
END IF
- IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
- $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
+ K1 = MAX( 1, KTOP-INCOL )
+ NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
+*
+* ==== Horizontal Multiply ====
*
-* ==== Updates not exploiting the 2-by-2 block
-* . structure of U. K1 and NU keep track of
-* . the location and size of U in the special
-* . cases of introducing bulges and chasing
-* . bulges off the bottom. In these special
-* . cases and in case the number of shifts
-* . is NS = 2, there is no 2-by-2 block
-* . structure to exploit. ====
-*
- K1 = MAX( 1, KTOP-INCOL )
- NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
-*
-* ==== Horizontal Multiply ====
-*
- DO 160 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
- JLEN = MIN( NH, JBOT-JCOL+1 )
- CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
+ DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
+ JLEN = MIN( NH, JBOT-JCOL+1 )
+ CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
$ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
$ LDWH )
- CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH,
+ CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH,
$ H( INCOL+K1, JCOL ), LDH )
- 160 CONTINUE
+ 150 CONTINUE
*
-* ==== Vertical multiply ====
+* ==== Vertical multiply ====
*
- DO 170 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
- JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
+ DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
+ JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
+ CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE,
+ $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
+ $ LDU, ZERO, WV, LDWV )
+ CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV,
+ $ H( JROW, INCOL+K1 ), LDH )
+ 160 CONTINUE
+*
+* ==== Z multiply (also vertical) ====
+*
+ IF( WANTZ ) THEN
+ DO 170 JROW = ILOZ, IHIZ, NV
+ JLEN = MIN( NV, IHIZ-JROW+1 )
CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE,
- $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
+ $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
$ LDU, ZERO, WV, LDWV )
CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV,
- $ H( JROW, INCOL+K1 ), LDH )
+ $ Z( JROW, INCOL+K1 ), LDZ )
170 CONTINUE
-*
-* ==== Z multiply (also vertical) ====
-*
- IF( WANTZ ) THEN
- DO 180 JROW = ILOZ, IHIZ, NV
- JLEN = MIN( NV, IHIZ-JROW+1 )
- CALL DGEMM( 'N', 'N', JLEN, NU, NU, ONE,
- $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
- $ LDU, ZERO, WV, LDWV )
- CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV,
- $ Z( JROW, INCOL+K1 ), LDZ )
- 180 CONTINUE
- END IF
- ELSE
-*
-* ==== Updates exploiting U's 2-by-2 block structure.
-* . (I2, I4, J2, J4 are the last rows and columns
-* . of the blocks.) ====
-*
- I2 = ( KDU+1 ) / 2
- I4 = KDU
- J2 = I4 - I2
- J4 = KDU
-*
-* ==== KZS and KNZ deal with the band of zeros
-* . along the diagonal of one of the triangular
-* . blocks. ====
-*
- KZS = ( J4-J2 ) - ( NS+1 )
- KNZ = NS + 1
-*
-* ==== Horizontal multiply ====
-*
- DO 190 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
- JLEN = MIN( NH, JBOT-JCOL+1 )
-*
-* ==== Copy bottom of H to top+KZS of scratch ====
-* (The first KZS rows get multiplied by zero.) ====
-*
- CALL DLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
- $ LDH, WH( KZS+1, 1 ), LDWH )
-*
-* ==== Multiply by U21' ====
-*
- CALL DLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
- CALL DTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
- $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
- $ LDWH )
-*
-* ==== Multiply top of H by U11' ====
-*
- CALL DGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
- $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
-*
-* ==== Copy top of H to bottom of WH ====
-*
- CALL DLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
- $ WH( I2+1, 1 ), LDWH )
-*
-* ==== Multiply by U21' ====
-*
- CALL DTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
- $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
-*
-* ==== Multiply by U22 ====
-*
- CALL DGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
- $ U( J2+1, I2+1 ), LDU,
- $ H( INCOL+1+J2, JCOL ), LDH, ONE,
- $ WH( I2+1, 1 ), LDWH )
-*
-* ==== Copy it back ====
-*
- CALL DLACPY( 'ALL', KDU, JLEN, WH, LDWH,
- $ H( INCOL+1, JCOL ), LDH )
- 190 CONTINUE
-*
-* ==== Vertical multiply ====
-*
- DO 200 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
- JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
-*
-* ==== Copy right of H to scratch (the first KZS
-* . columns get multiplied by zero) ====
-*
- CALL DLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
- $ LDH, WV( 1, 1+KZS ), LDWV )
-*
-* ==== Multiply by U21 ====
-*
- CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
- CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
- $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
- $ LDWV )
-*
-* ==== Multiply by U11 ====
-*
- CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE,
- $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
- $ LDWV )
-*
-* ==== Copy left of H to right of scratch ====
-*
- CALL DLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
- $ WV( 1, 1+I2 ), LDWV )
-*
-* ==== Multiply by U21 ====
-*
- CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
- $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
-*
-* ==== Multiply by U22 ====
-*
- CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
- $ H( JROW, INCOL+1+J2 ), LDH,
- $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
- $ LDWV )
-*
-* ==== Copy it back ====
-*
- CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV,
- $ H( JROW, INCOL+1 ), LDH )
- 200 CONTINUE
-*
-* ==== Multiply Z (also vertical) ====
-*
- IF( WANTZ ) THEN
- DO 210 JROW = ILOZ, IHIZ, NV
- JLEN = MIN( NV, IHIZ-JROW+1 )
-*
-* ==== Copy right of Z to left of scratch (first
-* . KZS columns get multiplied by zero) ====
-*
- CALL DLACPY( 'ALL', JLEN, KNZ,
- $ Z( JROW, INCOL+1+J2 ), LDZ,
- $ WV( 1, 1+KZS ), LDWV )
-*
-* ==== Multiply by U12 ====
-*
- CALL DLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
- $ LDWV )
- CALL DTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
- $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
- $ LDWV )
-*
-* ==== Multiply by U11 ====
-*
- CALL DGEMM( 'N', 'N', JLEN, I2, J2, ONE,
- $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
- $ WV, LDWV )
-*
-* ==== Copy left of Z to right of scratch ====
-*
- CALL DLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
- $ LDZ, WV( 1, 1+I2 ), LDWV )
-*
-* ==== Multiply by U21 ====
-*
- CALL DTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
- $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
- $ LDWV )
-*
-* ==== Multiply by U22 ====
-*
- CALL DGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
- $ Z( JROW, INCOL+1+J2 ), LDZ,
- $ U( J2+1, I2+1 ), LDU, ONE,
- $ WV( 1, 1+I2 ), LDWV )
-*
-* ==== Copy the result back to Z ====
-*
- CALL DLACPY( 'ALL', JLEN, KDU, WV, LDWV,
- $ Z( JROW, INCOL+1 ), LDZ )
- 210 CONTINUE
- END IF
END IF
END IF
- 220 CONTINUE
+ 180 CONTINUE
*
* ==== End of DLAQR5 ====
*