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

Diff for /rpl/lapack/lapack/zlahef.f between versions 1.13 and 1.14

version 1.13, 2012/12/14 14:22:50 version 1.14, 2014/01/27 09:24:36
Line 1 Line 1
 *> \brief \b ZLAHEF computes a partial factorization of a complex Hermitian indefinite matrix, using the diagonal pivoting method.  *> \brief \b ZLAHEF computes a partial factorization of a complex Hermitian indefinite matrix using the Bunch-Kaufman diagonal pivoting method (blocked algorithm, calling Level 3 BLAS).
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
Line 110 Line 110
 *> \verbatim  *> \verbatim
 *>          IPIV is INTEGER array, dimension (N)  *>          IPIV is INTEGER array, dimension (N)
 *>          Details of the interchanges and the block structure of D.  *>          Details of the interchanges and the block structure of D.
 *>          If UPLO = 'U', only the last KB elements of IPIV are set;  
 *>          if UPLO = 'L', only the first KB elements are set.  
 *>  *>
 *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were  *>          If UPLO = 'U':
 *>          interchanged and D(k,k) is a 1-by-1 diagonal block.  *>             Only the last KB elements of IPIV are set.
 *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and  *>
 *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)  *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
 *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =  *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
 *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were  *>
 *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.  *>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
   *>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
   *>             is a 2-by-2 diagonal block.
   *>
   *>          If UPLO = 'L':
   *>             Only the first KB elements of IPIV are set.
   *>
   *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
   *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
   *>
   *>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
   *>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
   *>             is a 2-by-2 diagonal block.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] W  *> \param[out] W
Line 150 Line 160
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd.   *> \author NAG Ltd. 
 *  *
 *> \date September 2012  *> \date November 2013
 *  *
 *> \ingroup complex16HEcomputational  *> \ingroup complex16HEcomputational
 *  *
   *> \par Contributors:
   *  ==================
   *>
   *> \verbatim
   *>
   *>  November 2013,  Igor Kozachenko,
   *>                  Computer Science Division,
   *>                  University of California, Berkeley
   *> \endverbatim
   *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )        SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.2) --  *  -- LAPACK computational routine (version 3.5.0) --
 *  -- 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..--
 *     September 2012  *     November 2013
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          UPLO        CHARACTER          UPLO
Line 231 Line 251
          IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )           IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
      $      GO TO 30       $      GO TO 30
 *  *
            KSTEP = 1
   *
 *        Copy column K of A to column KW of W and update it  *        Copy column K of A to column KW of W and update it
 *  *
          CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )           CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )
Line 241 Line 263
             W( K, KW ) = DBLE( W( K, KW ) )              W( K, KW ) = DBLE( W( K, KW ) )
          END IF           END IF
 *  *
          KSTEP = 1  
 *  
 *        Determine rows and columns to be interchanged and whether  *        Determine rows and columns to be interchanged and whether
 *        a 1-by-1 or 2-by-2 pivot block will be used  *        a 1-by-1 or 2-by-2 pivot block will be used
 *  *
          ABSAKK = ABS( DBLE( W( K, KW ) ) )           ABSAKK = ABS( DBLE( W( K, KW ) ) )
 *  *
 *        IMAX is the row-index of the largest off-diagonal element in  *        IMAX is the row-index of the largest off-diagonal element in
 *        column K, and COLMAX is its absolute value  *        column K, and COLMAX is its absolute value.
   *        Determine both COLMAX and IMAX.
 *  *
          IF( K.GT.1 ) THEN           IF( K.GT.1 ) THEN
             IMAX = IZAMAX( K-1, W( 1, KW ), 1 )              IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
Line 260 Line 281
 *  *
          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN           IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
 *  *
 *           Column K is zero: set INFO and continue  *           Column K is zero or underflow: set INFO and continue
 *  *
             IF( INFO.EQ.0 )              IF( INFO.EQ.0 )
      $         INFO = K       $         INFO = K
             KP = K              KP = K
             A( K, K ) = DBLE( A( K, K ) )              A( K, K ) = DBLE( A( K, K ) )
          ELSE           ELSE
   *
   *           ============================================================
   *
   *           BEGIN pivot search
   *
   *           Case(1)
             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN              IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
 *  *
 *              no interchange, use 1-by-1 pivot block  *              no interchange, use 1-by-1 pivot block
Line 274 Line 301
                KP = K                 KP = K
             ELSE              ELSE
 *  *
   *              BEGIN pivot search along IMAX row
   *
   *
 *              Copy column IMAX to column KW-1 of W and update it  *              Copy column IMAX to column KW-1 of W and update it
 *  *
                CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )                 CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
Line 289 Line 319
                END IF                 END IF
 *  *
 *              JMAX is the column-index of the largest off-diagonal  *              JMAX is the column-index of the largest off-diagonal
 *              element in row IMAX, and ROWMAX is its absolute value  *              element in row IMAX, and ROWMAX is its absolute value.
   *              Determine only ROWMAX.
 *  *
                JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )                 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
                ROWMAX = CABS1( W( JMAX, KW-1 ) )                 ROWMAX = CABS1( W( JMAX, KW-1 ) )
Line 298 Line 329
                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )                    ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
                END IF                 END IF
 *  *
   *              Case(2)
                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN                 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
 *  *
 *                 no interchange, use 1-by-1 pivot block  *                 no interchange, use 1-by-1 pivot block
 *  *
                   KP = K                    KP = K
   *
   *              Case(3)
                ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )                 ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )
      $                   THEN       $                   THEN
 *  *
Line 311 Line 345
 *  *
                   KP = IMAX                    KP = IMAX
 *  *
 *                 copy column KW-1 of W to column KW  *                 copy column KW-1 of W to column KW of W
 *  *
                   CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )                    CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
   *
   *              Case(4)
                ELSE                 ELSE
 *  *
 *                 interchange rows and columns K-1 and IMAX, use 2-by-2  *                 interchange rows and columns K-1 and IMAX, use 2-by-2
Line 322 Line 358
                   KP = IMAX                    KP = IMAX
                   KSTEP = 2                    KSTEP = 2
                END IF                 END IF
   *
   *
   *              END pivot search along IMAX row
   *
             END IF              END IF
 *  *
   *           END pivot search
   *
   *           ============================================================
   *
   *           KK is the column of A where pivoting step stopped
   *
             KK = K - KSTEP + 1              KK = K - KSTEP + 1
   *
   *           KKW is the column of W which corresponds to column KK of A
   *
             KKW = NB + KK - N              KKW = NB + KK - N
 *  *
 *           Updated column KP is already stored in column KKW of W  *           Interchange rows and columns KP and KK.
   *           Updated column KP is already stored in column KKW of W.
 *  *
             IF( KP.NE.KK ) THEN              IF( KP.NE.KK ) THEN
 *  *
 *              Copy non-updated column KK to column KP  *              Copy non-updated column KK to column KP of submatrix A
   *              at step K. No need to copy element into column K
   *              (or K and K-1 for 2-by-2 pivot) of A, since these columns
   *              will be later overwritten.
 *  *
                A( KP, KP ) = DBLE( A( KK, KK ) )                 A( KP, KP ) = DBLE( A( KK, KK ) )
                CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),                 CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
      $                     LDA )       $                     LDA )
                CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA )                 CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA )
                CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )                 IF( KP.GT.1 )
        $            CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
 *  *
 *              Interchange rows KK and KP in last KK columns of A and W  *              Interchange rows KK and KP in last K+1 to N columns of A
   *              (columns K (or K and K-1 for 2-by-2 pivot) of A will be
   *              later overwritten). Interchange rows KK and KP
   *              in last KKW to NB columns of W.
 *  *
                IF( KK.LT.N )                 IF( K.LT.N )
      $            CALL ZSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ),       $            CALL ZSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
      $                        LDA )       $                        LDA )
                CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),                 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
      $                     LDW )       $                     LDW )
Line 350 Line 407
 *  *
             IF( KSTEP.EQ.1 ) THEN              IF( KSTEP.EQ.1 ) THEN
 *  *
 *              1-by-1 pivot block D(k): column KW of W now holds  *              1-by-1 pivot block D(k): column kw of W now holds
 *  *
 *              W(k) = U(k)*D(k)  *              W(kw) = U(k)*D(k),
 *  *
 *              where U(k) is the k-th column of U  *              where U(k) is the k-th column of U
 *  *
 *              Store U(k) in column k of A  *              (1) Store subdiag. elements of column U(k)
 *  *              and 1-by-1 block D(k) in column k of A.
   *              (NOTE: Diagonal element U(k,k) is a UNIT element
   *              and not stored)
   *                 A(k,k) := D(k,k) = W(k,kw)
   *                 A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
   *
   *              (NOTE: No need to use for Hermitian matrix
   *              A( K, K ) = DBLE( W( K, K) ) to separately copy diagonal
   *              element D(k,k) from W (potentially saves only one load))
                CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )                 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
                R1 = ONE / DBLE( A( K, K ) )                 IF( K.GT.1 ) THEN
                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )  
 *  *
 *              Conjugate W(k)  *                 (NOTE: No need to check if A(k,k) is NOT ZERO,
   *                  since that was ensured earlier in pivot search:
   *                  case A(k,k) = 0 falls into 2x2 pivot case(4))
   *
                     R1 = ONE / DBLE( A( K, K ) )
                     CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
   *
   *                 (2) Conjugate column W(kw)
   *
                     CALL ZLACGV( K-1, W( 1, KW ), 1 )
                  END IF
 *  *
                CALL ZLACGV( K-1, W( 1, KW ), 1 )  
             ELSE              ELSE
 *  *
 *              2-by-2 pivot block D(k): columns KW and KW-1 of W now  *              2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
 *              hold  
 *  *
 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)  *              ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
 *  *
 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns  *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
 *              of U  *              of U
 *  *
   *              (1) Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
   *              block D(k-1:k,k-1:k) in columns k-1 and k of A.
   *              (NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
   *              block and not stored)
   *                 A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
   *                 A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
   *                 = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
   *
                IF( K.GT.2 ) THEN                 IF( K.GT.2 ) THEN
 *  *
 *                 Store U(k) and U(k-1) in columns k and k-1 of A  *                 Factor out the columns of the inverse of 2-by-2 pivot
   *                 block D, so that each column contains 1, to reduce the
   *                 number of FLOPS when we multiply panel
   *                 ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
   *
   *                 D**(-1) = ( d11 cj(d21) )**(-1) =
   *                           ( d21    d22 )
   *
   *                 = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
   *                                          ( (-d21) (     d11 ) )
   *
   *                 = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
   *
   *                   * ( d21*( d22/d21 ) conj(d21)*(           - 1 ) ) =
   *                     (     (      -1 )           ( d11/conj(d21) ) )
   *
   *                 = 1/(|d21|**2) * 1/(D22*D11-1) *
   *
   *                   * ( d21*( D11 ) conj(d21)*(  -1 ) ) =
   *                     (     (  -1 )           ( D22 ) )
   *
   *                 = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*(  -1 ) ) =
   *                                      (     (  -1 )           ( D22 ) )
   *
   *                 = ( (T/conj(d21))*( D11 ) (T/d21)*(  -1 ) ) =
   *                   (               (  -1 )         ( D22 ) )
   *
   *                 = ( conj(D21)*( D11 ) D21*(  -1 ) )
   *                   (           (  -1 )     ( D22 ) ),
   *
   *                 where D11 = d22/d21,
   *                       D22 = d11/conj(d21),
   *                       D21 = T/d21,
   *                       T = 1/(D22*D11-1).
   *
   *                 (NOTE: No need to check for division by ZERO,
   *                  since that was ensured earlier in pivot search:
   *                  (a) d21 != 0, since in 2x2 pivot case(4)
   *                      |d21| should be larger than |d11| and |d22|;
   *                  (b) (D22*D11 - 1) != 0, since from (a),
   *                      both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
 *  *
                   D21 = W( K-1, KW )                    D21 = W( K-1, KW )
                   D11 = W( K, KW ) / DCONJG( D21 )                    D11 = W( K, KW ) / DCONJG( D21 )
                   D22 = W( K-1, KW-1 ) / D21                    D22 = W( K-1, KW-1 ) / D21
                   T = ONE / ( DBLE( D11*D22 )-ONE )                    T = ONE / ( DBLE( D11*D22 )-ONE )
                   D21 = T / D21                    D21 = T / D21
   *
   *                 Update elements in columns A(k-1) and A(k) as
   *                 dot products of rows of ( W(kw-1) W(kw) ) and columns
   *                 of D**(-1)
   *
                   DO 20 J = 1, K - 2                    DO 20 J = 1, K - 2
                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )                       A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
                      A( J, K ) = DCONJG( D21 )*                       A( J, K ) = DCONJG( D21 )*
Line 397 Line 522
                A( K-1, K ) = W( K-1, KW )                 A( K-1, K ) = W( K-1, KW )
                A( K, K ) = W( K, KW )                 A( K, K ) = W( K, KW )
 *  *
 *              Conjugate W(k) and W(k-1)  *              (2) Conjugate columns W(kw) and W(kw-1)
 *  *
                CALL ZLACGV( K-1, W( 1, KW ), 1 )                 CALL ZLACGV( K-1, W( 1, KW ), 1 )
                CALL ZLACGV( K-2, W( 1, KW-1 ), 1 )                 CALL ZLACGV( K-2, W( 1, KW-1 ), 1 )
   *
             END IF              END IF
   *
          END IF           END IF
 *  *
 *        Store details of the interchanges in IPIV  *        Store details of the interchanges in IPIV
Line 448 Line 575
    50    CONTINUE     50    CONTINUE
 *  *
 *        Put U12 in standard form by partially undoing the interchanges  *        Put U12 in standard form by partially undoing the interchanges
 *        in columns k+1:n  *        in columns k+1:n looping backwards from k+1 to n
 *  *
          J = K + 1           J = K + 1
    60    CONTINUE     60    CONTINUE
          JJ = J  *
          JP = IPIV( J )  *           Undo the interchanges (if any) of rows JJ and JP at each
          IF( JP.LT.0 ) THEN  *           step J
             JP = -JP  *
   *           (Here, J is a diagonal index)
               JJ = J
               JP = IPIV( J )
               IF( JP.LT.0 ) THEN
                  JP = -JP
   *              (Here, J is a diagonal index)
                  J = J + 1
               END IF
   *           (NOTE: Here, J is used to determine row length. Length N-J+1
   *           of the rows to swap back doesn't include diagonal element)
             J = J + 1              J = J + 1
          END IF              IF( JP.NE.JJ .AND. J.LE.N )
          J = J + 1       $         CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
          IF( JP.NE.JJ .AND. J.LE.N )           IF( J.LT.N )
      $      CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )  
          IF( J.LE.N )  
      $      GO TO 60       $      GO TO 60
 *  *
 *        Set KB to the number of columns factorized  *        Set KB to the number of columns factorized
Line 484 Line 619
          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )           IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
      $      GO TO 90       $      GO TO 90
 *  *
            KSTEP = 1
   *
 *        Copy column K of A to column K of W and update it  *        Copy column K of A to column K of W and update it
 *  *
          W( K, K ) = DBLE( A( K, K ) )           W( K, K ) = DBLE( A( K, K ) )
Line 493 Line 630
      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )       $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
          W( K, K ) = DBLE( W( K, K ) )           W( K, K ) = DBLE( W( K, K ) )
 *  *
          KSTEP = 1  
 *  
 *        Determine rows and columns to be interchanged and whether  *        Determine rows and columns to be interchanged and whether
 *        a 1-by-1 or 2-by-2 pivot block will be used  *        a 1-by-1 or 2-by-2 pivot block will be used
 *  *
          ABSAKK = ABS( DBLE( W( K, K ) ) )           ABSAKK = ABS( DBLE( W( K, K ) ) )
 *  *
 *        IMAX is the row-index of the largest off-diagonal element in  *        IMAX is the row-index of the largest off-diagonal element in
 *        column K, and COLMAX is its absolute value  *        column K, and COLMAX is its absolute value.
   *        Determine both COLMAX and IMAX.
 *  *
          IF( K.LT.N ) THEN           IF( K.LT.N ) THEN
             IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )              IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
Line 512 Line 648
 *  *
          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN           IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
 *  *
 *           Column K is zero: set INFO and continue  *           Column K is zero or underflow: set INFO and continue
 *  *
             IF( INFO.EQ.0 )              IF( INFO.EQ.0 )
      $         INFO = K       $         INFO = K
             KP = K              KP = K
             A( K, K ) = DBLE( A( K, K ) )              A( K, K ) = DBLE( A( K, K ) )
          ELSE           ELSE
   *
   *           ============================================================
   *
   *           BEGIN pivot search
   *
   *           Case(1)
             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN              IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
 *  *
 *              no interchange, use 1-by-1 pivot block  *              no interchange, use 1-by-1 pivot block
Line 526 Line 668
                KP = K                 KP = K
             ELSE              ELSE
 *  *
   *              BEGIN pivot search along IMAX row
   *
   *
 *              Copy column IMAX to column K+1 of W and update it  *              Copy column IMAX to column K+1 of W and update it
 *  *
                CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )                 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
Line 540 Line 685
                W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) )                 W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) )
 *  *
 *              JMAX is the column-index of the largest off-diagonal  *              JMAX is the column-index of the largest off-diagonal
 *              element in row IMAX, and ROWMAX is its absolute value  *              element in row IMAX, and ROWMAX is its absolute value.
   *              Determine only ROWMAX.
 *  *
                JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )                 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
                ROWMAX = CABS1( W( JMAX, K+1 ) )                 ROWMAX = CABS1( W( JMAX, K+1 ) )
Line 549 Line 695
                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )                    ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
                END IF                 END IF
 *  *
   *              Case(2)
                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN                 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
 *  *
 *                 no interchange, use 1-by-1 pivot block  *                 no interchange, use 1-by-1 pivot block
 *  *
                   KP = K                    KP = K
   *
   *              Case(3)
                ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )                 ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )
      $                   THEN       $                   THEN
 *  *
Line 562 Line 711
 *  *
                   KP = IMAX                    KP = IMAX
 *  *
 *                 copy column K+1 of W to column K  *                 copy column K+1 of W to column K of W
 *  *
                   CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )                    CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
   *
   *              Case(4)
                ELSE                 ELSE
 *  *
 *                 interchange rows and columns K+1 and IMAX, use 2-by-2  *                 interchange rows and columns K+1 and IMAX, use 2-by-2
Line 573 Line 724
                   KP = IMAX                    KP = IMAX
                   KSTEP = 2                    KSTEP = 2
                END IF                 END IF
   *
   *
   *              END pivot search along IMAX row
   *
             END IF              END IF
 *  *
   *           END pivot search
   *
   *           ============================================================
   *
   *           KK is the column of A where pivoting step stopped
   *
             KK = K + KSTEP - 1              KK = K + KSTEP - 1
 *  *
 *           Updated column KP is already stored in column KK of W  *           Interchange rows and columns KP and KK.
   *           Updated column KP is already stored in column KK of W.
 *  *
             IF( KP.NE.KK ) THEN              IF( KP.NE.KK ) THEN
 *  *
 *              Copy non-updated column KK to column KP  *              Copy non-updated column KK to column KP of submatrix A
   *              at step K. No need to copy element into column K
   *              (or K and K+1 for 2-by-2 pivot) of A, since these columns
   *              will be later overwritten.
 *  *
                A( KP, KP ) = DBLE( A( KK, KK ) )                 A( KP, KP ) = DBLE( A( KK, KK ) )
                CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),                 CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
Line 590 Line 755
                IF( KP.LT.N )                 IF( KP.LT.N )
      $            CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )       $            CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
 *  *
 *              Interchange rows KK and KP in first KK columns of A and W  *              Interchange rows KK and KP in first K-1 columns of A
   *              (columns K (or K and K+1 for 2-by-2 pivot) of A will be
   *              later overwritten). Interchange rows KK and KP
   *              in first KK columns of W.
 *  *
                CALL ZSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )                 IF( K.GT.1 )
        $            CALL ZSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
                CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )                 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
             END IF              END IF
 *  *
Line 600 Line 769
 *  *
 *              1-by-1 pivot block D(k): column k of W now holds  *              1-by-1 pivot block D(k): column k of W now holds
 *  *
 *              W(k) = L(k)*D(k)  *              W(k) = L(k)*D(k),
 *  *
 *              where L(k) is the k-th column of L  *              where L(k) is the k-th column of L
 *  *
 *              Store L(k) in column k of A  *              (1) Store subdiag. elements of column L(k)
 *  *              and 1-by-1 block D(k) in column k of A.
   *              (NOTE: Diagonal element L(k,k) is a UNIT element
   *              and not stored)
   *                 A(k,k) := D(k,k) = W(k,k)
   *                 A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
   *
   *              (NOTE: No need to use for Hermitian matrix
   *              A( K, K ) = DBLE( W( K, K) ) to separately copy diagonal
   *              element D(k,k) from W (potentially saves only one load))
                CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )                 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
                IF( K.LT.N ) THEN                 IF( K.LT.N ) THEN
   *
   *                 (NOTE: No need to check if A(k,k) is NOT ZERO,
   *                  since that was ensured earlier in pivot search:
   *                  case A(k,k) = 0 falls into 2x2 pivot case(4))
   *
                   R1 = ONE / DBLE( A( K, K ) )                    R1 = ONE / DBLE( A( K, K ) )
                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )                    CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
 *  *
 *                 Conjugate W(k)  *                 (2) Conjugate column W(k)
 *  *
                   CALL ZLACGV( N-K, W( K+1, K ), 1 )                    CALL ZLACGV( N-K, W( K+1, K ), 1 )
                END IF                 END IF
   *
             ELSE              ELSE
 *  *
 *              2-by-2 pivot block D(k): columns k and k+1 of W now hold  *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
Line 624 Line 807
 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns  *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
 *              of L  *              of L
 *  *
   *              (1) Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
   *              block D(k:k+1,k:k+1) in columns k and k+1 of A.
   *              (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
   *              block and not stored)
   *                 A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
   *                 A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
   *                 = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
   *
                IF( K.LT.N-1 ) THEN                 IF( K.LT.N-1 ) THEN
 *  *
 *                 Store L(k) and L(k+1) in columns k and k+1 of A  *                 Factor out the columns of the inverse of 2-by-2 pivot
   *                 block D, so that each column contains 1, to reduce the
   *                 number of FLOPS when we multiply panel
   *                 ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
   *
   *                 D**(-1) = ( d11 cj(d21) )**(-1) =
   *                           ( d21    d22 )
   *
   *                 = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
   *                                          ( (-d21) (     d11 ) )
   *
   *                 = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
   *
   *                   * ( d21*( d22/d21 ) conj(d21)*(           - 1 ) ) =
   *                     (     (      -1 )           ( d11/conj(d21) ) )
   *
   *                 = 1/(|d21|**2) * 1/(D22*D11-1) *
   *
   *                   * ( d21*( D11 ) conj(d21)*(  -1 ) ) =
   *                     (     (  -1 )           ( D22 ) )
   *
   *                 = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*(  -1 ) ) =
   *                                      (     (  -1 )           ( D22 ) )
   *
   *                 = ( (T/conj(d21))*( D11 ) (T/d21)*(  -1 ) ) =
   *                   (               (  -1 )         ( D22 ) )
   *
   *                 = ( conj(D21)*( D11 ) D21*(  -1 ) )
   *                   (           (  -1 )     ( D22 ) ),
   *
   *                 where D11 = d22/d21,
   *                       D22 = d11/conj(d21),
   *                       D21 = T/d21,
   *                       T = 1/(D22*D11-1).
   *
   *                 (NOTE: No need to check for division by ZERO,
   *                  since that was ensured earlier in pivot search:
   *                  (a) d21 != 0, since in 2x2 pivot case(4)
   *                      |d21| should be larger than |d11| and |d22|;
   *                  (b) (D22*D11 - 1) != 0, since from (a),
   *                      both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
 *  *
                   D21 = W( K+1, K )                    D21 = W( K+1, K )
                   D11 = W( K+1, K+1 ) / D21                    D11 = W( K+1, K+1 ) / D21
                   D22 = W( K, K ) / DCONJG( D21 )                    D22 = W( K, K ) / DCONJG( D21 )
                   T = ONE / ( DBLE( D11*D22 )-ONE )                    T = ONE / ( DBLE( D11*D22 )-ONE )
                   D21 = T / D21                    D21 = T / D21
   *
   *                 Update elements in columns A(k) and A(k+1) as
   *                 dot products of rows of ( W(k) W(k+1) ) and columns
   *                 of D**(-1)
   *
                   DO 80 J = K + 2, N                    DO 80 J = K + 2, N
                      A( J, K ) = DCONJG( D21 )*                       A( J, K ) = DCONJG( D21 )*
      $                           ( D11*W( J, K )-W( J, K+1 ) )       $                           ( D11*W( J, K )-W( J, K+1 ) )
Line 646 Line 882
                A( K+1, K ) = W( K+1, K )                 A( K+1, K ) = W( K+1, K )
                A( K+1, K+1 ) = W( K+1, K+1 )                 A( K+1, K+1 ) = W( K+1, K+1 )
 *  *
 *              Conjugate W(k) and W(k+1)  *              (2) Conjugate columns W(k) and W(k+1)
 *  *
                CALL ZLACGV( N-K, W( K+1, K ), 1 )                 CALL ZLACGV( N-K, W( K+1, K ), 1 )
                CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 )                 CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 )
   *
             END IF              END IF
   *
          END IF           END IF
 *  *
 *        Store details of the interchanges in IPIV  *        Store details of the interchanges in IPIV
Line 698 Line 936
   110    CONTINUE    110    CONTINUE
 *  *
 *        Put L21 in standard form by partially undoing the interchanges  *        Put L21 in standard form by partially undoing the interchanges
 *        in columns 1:k-1  *        of rows in columns 1:k-1 looping backwards from k-1 to 1
 *  *
          J = K - 1           J = K - 1
   120    CONTINUE    120    CONTINUE
          JJ = J  *
          JP = IPIV( J )  *           Undo the interchanges (if any) of rows JJ and JP at each
          IF( JP.LT.0 ) THEN  *           step J
             JP = -JP  *
   *           (Here, J is a diagonal index)
               JJ = J
               JP = IPIV( J )
               IF( JP.LT.0 ) THEN
                  JP = -JP
   *              (Here, J is a diagonal index)
                  J = J - 1
               END IF
   *           (NOTE: Here, J is used to determine row length. Length J
   *           of the rows to swap back doesn't include diagonal element)
             J = J - 1              J = J - 1
          END IF              IF( JP.NE.JJ .AND. J.GE.1 )
          J = J - 1       $         CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
          IF( JP.NE.JJ .AND. J.GE.1 )           IF( J.GT.1 )
      $      CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )  
          IF( J.GE.1 )  
      $      GO TO 120       $      GO TO 120
 *  *
 *        Set KB to the number of columns factorized  *        Set KB to the number of columns factorized
Removed from v.1.13  
changed lines
  Added in v.1.14

CVSweb interface <joel.bertrand@systella.fr>