Annotation of rpl/lapack/lapack/zsytf2.f, revision 1.14

1.12      bertrand    1: *> \brief \b ZSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
1.9       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download ZSYTF2 + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO )
                     22: * 
                     23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          UPLO
                     25: *       INTEGER            INFO, LDA, N
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       INTEGER            IPIV( * )
                     29: *       COMPLEX*16         A( LDA, * )
                     30: *       ..
                     31: *  
                     32: *
                     33: *> \par Purpose:
                     34: *  =============
                     35: *>
                     36: *> \verbatim
                     37: *>
                     38: *> ZSYTF2 computes the factorization of a complex symmetric matrix A
                     39: *> using the Bunch-Kaufman diagonal pivoting method:
                     40: *>
                     41: *>    A = U*D*U**T  or  A = L*D*L**T
                     42: *>
                     43: *> where U (or L) is a product of permutation and unit upper (lower)
                     44: *> triangular matrices, U**T is the transpose of U, and D is symmetric and
                     45: *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
                     46: *>
                     47: *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
                     48: *> \endverbatim
                     49: *
                     50: *  Arguments:
                     51: *  ==========
                     52: *
                     53: *> \param[in] UPLO
                     54: *> \verbatim
                     55: *>          UPLO is CHARACTER*1
                     56: *>          Specifies whether the upper or lower triangular part of the
                     57: *>          symmetric matrix A is stored:
                     58: *>          = 'U':  Upper triangular
                     59: *>          = 'L':  Lower triangular
                     60: *> \endverbatim
                     61: *>
                     62: *> \param[in] N
                     63: *> \verbatim
                     64: *>          N is INTEGER
                     65: *>          The order of the matrix A.  N >= 0.
                     66: *> \endverbatim
                     67: *>
                     68: *> \param[in,out] A
                     69: *> \verbatim
                     70: *>          A is COMPLEX*16 array, dimension (LDA,N)
                     71: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
                     72: *>          n-by-n upper triangular part of A contains the upper
                     73: *>          triangular part of the matrix A, and the strictly lower
                     74: *>          triangular part of A is not referenced.  If UPLO = 'L', the
                     75: *>          leading n-by-n lower triangular part of A contains the lower
                     76: *>          triangular part of the matrix A, and the strictly upper
                     77: *>          triangular part of A is not referenced.
                     78: *>
                     79: *>          On exit, the block diagonal matrix D and the multipliers used
                     80: *>          to obtain the factor U or L (see below for further details).
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[in] LDA
                     84: *> \verbatim
                     85: *>          LDA is INTEGER
                     86: *>          The leading dimension of the array A.  LDA >= max(1,N).
                     87: *> \endverbatim
                     88: *>
                     89: *> \param[out] IPIV
                     90: *> \verbatim
                     91: *>          IPIV is INTEGER array, dimension (N)
                     92: *>          Details of the interchanges and the block structure of D.
1.14    ! bertrand   93: *>
        !            94: *>          If UPLO = 'U':
        !            95: *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
        !            96: *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
        !            97: *>
        !            98: *>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
        !            99: *>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
        !           100: *>             is a 2-by-2 diagonal block.
        !           101: *>
        !           102: *>          If UPLO = 'L':
        !           103: *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
        !           104: *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
        !           105: *>
        !           106: *>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
        !           107: *>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
        !           108: *>             is a 2-by-2 diagonal block.
1.9       bertrand  109: *> \endverbatim
                    110: *>
                    111: *> \param[out] INFO
                    112: *> \verbatim
                    113: *>          INFO is INTEGER
                    114: *>          = 0: successful exit
                    115: *>          < 0: if INFO = -k, the k-th argument had an illegal value
                    116: *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
                    117: *>               has been completed, but the block diagonal matrix D is
                    118: *>               exactly singular, and division by zero will occur if it
                    119: *>               is used to solve a system of equations.
                    120: *> \endverbatim
                    121: *
                    122: *  Authors:
                    123: *  ========
                    124: *
                    125: *> \author Univ. of Tennessee 
                    126: *> \author Univ. of California Berkeley 
                    127: *> \author Univ. of Colorado Denver 
                    128: *> \author NAG Ltd. 
                    129: *
1.14    ! bertrand  130: *> \date November 2013
1.9       bertrand  131: *
                    132: *> \ingroup complex16SYcomputational
                    133: *
                    134: *> \par Further Details:
                    135: *  =====================
                    136: *>
                    137: *> \verbatim
                    138: *>
                    139: *>  If UPLO = 'U', then A = U*D*U**T, where
                    140: *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
                    141: *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
                    142: *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
                    143: *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
                    144: *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
                    145: *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
                    146: *>
                    147: *>             (   I    v    0   )   k-s
                    148: *>     U(k) =  (   0    I    0   )   s
                    149: *>             (   0    0    I   )   n-k
                    150: *>                k-s   s   n-k
                    151: *>
                    152: *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
                    153: *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
                    154: *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
                    155: *>
                    156: *>  If UPLO = 'L', then A = L*D*L**T, where
                    157: *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
                    158: *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
                    159: *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
                    160: *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
                    161: *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
                    162: *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
                    163: *>
                    164: *>             (   I    0     0   )  k-1
                    165: *>     L(k) =  (   0    I     0   )  s
                    166: *>             (   0    v     I   )  n-k-s+1
                    167: *>                k-1   s  n-k-s+1
                    168: *>
                    169: *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
                    170: *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
                    171: *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
                    172: *> \endverbatim
                    173: *
                    174: *> \par Contributors:
                    175: *  ==================
                    176: *>
                    177: *> \verbatim
                    178: *>
                    179: *>  09-29-06 - patch from
                    180: *>    Bobby Cheng, MathWorks
                    181: *>
                    182: *>    Replace l.209 and l.377
                    183: *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
                    184: *>    by
                    185: *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
                    186: *>
                    187: *>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
                    188: *>         Company
                    189: *> \endverbatim
                    190: *
                    191: *  =====================================================================
1.1       bertrand  192:       SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO )
                    193: *
1.14    ! bertrand  194: *  -- LAPACK computational routine (version 3.5.0) --
1.1       bertrand  195: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    196: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.14    ! bertrand  197: *     November 2013
1.1       bertrand  198: *
                    199: *     .. Scalar Arguments ..
                    200:       CHARACTER          UPLO
                    201:       INTEGER            INFO, LDA, N
                    202: *     ..
                    203: *     .. Array Arguments ..
                    204:       INTEGER            IPIV( * )
                    205:       COMPLEX*16         A( LDA, * )
                    206: *     ..
                    207: *
                    208: *  =====================================================================
                    209: *
                    210: *     .. Parameters ..
                    211:       DOUBLE PRECISION   ZERO, ONE
                    212:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
                    213:       DOUBLE PRECISION   EIGHT, SEVTEN
                    214:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
                    215:       COMPLEX*16         CONE
                    216:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
                    217: *     ..
                    218: *     .. Local Scalars ..
                    219:       LOGICAL            UPPER
                    220:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
                    221:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, ROWMAX
                    222:       COMPLEX*16         D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
                    223: *     ..
                    224: *     .. External Functions ..
                    225:       LOGICAL            DISNAN, LSAME
                    226:       INTEGER            IZAMAX
                    227:       EXTERNAL           DISNAN, LSAME, IZAMAX
                    228: *     ..
                    229: *     .. External Subroutines ..
                    230:       EXTERNAL           XERBLA, ZSCAL, ZSWAP, ZSYR
                    231: *     ..
                    232: *     .. Intrinsic Functions ..
                    233:       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
                    234: *     ..
                    235: *     .. Statement Functions ..
                    236:       DOUBLE PRECISION   CABS1
                    237: *     ..
                    238: *     .. Statement Function definitions ..
                    239:       CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
                    240: *     ..
                    241: *     .. Executable Statements ..
                    242: *
                    243: *     Test the input parameters.
                    244: *
                    245:       INFO = 0
                    246:       UPPER = LSAME( UPLO, 'U' )
                    247:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                    248:          INFO = -1
                    249:       ELSE IF( N.LT.0 ) THEN
                    250:          INFO = -2
                    251:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    252:          INFO = -4
                    253:       END IF
                    254:       IF( INFO.NE.0 ) THEN
                    255:          CALL XERBLA( 'ZSYTF2', -INFO )
                    256:          RETURN
                    257:       END IF
                    258: *
                    259: *     Initialize ALPHA for use in choosing pivot block size.
                    260: *
                    261:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
                    262: *
                    263:       IF( UPPER ) THEN
                    264: *
1.8       bertrand  265: *        Factorize A as U*D*U**T using the upper triangle of A
1.1       bertrand  266: *
                    267: *        K is the main loop index, decreasing from N to 1 in steps of
                    268: *        1 or 2
                    269: *
                    270:          K = N
                    271:    10    CONTINUE
                    272: *
                    273: *        If K < 1, exit from loop
                    274: *
                    275:          IF( K.LT.1 )
                    276:      $      GO TO 70
                    277:          KSTEP = 1
                    278: *
                    279: *        Determine rows and columns to be interchanged and whether
                    280: *        a 1-by-1 or 2-by-2 pivot block will be used
                    281: *
                    282:          ABSAKK = CABS1( A( K, K ) )
                    283: *
                    284: *        IMAX is the row-index of the largest off-diagonal element in
1.14    ! bertrand  285: *        column K, and COLMAX is its absolute value.
        !           286: *        Determine both COLMAX and IMAX.
1.1       bertrand  287: *
                    288:          IF( K.GT.1 ) THEN
                    289:             IMAX = IZAMAX( K-1, A( 1, K ), 1 )
                    290:             COLMAX = CABS1( A( IMAX, K ) )
                    291:          ELSE
                    292:             COLMAX = ZERO
                    293:          END IF
                    294: *
                    295:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN
                    296: *
1.14    ! bertrand  297: *           Column K is zero or underflow, or contains a NaN:
        !           298: *           set INFO and continue
1.1       bertrand  299: *
                    300:             IF( INFO.EQ.0 )
                    301:      $         INFO = K
                    302:             KP = K
                    303:          ELSE
                    304:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
                    305: *
                    306: *              no interchange, use 1-by-1 pivot block
                    307: *
                    308:                KP = K
                    309:             ELSE
                    310: *
                    311: *              JMAX is the column-index of the largest off-diagonal
                    312: *              element in row IMAX, and ROWMAX is its absolute value
                    313: *
                    314:                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
                    315:                ROWMAX = CABS1( A( IMAX, JMAX ) )
                    316:                IF( IMAX.GT.1 ) THEN
                    317:                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
                    318:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
                    319:                END IF
                    320: *
                    321:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
                    322: *
                    323: *                 no interchange, use 1-by-1 pivot block
                    324: *
                    325:                   KP = K
                    326:                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
                    327: *
                    328: *                 interchange rows and columns K and IMAX, use 1-by-1
                    329: *                 pivot block
                    330: *
                    331:                   KP = IMAX
                    332:                ELSE
                    333: *
                    334: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
                    335: *                 pivot block
                    336: *
                    337:                   KP = IMAX
                    338:                   KSTEP = 2
                    339:                END IF
                    340:             END IF
                    341: *
                    342:             KK = K - KSTEP + 1
                    343:             IF( KP.NE.KK ) THEN
                    344: *
                    345: *              Interchange rows and columns KK and KP in the leading
                    346: *              submatrix A(1:k,1:k)
                    347: *
                    348:                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
                    349:                CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
                    350:      $                     LDA )
                    351:                T = A( KK, KK )
                    352:                A( KK, KK ) = A( KP, KP )
                    353:                A( KP, KP ) = T
                    354:                IF( KSTEP.EQ.2 ) THEN
                    355:                   T = A( K-1, K )
                    356:                   A( K-1, K ) = A( KP, K )
                    357:                   A( KP, K ) = T
                    358:                END IF
                    359:             END IF
                    360: *
                    361: *           Update the leading submatrix
                    362: *
                    363:             IF( KSTEP.EQ.1 ) THEN
                    364: *
                    365: *              1-by-1 pivot block D(k): column k now holds
                    366: *
                    367: *              W(k) = U(k)*D(k)
                    368: *
                    369: *              where U(k) is the k-th column of U
                    370: *
                    371: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
                    372: *
1.8       bertrand  373: *              A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
1.1       bertrand  374: *
                    375:                R1 = CONE / A( K, K )
                    376:                CALL ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
                    377: *
                    378: *              Store U(k) in column k
                    379: *
                    380:                CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
                    381:             ELSE
                    382: *
                    383: *              2-by-2 pivot block D(k): columns k and k-1 now hold
                    384: *
                    385: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
                    386: *
                    387: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
                    388: *              of U
                    389: *
                    390: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
                    391: *
1.8       bertrand  392: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
                    393: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
1.1       bertrand  394: *
                    395:                IF( K.GT.2 ) THEN
                    396: *
                    397:                   D12 = A( K-1, K )
                    398:                   D22 = A( K-1, K-1 ) / D12
                    399:                   D11 = A( K, K ) / D12
                    400:                   T = CONE / ( D11*D22-CONE )
                    401:                   D12 = T / D12
                    402: *
                    403:                   DO 30 J = K - 2, 1, -1
                    404:                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
                    405:                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
                    406:                      DO 20 I = J, 1, -1
                    407:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
                    408:      $                              A( I, K-1 )*WKM1
                    409:    20                CONTINUE
                    410:                      A( J, K ) = WK
                    411:                      A( J, K-1 ) = WKM1
                    412:    30             CONTINUE
                    413: *
                    414:                END IF
                    415: *
                    416:             END IF
                    417:          END IF
                    418: *
                    419: *        Store details of the interchanges in IPIV
                    420: *
                    421:          IF( KSTEP.EQ.1 ) THEN
                    422:             IPIV( K ) = KP
                    423:          ELSE
                    424:             IPIV( K ) = -KP
                    425:             IPIV( K-1 ) = -KP
                    426:          END IF
                    427: *
                    428: *        Decrease K and return to the start of the main loop
                    429: *
                    430:          K = K - KSTEP
                    431:          GO TO 10
                    432: *
                    433:       ELSE
                    434: *
1.8       bertrand  435: *        Factorize A as L*D*L**T using the lower triangle of A
1.1       bertrand  436: *
                    437: *        K is the main loop index, increasing from 1 to N in steps of
                    438: *        1 or 2
                    439: *
                    440:          K = 1
                    441:    40    CONTINUE
                    442: *
                    443: *        If K > N, exit from loop
                    444: *
                    445:          IF( K.GT.N )
                    446:      $      GO TO 70
                    447:          KSTEP = 1
                    448: *
                    449: *        Determine rows and columns to be interchanged and whether
                    450: *        a 1-by-1 or 2-by-2 pivot block will be used
                    451: *
                    452:          ABSAKK = CABS1( A( K, K ) )
                    453: *
                    454: *        IMAX is the row-index of the largest off-diagonal element in
1.14    ! bertrand  455: *        column K, and COLMAX is its absolute value.
        !           456: *        Determine both COLMAX and IMAX.
1.1       bertrand  457: *
                    458:          IF( K.LT.N ) THEN
                    459:             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
                    460:             COLMAX = CABS1( A( IMAX, K ) )
                    461:          ELSE
                    462:             COLMAX = ZERO
                    463:          END IF
                    464: *
                    465:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN
                    466: *
1.14    ! bertrand  467: *           Column K is zero or underflow, or contains a NaN:
        !           468: *           set INFO and continue
1.1       bertrand  469: *
                    470:             IF( INFO.EQ.0 )
                    471:      $         INFO = K
                    472:             KP = K
                    473:          ELSE
                    474:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
                    475: *
                    476: *              no interchange, use 1-by-1 pivot block
                    477: *
                    478:                KP = K
                    479:             ELSE
                    480: *
                    481: *              JMAX is the column-index of the largest off-diagonal
                    482: *              element in row IMAX, and ROWMAX is its absolute value
                    483: *
                    484:                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
                    485:                ROWMAX = CABS1( A( IMAX, JMAX ) )
                    486:                IF( IMAX.LT.N ) THEN
                    487:                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
                    488:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
                    489:                END IF
                    490: *
                    491:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
                    492: *
                    493: *                 no interchange, use 1-by-1 pivot block
                    494: *
                    495:                   KP = K
                    496:                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
                    497: *
                    498: *                 interchange rows and columns K and IMAX, use 1-by-1
                    499: *                 pivot block
                    500: *
                    501:                   KP = IMAX
                    502:                ELSE
                    503: *
                    504: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
                    505: *                 pivot block
                    506: *
                    507:                   KP = IMAX
                    508:                   KSTEP = 2
                    509:                END IF
                    510:             END IF
                    511: *
                    512:             KK = K + KSTEP - 1
                    513:             IF( KP.NE.KK ) THEN
                    514: *
                    515: *              Interchange rows and columns KK and KP in the trailing
                    516: *              submatrix A(k:n,k:n)
                    517: *
                    518:                IF( KP.LT.N )
                    519:      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
                    520:                CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
                    521:      $                     LDA )
                    522:                T = A( KK, KK )
                    523:                A( KK, KK ) = A( KP, KP )
                    524:                A( KP, KP ) = T
                    525:                IF( KSTEP.EQ.2 ) THEN
                    526:                   T = A( K+1, K )
                    527:                   A( K+1, K ) = A( KP, K )
                    528:                   A( KP, K ) = T
                    529:                END IF
                    530:             END IF
                    531: *
                    532: *           Update the trailing submatrix
                    533: *
                    534:             IF( KSTEP.EQ.1 ) THEN
                    535: *
                    536: *              1-by-1 pivot block D(k): column k now holds
                    537: *
                    538: *              W(k) = L(k)*D(k)
                    539: *
                    540: *              where L(k) is the k-th column of L
                    541: *
                    542:                IF( K.LT.N ) THEN
                    543: *
                    544: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
                    545: *
1.8       bertrand  546: *                 A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
1.1       bertrand  547: *
                    548:                   R1 = CONE / A( K, K )
                    549:                   CALL ZSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
                    550:      $                       A( K+1, K+1 ), LDA )
                    551: *
                    552: *                 Store L(k) in column K
                    553: *
                    554:                   CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
                    555:                END IF
                    556:             ELSE
                    557: *
                    558: *              2-by-2 pivot block D(k)
                    559: *
                    560:                IF( K.LT.N-1 ) THEN
                    561: *
                    562: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
                    563: *
1.8       bertrand  564: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
                    565: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
1.1       bertrand  566: *
                    567: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
                    568: *                 columns of L
                    569: *
                    570:                   D21 = A( K+1, K )
                    571:                   D11 = A( K+1, K+1 ) / D21
                    572:                   D22 = A( K, K ) / D21
                    573:                   T = CONE / ( D11*D22-CONE )
                    574:                   D21 = T / D21
                    575: *
                    576:                   DO 60 J = K + 2, N
                    577:                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
                    578:                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
                    579:                      DO 50 I = J, N
                    580:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
                    581:      $                              A( I, K+1 )*WKP1
                    582:    50                CONTINUE
                    583:                      A( J, K ) = WK
                    584:                      A( J, K+1 ) = WKP1
                    585:    60             CONTINUE
                    586:                END IF
                    587:             END IF
                    588:          END IF
                    589: *
                    590: *        Store details of the interchanges in IPIV
                    591: *
                    592:          IF( KSTEP.EQ.1 ) THEN
                    593:             IPIV( K ) = KP
                    594:          ELSE
                    595:             IPIV( K ) = -KP
                    596:             IPIV( K+1 ) = -KP
                    597:          END IF
                    598: *
                    599: *        Increase K and return to the start of the main loop
                    600: *
                    601:          K = K + KSTEP
                    602:          GO TO 40
                    603: *
                    604:       END IF
                    605: *
                    606:    70 CONTINUE
                    607:       RETURN
                    608: *
                    609: *     End of ZSYTF2
                    610: *
                    611:       END

CVSweb interface <joel.bertrand@systella.fr>