[BACK]Return to dlaqr5.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dlaqr5.f between versions 1.22 and 1.23

version 1.22, 2020/05/21 21:46:00 version 1.23, 2023/08/07 08:38:56
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
Line 178 Line 177
 *>  *>
 *> \param[out] U  *> \param[out] U
 *> \verbatim  *> \verbatim
 *>          U is DOUBLE PRECISION array, dimension (LDU,3*NSHFTS-3)  *>          U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS)
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDU  *> \param[in] LDU
 *> \verbatim  *> \verbatim
 *>          LDU is INTEGER  *>          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 >= 3*NSHFTS-3.  *>             in the calling subroutine.  LDU >= 2*NSHFTS.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] NV  *> \param[in] NV
Line 197 Line 196
 *>  *>
 *> \param[out] WV  *> \param[out] WV
 *> \verbatim  *> \verbatim
 *>          WV is DOUBLE PRECISION array, dimension (LDWV,3*NSHFTS-3)  *>          WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS)
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDWV  *> \param[in] LDWV
Line 223 Line 222
 *> \verbatim  *> \verbatim
 *>          LDWH is INTEGER  *>          LDWH is INTEGER
 *>             Leading dimension of WH just as declared in the  *>             Leading dimension of WH just as declared in the
 *>             calling procedure.  LDWH >= 3*NSHFTS-3.  *>             calling procedure.  LDWH >= 2*NSHFTS.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *  Authors:  *  Authors:
Line 234 Line 233
 *> \author Univ. of Colorado Denver  *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.  *> \author NAG Ltd.
 *  *
 *> \date June 2016  
 *  
 *> \ingroup doubleOTHERauxiliary  *> \ingroup doubleOTHERauxiliary
 *  *
 *> \par Contributors:  *> \par Contributors:
Line 243 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 252 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.7.1) --  *  -- 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..--
 *     June 2016  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,        INTEGER            IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
Line 280 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 356 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 371 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 401 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 419 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 467 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 481 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 639 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 667 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
   *
   *           ==== Multiply H by reflections from the left ====
 *  *
 *           ==== Fill in the last row of each bulge. ====              IF( ACCUM ) THEN
                  JBOT = MIN( NDCOL, KBOT )
               ELSE IF( WANTT ) THEN
                  JBOT = N
               ELSE
                  JBOT = KBOT
               END IF
 *  *
             MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )              DO 100 M = MBOT, MTOP, -1
             DO 140 M = MTOP, MEND                 K = KRCOL + 2*( M-1 )
                K = KRCOL + 3*( M-1 )                 T1 = V( 1, M )
                REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )                 T2 = T1*V( 2, M )
                H( K+4, K+1 ) = -REFSUM                 T3 = T1*V( 3, M )
                H( K+4, K+2 ) = -REFSUM*V( 2, M )                 DO 90 J = MAX( KTOP, KRCOL + 2*M ), JBOT
                H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M )                    REFSUM = H( K+1, J ) + V( 2, M )*H( K+2, J )
   140       CONTINUE       $                     + 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 699 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 ====
 *  *
Removed from v.1.22  
changed lines
  Added in v.1.23

CVSweb interface <joel.bertrand@systella.fr>