Return to dsytf2.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DSYTF2 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DSYTF2( 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: *       DOUBLE PRECISION   A( LDA, * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DSYTF2 computes the factorization of a real symmetric matrix A using
   39: *> 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 DOUBLE PRECISION 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.
   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.
  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: *
  130: *> \date November 2013
  131: *
  132: *> \ingroup doubleSYcomputational
  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.204 and l.372
  183: *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  184: *>    by
  185: *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  186: *>
  187: *>  01-01-96 - Based on modifications by
  188: *>    J. Lewis, Boeing Computer Services Company
  189: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  190: *>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
  191: *>         Company
  192: *> \endverbatim
  193: *
  194: *  =====================================================================
  195:       SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO )
  196: *
  197: *  -- LAPACK computational routine (version 3.5.0) --
  198: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  199: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  200: *     November 2013
  201: *
  202: *     .. Scalar Arguments ..
  203:       CHARACTER          UPLO
  204:       INTEGER            INFO, LDA, N
  205: *     ..
  206: *     .. Array Arguments ..
  207:       INTEGER            IPIV( * )
  208:       DOUBLE PRECISION   A( LDA, * )
  209: *     ..
  210: *
  211: *  =====================================================================
  212: *
  213: *     .. Parameters ..
  214:       DOUBLE PRECISION   ZERO, ONE
  215:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  216:       DOUBLE PRECISION   EIGHT, SEVTEN
  217:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
  218: *     ..
  219: *     .. Local Scalars ..
  220:       LOGICAL            UPPER
  221:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
  222:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
  223:      $                   ROWMAX, T, WK, WKM1, WKP1
  224: *     ..
  225: *     .. External Functions ..
  226:       LOGICAL            LSAME, DISNAN
  227:       INTEGER            IDAMAX
  228:       EXTERNAL           LSAME, IDAMAX, DISNAN
  229: *     ..
  230: *     .. External Subroutines ..
  231:       EXTERNAL           DSCAL, DSWAP, DSYR, XERBLA
  232: *     ..
  233: *     .. Intrinsic Functions ..
  234:       INTRINSIC          ABS, MAX, SQRT
  235: *     ..
  236: *     .. Executable Statements ..
  237: *
  238: *     Test the input parameters.
  239: *
  240:       INFO = 0
  241:       UPPER = LSAME( UPLO, 'U' )
  242:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  243:          INFO = -1
  244:       ELSE IF( N.LT.0 ) THEN
  245:          INFO = -2
  246:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  247:          INFO = -4
  248:       END IF
  249:       IF( INFO.NE.0 ) THEN
  250:          CALL XERBLA( 'DSYTF2', -INFO )
  251:          RETURN
  252:       END IF
  253: *
  254: *     Initialize ALPHA for use in choosing pivot block size.
  255: *
  256:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
  257: *
  258:       IF( UPPER ) THEN
  259: *
  260: *        Factorize A as U*D*U**T using the upper triangle of A
  261: *
  262: *        K is the main loop index, decreasing from N to 1 in steps of
  263: *        1 or 2
  264: *
  265:          K = N
  266:    10    CONTINUE
  267: *
  268: *        If K < 1, exit from loop
  269: *
  270:          IF( K.LT.1 )
  271:      $      GO TO 70
  272:          KSTEP = 1
  273: *
  274: *        Determine rows and columns to be interchanged and whether
  275: *        a 1-by-1 or 2-by-2 pivot block will be used
  276: *
  277:          ABSAKK = ABS( A( K, K ) )
  278: *
  279: *        IMAX is the row-index of the largest off-diagonal element in
  280: *        column K, and COLMAX is its absolute value.
  281: *        Determine both COLMAX and IMAX.
  282: *
  283:          IF( K.GT.1 ) THEN
  284:             IMAX = IDAMAX( K-1, A( 1, K ), 1 )
  285:             COLMAX = ABS( A( IMAX, K ) )
  286:          ELSE
  287:             COLMAX = ZERO
  288:          END IF
  289: *
  290:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  291: *
  292: *           Column K is zero or underflow, or contains a NaN:
  293: *           set INFO and continue
  294: *
  295:             IF( INFO.EQ.0 )
  296:      $         INFO = K
  297:             KP = K
  298:          ELSE
  299:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  300: *
  301: *              no interchange, use 1-by-1 pivot block
  302: *
  303:                KP = K
  304:             ELSE
  305: *
  306: *              JMAX is the column-index of the largest off-diagonal
  307: *              element in row IMAX, and ROWMAX is its absolute value
  308: *
  309:                JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
  310:                ROWMAX = ABS( A( IMAX, JMAX ) )
  311:                IF( IMAX.GT.1 ) THEN
  312:                   JMAX = IDAMAX( IMAX-1, A( 1, IMAX ), 1 )
  313:                   ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
  314:                END IF
  315: *
  316:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  317: *
  318: *                 no interchange, use 1-by-1 pivot block
  319: *
  320:                   KP = K
  321:                ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
  322: *
  323: *                 interchange rows and columns K and IMAX, use 1-by-1
  324: *                 pivot block
  325: *
  326:                   KP = IMAX
  327:                ELSE
  328: *
  329: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
  330: *                 pivot block
  331: *
  332:                   KP = IMAX
  333:                   KSTEP = 2
  334:                END IF
  335:             END IF
  336: *
  337:             KK = K - KSTEP + 1
  338:             IF( KP.NE.KK ) THEN
  339: *
  340: *              Interchange rows and columns KK and KP in the leading
  341: *              submatrix A(1:k,1:k)
  342: *
  343:                CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
  344:                CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
  345:      $                     LDA )
  346:                T = A( KK, KK )
  347:                A( KK, KK ) = A( KP, KP )
  348:                A( KP, KP ) = T
  349:                IF( KSTEP.EQ.2 ) THEN
  350:                   T = A( K-1, K )
  351:                   A( K-1, K ) = A( KP, K )
  352:                   A( KP, K ) = T
  353:                END IF
  354:             END IF
  355: *
  356: *           Update the leading submatrix
  357: *
  358:             IF( KSTEP.EQ.1 ) THEN
  359: *
  360: *              1-by-1 pivot block D(k): column k now holds
  361: *
  362: *              W(k) = U(k)*D(k)
  363: *
  364: *              where U(k) is the k-th column of U
  365: *
  366: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
  367: *
  368: *              A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
  369: *
  370:                R1 = ONE / A( K, K )
  371:                CALL DSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
  372: *
  373: *              Store U(k) in column k
  374: *
  375:                CALL DSCAL( K-1, R1, A( 1, K ), 1 )
  376:             ELSE
  377: *
  378: *              2-by-2 pivot block D(k): columns k and k-1 now hold
  379: *
  380: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
  381: *
  382: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
  383: *              of U
  384: *
  385: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
  386: *
  387: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
  388: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
  389: *
  390:                IF( K.GT.2 ) THEN
  391: *
  392:                   D12 = A( K-1, K )
  393:                   D22 = A( K-1, K-1 ) / D12
  394:                   D11 = A( K, K ) / D12
  395:                   T = ONE / ( D11*D22-ONE )
  396:                   D12 = T / D12
  397: *
  398:                   DO 30 J = K - 2, 1, -1
  399:                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
  400:                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
  401:                      DO 20 I = J, 1, -1
  402:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
  403:      $                              A( I, K-1 )*WKM1
  404:    20                CONTINUE
  405:                      A( J, K ) = WK
  406:                      A( J, K-1 ) = WKM1
  407:    30             CONTINUE
  408: *
  409:                END IF
  410: *
  411:             END IF
  412:          END IF
  413: *
  414: *        Store details of the interchanges in IPIV
  415: *
  416:          IF( KSTEP.EQ.1 ) THEN
  417:             IPIV( K ) = KP
  418:          ELSE
  419:             IPIV( K ) = -KP
  420:             IPIV( K-1 ) = -KP
  421:          END IF
  422: *
  423: *        Decrease K and return to the start of the main loop
  424: *
  425:          K = K - KSTEP
  426:          GO TO 10
  427: *
  428:       ELSE
  429: *
  430: *        Factorize A as L*D*L**T using the lower triangle of A
  431: *
  432: *        K is the main loop index, increasing from 1 to N in steps of
  433: *        1 or 2
  434: *
  435:          K = 1
  436:    40    CONTINUE
  437: *
  438: *        If K > N, exit from loop
  439: *
  440:          IF( K.GT.N )
  441:      $      GO TO 70
  442:          KSTEP = 1
  443: *
  444: *        Determine rows and columns to be interchanged and whether
  445: *        a 1-by-1 or 2-by-2 pivot block will be used
  446: *
  447:          ABSAKK = ABS( A( K, K ) )
  448: *
  449: *        IMAX is the row-index of the largest off-diagonal element in
  450: *        column K, and COLMAX is its absolute value.
  451: *        Determine both COLMAX and IMAX.
  452: *
  453:          IF( K.LT.N ) THEN
  454:             IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 )
  455:             COLMAX = ABS( A( IMAX, K ) )
  456:          ELSE
  457:             COLMAX = ZERO
  458:          END IF
  459: *
  460:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  461: *
  462: *           Column K is zero or underflow, or contains a NaN:
  463: *           set INFO and continue
  464: *
  465:             IF( INFO.EQ.0 )
  466:      $         INFO = K
  467:             KP = K
  468:          ELSE
  469:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  470: *
  471: *              no interchange, use 1-by-1 pivot block
  472: *
  473:                KP = K
  474:             ELSE
  475: *
  476: *              JMAX is the column-index of the largest off-diagonal
  477: *              element in row IMAX, and ROWMAX is its absolute value
  478: *
  479:                JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA )
  480:                ROWMAX = ABS( A( IMAX, JMAX ) )
  481:                IF( IMAX.LT.N ) THEN
  482:                   JMAX = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
  483:                   ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
  484:                END IF
  485: *
  486:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  487: *
  488: *                 no interchange, use 1-by-1 pivot block
  489: *
  490:                   KP = K
  491:                ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
  492: *
  493: *                 interchange rows and columns K and IMAX, use 1-by-1
  494: *                 pivot block
  495: *
  496:                   KP = IMAX
  497:                ELSE
  498: *
  499: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
  500: *                 pivot block
  501: *
  502:                   KP = IMAX
  503:                   KSTEP = 2
  504:                END IF
  505:             END IF
  506: *
  507:             KK = K + KSTEP - 1
  508:             IF( KP.NE.KK ) THEN
  509: *
  510: *              Interchange rows and columns KK and KP in the trailing
  511: *              submatrix A(k:n,k:n)
  512: *
  513:                IF( KP.LT.N )
  514:      $            CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
  515:                CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
  516:      $                     LDA )
  517:                T = A( KK, KK )
  518:                A( KK, KK ) = A( KP, KP )
  519:                A( KP, KP ) = T
  520:                IF( KSTEP.EQ.2 ) THEN
  521:                   T = A( K+1, K )
  522:                   A( K+1, K ) = A( KP, K )
  523:                   A( KP, K ) = T
  524:                END IF
  525:             END IF
  526: *
  527: *           Update the trailing submatrix
  528: *
  529:             IF( KSTEP.EQ.1 ) THEN
  530: *
  531: *              1-by-1 pivot block D(k): column k now holds
  532: *
  533: *              W(k) = L(k)*D(k)
  534: *
  535: *              where L(k) is the k-th column of L
  536: *
  537:                IF( K.LT.N ) THEN
  538: *
  539: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
  540: *
  541: *                 A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
  542: *
  543:                   D11 = ONE / A( K, K )
  544:                   CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
  545:      $                       A( K+1, K+1 ), LDA )
  546: *
  547: *                 Store L(k) in column K
  548: *
  549:                   CALL DSCAL( N-K, D11, A( K+1, K ), 1 )
  550:                END IF
  551:             ELSE
  552: *
  553: *              2-by-2 pivot block D(k)
  554: *
  555:                IF( K.LT.N-1 ) THEN
  556: *
  557: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
  558: *
  559: *                 A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))**T
  560: *
  561: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
  562: *                 columns of L
  563: *
  564:                   D21 = A( K+1, K )
  565:                   D11 = A( K+1, K+1 ) / D21
  566:                   D22 = A( K, K ) / D21
  567:                   T = ONE / ( D11*D22-ONE )
  568:                   D21 = T / D21
  569: *
  570:                   DO 60 J = K + 2, N
  571: *
  572:                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
  573:                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
  574: *
  575:                      DO 50 I = J, N
  576:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
  577:      $                              A( I, K+1 )*WKP1
  578:    50                CONTINUE
  579: *
  580:                      A( J, K ) = WK
  581:                      A( J, K+1 ) = WKP1
  582: *
  583:    60             CONTINUE
  584:                END IF
  585:             END IF
  586:          END IF
  587: *
  588: *        Store details of the interchanges in IPIV
  589: *
  590:          IF( KSTEP.EQ.1 ) THEN
  591:             IPIV( K ) = KP
  592:          ELSE
  593:             IPIV( K ) = -KP
  594:             IPIV( K+1 ) = -KP
  595:          END IF
  596: *
  597: *        Increase K and return to the start of the main loop
  598: *
  599:          K = K + KSTEP
  600:          GO TO 40
  601: *
  602:       END IF
  603: *
  604:    70 CONTINUE
  605: *
  606:       RETURN
  607: *
  608: *     End of DSYTF2
  609: *
  610:       END