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

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

version 1.13, 2012/12/14 14:22:52 version 1.14, 2014/01/27 09:24:36
Line 1 Line 1
 *> \brief \b ZLASYF computes a partial factorization of a complex symmetric matrix, using the diagonal pivoting method.  *> \brief \b ZLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagonal pivoting method.
 *  *
 *  =========== 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 ZLASYF + dependencies   *> Download ZLASYF + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )  *       SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       CHARACTER          UPLO  *       CHARACTER          UPLO
 *       INTEGER            INFO, KB, LDA, LDW, N, NB  *       INTEGER            INFO, KB, LDA, LDW, N, NB
Line 28 Line 28
 *       INTEGER            IPIV( * )  *       INTEGER            IPIV( * )
 *       COMPLEX*16         A( LDA, * ), W( LDW, * )  *       COMPLEX*16         A( LDA, * ), W( LDW, * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
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 145 Line 155
 *  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 September 2012  *> \date November 2013
 *  *
 *> \ingroup complex16SYcomputational  *> \ingroup complex16SYcomputational
 *  *
   *> \par Contributors:
   *  ==================
   *>
   *> \verbatim
   *>
   *>  November 2013,  Igor Kozachenko,
   *>                  Computer Science Division,
   *>                  University of California, Berkeley
   *> \endverbatim
   *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )        SUBROUTINE ZLASYF( 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 246 Line 266
          ABSAKK = CABS1( W( K, KW ) )           ABSAKK = CABS1( 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  
 *  *
          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 257 Line 277
 *  *
          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
Line 302 Line 322
 *  *
                   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 )
                ELSE                 ELSE
Line 315 Line 335
                END IF                 END IF
             END IF              END IF
 *  *
   *           ============================================================
   *
   *           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, K ) = A( KK, K )                 A( KP, KP ) = A( KK, KK )
                CALL ZCOPY( K-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 ZCOPY( KP, 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.
 *  *
                CALL ZSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )                 IF( K.LT.N )
        $            CALL ZSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
        $                        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 )
             END IF              END IF
 *  *
             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  *              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)
 *  *
                CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )                 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
                R1 = CONE / A( K, K )                 R1 = CONE / A( K, K )
                CALL ZSCAL( K-1, R1, A( 1, K ), 1 )                 CALL ZSCAL( K-1, R1, A( 1, K ), 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
 *  *
   *              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  *                 Compose the columns of the inverse of 2-by-2 pivot
   *                 block D in the following way to reduce the number
   *                 of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
   *                 this inverse
   *
   *                 D**(-1) = ( d11 d21 )**(-1) =
   *                           ( d21 d22 )
   *
   *                 = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
   *                                        ( (-d21 ) ( d11 ) )
   *
   *                 = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
   *
   *                   * ( ( d22/d21 ) (      -1 ) ) =
   *                     ( (      -1 ) ( d11/d21 ) )
   *
   *                 = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) (  -1 ) ) =
   *                                           ( ( -1  ) ( D22 ) )
   *
   *                 = 1/d21 * T * ( ( D11 ) (  -1 ) )
   *                               ( (  -1 ) ( D22 ) )
   *
   *                 = D21 * ( ( D11 ) (  -1 ) )
   *                         ( (  -1 ) ( D22 ) )
 *  *
                   D21 = W( K-1, KW )                    D21 = W( K-1, KW )
                   D11 = W( K, KW ) / D21                    D11 = W( K, KW ) / D21
                   D22 = W( K-1, KW-1 ) / D21                    D22 = W( K-1, KW-1 ) / D21
                   T = CONE / ( D11*D22-CONE )                    T = CONE / ( D11*D22-CONE )
                   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 ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )                       A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
Line 379 Line 457
                A( K-1, K-1 ) = W( K-1, KW-1 )                 A( K-1, K-1 ) = W( K-1, KW-1 )
                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 )
   *
             END IF              END IF
   *
          END IF           END IF
 *  *
 *        Store details of the interchanges in IPIV  *        Store details of the interchanges in IPIV
Line 423 Line 503
    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 473 Line 561
          ABSAKK = CABS1( W( K, K ) )           ABSAKK = CABS1( 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  
 *  *
          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 484 Line 572
 *  *
          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
Line 528 Line 616
 *  *
                   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 )
                ELSE                 ELSE
Line 541 Line 629
                END IF                 END IF
             END IF              END IF
 *  *
   *           ============================================================
   *
   *           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, K ) = A( KK, K )                 A( KP, KP ) = A( KK, KK )
                CALL ZCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )                 CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
                CALL ZCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )       $                     LDA )
                  IF( KP.LT.N )
        $            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, 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 563 Line 665
 *  *
 *              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  *              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)
 *  *
                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
                   R1 = CONE / A( K, K )                    R1 = CONE / A( K, K )
                   CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )                    CALL ZSCAL( N-K, R1, A( 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 583 Line 691
 *              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
 *  *
   *              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  *                 Compose the columns of the inverse of 2-by-2 pivot
   *                 block D in the following way to reduce the number
   *                 of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
   *                 this inverse
   *
   *                 D**(-1) = ( d11 d21 )**(-1) =
   *                           ( d21 d22 )
   *
   *                 = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
   *                                        ( (-d21 ) ( d11 ) )
   *
   *                 = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
   *
   *                   * ( ( d22/d21 ) (      -1 ) ) =
   *                     ( (      -1 ) ( d11/d21 ) )
   *
   *                 = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) (  -1 ) ) =
   *                                           ( ( -1  ) ( D22 ) )
   *
   *                 = 1/d21 * T * ( ( D11 ) (  -1 ) )
   *                               ( (  -1 ) ( D22 ) )
   *
   *                 = D21 * ( ( D11 ) (  -1 ) )
   *                         ( (  -1 ) ( D22 ) )
 *  *
                   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 ) / D21                    D22 = W( K, K ) / D21
                   T = CONE / ( D11*D22-CONE )                    T = CONE / ( D11*D22-CONE )
                   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 ) = D21*( D11*W( J, K )-W( J, K+1 ) )                       A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )                       A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
Line 603 Line 747
                A( K, K ) = W( K, K )                 A( K, K ) = W( K, K )
                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 )
   *
             END IF              END IF
   *
          END IF           END IF
 *  *
 *        Store details of the interchanges in IPIV  *        Store details of the interchanges in IPIV
Line 648 Line 794
   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>