version 1.12, 2012/08/22 09:48:18
|
version 1.23, 2023/08/07 08:38:56
|
Line 1
|
Line 1
|
*> \brief \b DLAQR5 |
*> \brief \b DLAQR5 performs a single small-bulge multi-shift QR sweep. |
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DLAQR5 + dependencies |
*> Download DLAQR5 + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr5.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
Line 21
|
Line 21
|
* SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, |
* SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, |
* SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, |
* SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, |
* LDU, NV, WV, LDWV, NH, WH, LDWH ) |
* LDU, NV, WV, LDWV, NH, WH, LDWH ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, |
* INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, |
* $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV |
* $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV |
Line 32
|
Line 32
|
* $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), |
* $ V( LDV, * ), WH( LDWH, * ), WV( LDWV, * ), |
* $ Z( LDZ, * ) |
* $ Z( LDZ, * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 48
|
Line 48
|
* |
* |
*> \param[in] WANTT |
*> \param[in] WANTT |
*> \verbatim |
*> \verbatim |
*> WANTT is logical scalar |
*> WANTT is LOGICAL |
*> WANTT = .true. if the quasi-triangular Schur factor |
*> WANTT = .true. if the quasi-triangular Schur factor |
*> is being computed. WANTT is set to .false. otherwise. |
*> is being computed. WANTT is set to .false. otherwise. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] WANTZ |
*> \param[in] WANTZ |
*> \verbatim |
*> \verbatim |
*> WANTZ is logical scalar |
*> WANTZ is LOGICAL |
*> WANTZ = .true. if the orthogonal Schur factor is being |
*> WANTZ = .true. if the orthogonal Schur factor is being |
*> computed. WANTZ is set to .false. otherwise. |
*> computed. WANTZ is set to .false. otherwise. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] KACC22 |
*> \param[in] KACC22 |
*> \verbatim |
*> \verbatim |
*> KACC22 is integer with value 0, 1, or 2. |
*> KACC22 is INTEGER with value 0, 1, or 2. |
*> Specifies the computation mode of far-from-diagonal |
*> Specifies the computation mode of far-from-diagonal |
*> orthogonal updates. |
*> orthogonal updates. |
*> = 0: DLAQR5 does not accumulate reflections and does not |
*> = 0: DLAQR5 does not accumulate reflections and does not |
Line 70
|
Line 70
|
*> matrix entries. |
*> matrix entries. |
*> = 1: DLAQR5 accumulates reflections and uses matrix-matrix |
*> = 1: DLAQR5 accumulates reflections and uses matrix-matrix |
*> multiply to update the far-from-diagonal matrix entries. |
*> multiply to update the far-from-diagonal matrix entries. |
*> = 2: DLAQR5 accumulates reflections, uses matrix-matrix |
*> = 2: Same as KACC22 = 1. This option used to enable exploiting |
*> multiply to update the far-from-diagonal matrix entries, |
*> the 2-by-2 structure during matrix multiplications, but |
*> and takes advantage of 2-by-2 block structure during |
*> this is no longer supported. |
*> matrix multiplies. |
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] N |
*> \param[in] N |
*> \verbatim |
*> \verbatim |
*> N is integer scalar |
*> N is INTEGER |
*> N is the order of the Hessenberg matrix H upon which this |
*> N is the order of the Hessenberg matrix H upon which this |
*> subroutine operates. |
*> subroutine operates. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] KTOP |
*> \param[in] KTOP |
*> \verbatim |
*> \verbatim |
*> KTOP is integer scalar |
*> KTOP is INTEGER |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] KBOT |
*> \param[in] KBOT |
*> \verbatim |
*> \verbatim |
*> KBOT is integer scalar |
*> KBOT is INTEGER |
*> These are the first and last rows and columns of an |
*> These are the first and last rows and columns of an |
*> isolated diagonal block upon which the QR sweep is to be |
*> isolated diagonal block upon which the QR sweep is to be |
*> applied. It is assumed without a check that |
*> applied. It is assumed without a check that |
Line 101
|
Line 100
|
*> |
*> |
*> \param[in] NSHFTS |
*> \param[in] NSHFTS |
*> \verbatim |
*> \verbatim |
*> NSHFTS is integer scalar |
*> NSHFTS is INTEGER |
*> NSHFTS gives the number of simultaneous shifts. NSHFTS |
*> NSHFTS gives the number of simultaneous shifts. NSHFTS |
*> must be positive and even. |
*> must be positive and even. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] SR |
*> \param[in,out] SR |
*> \verbatim |
*> \verbatim |
*> SR is DOUBLE PRECISION array of size (NSHFTS) |
*> SR is DOUBLE PRECISION array, dimension (NSHFTS) |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] SI |
*> \param[in,out] SI |
*> \verbatim |
*> \verbatim |
*> SI is DOUBLE PRECISION array of size (NSHFTS) |
*> SI is DOUBLE PRECISION array, dimension (NSHFTS) |
*> SR contains the real parts and SI contains the imaginary |
*> SR contains the real parts and SI contains the imaginary |
*> parts of the NSHFTS shifts of origin that define the |
*> parts of the NSHFTS shifts of origin that define the |
*> multi-shift QR sweep. On output SR and SI may be |
*> multi-shift QR sweep. On output SR and SI may be |
Line 122
|
Line 121
|
*> |
*> |
*> \param[in,out] H |
*> \param[in,out] H |
*> \verbatim |
*> \verbatim |
*> H is DOUBLE PRECISION array of size (LDH,N) |
*> H is DOUBLE PRECISION array, dimension (LDH,N) |
*> On input H contains a Hessenberg matrix. On output a |
*> On input H contains a Hessenberg matrix. On output a |
*> multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied |
*> multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied |
*> to the isolated diagonal block in rows and columns KTOP |
*> to the isolated diagonal block in rows and columns KTOP |
Line 131
|
Line 130
|
*> |
*> |
*> \param[in] LDH |
*> \param[in] LDH |
*> \verbatim |
*> \verbatim |
*> LDH is integer scalar |
*> LDH is INTEGER |
*> LDH is the leading dimension of H just as declared in the |
*> LDH is the leading dimension of H just as declared in the |
*> calling procedure. LDH.GE.MAX(1,N). |
*> calling procedure. LDH >= MAX(1,N). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] ILOZ |
*> \param[in] ILOZ |
Line 145
|
Line 144
|
*> \verbatim |
*> \verbatim |
*> IHIZ is INTEGER |
*> IHIZ is INTEGER |
*> Specify the rows of Z to which transformations must be |
*> Specify the rows of Z to which transformations must be |
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N |
*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] Z |
*> \param[in,out] Z |
*> \verbatim |
*> \verbatim |
*> Z is DOUBLE PRECISION array of size (LDZ,IHI) |
*> Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ) |
*> If WANTZ = .TRUE., then the QR Sweep orthogonal |
*> If WANTZ = .TRUE., then the QR Sweep orthogonal |
*> similarity transformation is accumulated into |
*> similarity transformation is accumulated into |
*> Z(ILOZ:IHIZ,ILO:IHI) from the right. |
*> Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. |
*> If WANTZ = .FALSE., then Z is unreferenced. |
*> If WANTZ = .FALSE., then Z is unreferenced. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDZ |
*> \param[in] LDZ |
*> \verbatim |
*> \verbatim |
*> LDZ is integer scalar |
*> LDZ is INTEGER |
*> LDA is the leading dimension of Z just as declared in |
*> LDA is the leading dimension of Z just as declared in |
*> the calling procedure. LDZ.GE.N. |
*> the calling procedure. LDZ >= N. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] V |
*> \param[out] V |
*> \verbatim |
*> \verbatim |
*> V is DOUBLE PRECISION array of size (LDV,NSHFTS/2) |
*> V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2) |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDV |
*> \param[in] LDV |
*> \verbatim |
*> \verbatim |
*> LDV is integer scalar |
*> LDV is INTEGER |
*> LDV is the leading dimension of V as declared in the |
*> LDV is the leading dimension of V as declared in the |
*> calling procedure. LDV.GE.3. |
*> calling procedure. LDV >= 3. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] U |
*> \param[out] U |
*> \verbatim |
*> \verbatim |
*> U is DOUBLE PRECISION array of size |
*> U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS) |
*> (LDU,3*NSHFTS-3) |
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDU |
*> \param[in] LDU |
*> \verbatim |
*> \verbatim |
*> LDU is integer scalar |
*> LDU is INTEGER |
*> LDU is the leading dimension of U just as declared in the |
*> LDU is the leading dimension of U just as declared in the |
*> in the calling subroutine. LDU.GE.3*NSHFTS-3. |
*> in the calling subroutine. LDU >= 2*NSHFTS. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] NH |
*> \param[in] NV |
*> \verbatim |
*> \verbatim |
*> NH is integer scalar |
*> NV is INTEGER |
*> NH is the number of columns in array WH available for |
*> NV is the number of rows in WV agailable for workspace. |
*> workspace. NH.GE.1. |
*> NV >= 1. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] WH |
*> \param[out] WV |
*> \verbatim |
*> \verbatim |
*> WH is DOUBLE PRECISION array of size (LDWH,NH) |
*> WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS) |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDWH |
*> \param[in] LDWV |
*> \verbatim |
*> \verbatim |
*> LDWH is integer scalar |
*> LDWV is INTEGER |
*> Leading dimension of WH just as declared in the |
*> LDWV is the leading dimension of WV as declared in the |
*> calling procedure. LDWH.GE.3*NSHFTS-3. |
*> in the calling subroutine. LDWV >= NV. |
*> \endverbatim |
*> \endverbatim |
*> |
* |
*> \param[in] NV |
*> \param[in] NH |
*> \verbatim |
*> \verbatim |
*> NV is integer scalar |
*> NH is INTEGER |
*> NV is the number of rows in WV agailable for workspace. |
*> NH is the number of columns in array WH available for |
*> NV.GE.1. |
*> workspace. NH >= 1. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] WV |
*> \param[out] WH |
*> \verbatim |
*> \verbatim |
*> WV is DOUBLE PRECISION array of size |
*> WH is DOUBLE PRECISION array, dimension (LDWH,NH) |
*> (LDWV,3*NSHFTS-3) |
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDWV |
*> \param[in] LDWH |
*> \verbatim |
*> \verbatim |
*> LDWV is integer scalar |
*> LDWH is INTEGER |
*> LDWV is the leading dimension of WV as declared in the |
*> Leading dimension of WH just as declared in the |
*> in the calling subroutine. LDWV.GE.NV. |
*> calling procedure. LDWH >= 2*NSHFTS. |
*> \endverbatim |
*> \endverbatim |
* |
*> |
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
|
*> \date November 2011 |
|
* |
* |
*> \ingroup doubleOTHERauxiliary |
*> \ingroup doubleOTHERauxiliary |
* |
* |
Line 245
|
Line 240
|
*> |
*> |
*> Karen Braman and Ralph Byers, Department of Mathematics, |
*> Karen Braman and Ralph Byers, Department of Mathematics, |
*> University of Kansas, USA |
*> University of Kansas, USA |
|
*> |
|
*> Lars Karlsson, Daniel Kressner, and Bruno Lang |
|
*> |
|
*> Thijs Steel, Department of Computer science, |
|
*> KU Leuven, Belgium |
* |
* |
*> \par References: |
*> \par References: |
* ================ |
* ================ |
Line 254
|
Line 254
|
*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages |
*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages |
*> 929--947, 2002. |
*> 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, |
SUBROUTINE DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, |
$ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, |
$ SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, |
$ LDU, NV, WV, LDWV, NH, WH, LDWH ) |
$ LDU, NV, WV, LDWV, NH, WH, LDWH ) |
|
IMPLICIT NONE |
* |
* |
* -- LAPACK auxiliary routine (version 3.4.0) -- |
* -- LAPACK auxiliary routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2011 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, |
INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV, |
Line 282
|
Line 286
|
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM, |
DOUBLE PRECISION ALPHA, BETA, H11, H12, H21, H22, REFSUM, |
$ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, TST1, TST2, |
$ SAFMAX, SAFMIN, SCL, SMLNUM, SWAP, T1, T2, |
$ ULP |
$ T3, TST1, TST2, ULP |
INTEGER I, I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN, |
INTEGER I, I2, I4, INCOL, J, JBOT, JCOL, JLEN, |
$ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS, |
$ JROW, JTOP, K, K1, KDU, KMS, KRCOL, |
$ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL, |
$ M, M22, MBOT, MTOP, NBMPS, NDCOL, |
$ NS, NU |
$ NS, NU |
LOGICAL ACCUM, BLK22, BMP22 |
LOGICAL ACCUM, BMP22 |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
DOUBLE PRECISION DLAMCH |
DOUBLE PRECISION DLAMCH |
Line 358
|
Line 362
|
* |
* |
ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 ) |
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 ==== |
* ==== clear trash ==== |
* |
* |
IF( KTOP+2.LE.KBOT ) |
IF( KTOP+2.LE.KBOT ) |
Line 373
|
Line 373
|
* |
* |
* ==== KDU = width of slab ==== |
* ==== KDU = width of slab ==== |
* |
* |
KDU = 6*NBMPS - 3 |
KDU = 4*NBMPS |
* |
* |
* ==== Create and chase chains of NBMPS bulges ==== |
* ==== 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 |
NDCOL = INCOL + KDU |
IF( ACCUM ) |
IF( ACCUM ) |
$ CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU ) |
$ CALL DLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU ) |
* |
* |
* ==== Near-the-diagonal bulge chase. The following loop |
* ==== Near-the-diagonal bulge chase. The following loop |
* . performs the near-the-diagonal part of a small bulge |
* . 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 |
* . chunk extends from column INCOL to column NDCOL |
* . (including both column INCOL and column NDCOL). The |
* . (including both column INCOL and column NDCOL). The |
* . following loop chases a 3*NBMPS column long chain of |
* . following loop chases a 2*NBMPS+1 column long chain of |
* . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL |
* . NBMPS bulges 2*NBMPS columns to the right. (INCOL |
* . may be less than KTOP and and NDCOL may be greater than |
* . may be less than KTOP and and NDCOL may be greater than |
* . KBOT indicating phantom columns from which to chase |
* . KBOT indicating phantom columns from which to chase |
* . bulges before they are actually introduced or to which |
* . bulges before they are actually introduced or to which |
* . to chase bulges beyond column KBOT.) ==== |
* . 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 |
* ==== Bulges number MTOP to MBOT are active double implicit |
* . shift bulges. There may or may not also be small |
* . shift bulges. There may or may not also be small |
Line 403
|
Line 414
|
* . down the diagonal to make room. The phantom matrix |
* . down the diagonal to make room. The phantom matrix |
* . paradigm described above helps keep track. ==== |
* . paradigm described above helps keep track. ==== |
* |
* |
MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 ) |
MTOP = MAX( 1, ( KTOP-KRCOL ) / 2+1 ) |
MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 ) |
MBOT = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 2 ) |
M22 = MBOT + 1 |
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 ) |
$ ( KBOT-2 ) |
* |
* |
* ==== Generate reflections to chase the chain right |
* ==== Generate reflections to chase the chain right |
* . one column. (The minimum value of K is KTOP-1.) ==== |
* . one column. (The minimum value of K is KTOP-1.) ==== |
* |
* |
DO 20 M = MTOP, MBOT |
IF ( BMP22 ) THEN |
K = KRCOL + 3*( M-1 ) |
* |
|
* ==== 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 |
IF( K.EQ.KTOP-1 ) THEN |
CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ), |
CALL DLAQR1( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ), |
$ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), |
$ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ), |
Line 421
|
Line 553
|
ALPHA = V( 1, M ) |
ALPHA = V( 1, M ) |
CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) ) |
CALL DLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) ) |
ELSE |
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( 2, M ) = H( K+2, K ) |
V( 3, M ) = H( K+3, K ) |
V( 3, M ) = H( K+3, K ) |
CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) ) |
CALL DLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) ) |
Line 469
|
Line 614
|
H( K+3, K ) = ZERO |
H( K+3, K ) = ZERO |
ELSE |
ELSE |
* |
* |
* ==== Stating a new bulge here would |
* ==== Starting a new bulge here would |
* . create only negligible fill. |
* . create only negligible fill. |
* . Replace the old reflector with |
* . Replace the old reflector with |
* . the new one. ==== |
* . the new one. ==== |
Line 483
|
Line 628
|
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
20 CONTINUE |
|
* |
* |
* ==== Generate a 2-by-2 reflection, if needed. ==== |
* ==== Apply reflection from the right and |
* |
* . the first column of update from the left. |
K = KRCOL + 3*( M22-1 ) |
* . These updates are required for the vigilant |
IF( BMP22 ) THEN |
* . deflation check. We still delay most of the |
IF( K.EQ.KTOP-1 ) THEN |
* . updates from the left for efficiency. ==== |
CALL DLAQR1( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ), |
* |
$ SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ), |
T1 = V( 1, M ) |
$ V( 1, M22 ) ) |
T2 = T1*V( 2, M ) |
BETA = V( 1, M22 ) |
T3 = T1*V( 3, M ) |
CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) ) |
DO 70 J = JTOP, MIN( KBOT, K+3 ) |
ELSE |
REFSUM = H( J, K+1 ) + V( 2, M )*H( J, K+2 ) |
BETA = H( K+1, K ) |
$ + V( 3, M )*H( J, K+3 ) |
V( 2, M22 ) = H( K+2, K ) |
H( J, K+1 ) = H( J, K+1 ) - REFSUM*T1 |
CALL DLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) ) |
H( J, K+2 ) = H( J, K+2 ) - REFSUM*T2 |
H( K+1, K ) = BETA |
H( J, K+3 ) = H( J, K+3 ) - REFSUM*T3 |
H( K+2, K ) = ZERO |
70 CONTINUE |
END IF |
* |
END IF |
* ==== Perform update from left for subsequent |
* |
* . column. ==== |
* ==== Multiply H by reflections from the left ==== |
* |
* |
REFSUM = H( K+1, K+1 ) + V( 2, M )*H( K+2, K+1 ) |
IF( ACCUM ) THEN |
$ + V( 3, M )*H( K+3, K+1 ) |
JBOT = MIN( NDCOL, KBOT ) |
H( K+1, K+1 ) = H( K+1, K+1 ) - REFSUM*T1 |
ELSE IF( WANTT ) THEN |
H( K+2, K+1 ) = H( K+2, K+1 ) - REFSUM*T2 |
JBOT = N |
H( K+3, K+1 ) = H( K+3, K+1 ) - REFSUM*T3 |
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 ) ) |
|
* |
* |
* ==== The following convergence test requires that |
* ==== The following convergence test requires that |
* . the tradition small-compared-to-nearby-diagonals |
* . the tradition small-compared-to-nearby-diagonals |
Line 641
|
Line 664
|
* . is zero (as done here) is traditional but probably |
* . is zero (as done here) is traditional but probably |
* . unnecessary. ==== |
* . unnecessary. ==== |
* |
* |
|
IF( K.LT.KTOP) |
|
$ CYCLE |
IF( H( K+1, K ).NE.ZERO ) THEN |
IF( H( K+1, K ).NE.ZERO ) THEN |
TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) |
TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) |
IF( TST1.EQ.ZERO ) THEN |
IF( TST1.EQ.ZERO ) THEN |
Line 669
|
Line 694
|
TST2 = H22*( H11 / SCL ) |
TST2 = H22*( H11 / SCL ) |
* |
* |
IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. |
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 |
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 ) |
IF( ACCUM ) THEN |
DO 140 M = MTOP, MEND |
JBOT = MIN( NDCOL, KBOT ) |
K = KRCOL + 3*( M-1 ) |
ELSE IF( WANTT ) THEN |
REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 ) |
JBOT = N |
H( K+4, K+1 ) = -REFSUM |
ELSE |
H( K+4, K+2 ) = -REFSUM*V( 2, M ) |
JBOT = KBOT |
H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M ) |
END IF |
140 CONTINUE |
* |
|
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. ==== |
* ==== End of near-the-diagonal bulge chase. ==== |
* |
* |
150 CONTINUE |
145 CONTINUE |
* |
* |
* ==== Use U (if accumulated) to update far-from-diagonal |
* ==== Use U (if accumulated) to update far-from-diagonal |
* . entries in H. If required, use U to update Z as |
* . entries in H. If required, use U to update Z as |
Line 701
|
Line 787
|
JTOP = KTOP |
JTOP = KTOP |
JBOT = KBOT |
JBOT = KBOT |
END IF |
END IF |
IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR. |
K1 = MAX( 1, KTOP-INCOL ) |
$ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN |
NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1 |
|
* |
|
* ==== Horizontal Multiply ==== |
* |
* |
* ==== Updates not exploiting the 2-by-2 block |
DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH |
* . structure of U. K1 and NU keep track of |
JLEN = MIN( NH, JBOT-JCOL+1 ) |
* . the location and size of U in the special |
CALL DGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ), |
* . 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 ), |
|
$ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH, |
$ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH, |
$ LDWH ) |
$ LDWH ) |
CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH, |
CALL DLACPY( 'ALL', NU, JLEN, WH, LDWH, |
$ H( INCOL+K1, JCOL ), LDH ) |
$ H( INCOL+K1, JCOL ), LDH ) |
160 CONTINUE |
150 CONTINUE |
* |
* |
* ==== Vertical multiply ==== |
* ==== Vertical multiply ==== |
* |
* |
DO 170 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV |
DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV |
JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW ) |
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, |
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 ) |
$ LDU, ZERO, WV, LDWV ) |
CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, |
CALL DLACPY( 'ALL', JLEN, NU, WV, LDWV, |
$ H( JROW, INCOL+K1 ), LDH ) |
$ Z( JROW, INCOL+K1 ), LDZ ) |
170 CONTINUE |
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**T ==== |
|
* |
|
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**T ==== |
|
* |
|
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**T ==== |
|
* |
|
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 |
END IF |
END IF |
220 CONTINUE |
180 CONTINUE |
* |
* |
* ==== End of DLAQR5 ==== |
* ==== End of DLAQR5 ==== |
* |
* |