Annotation of rpl/lapack/lapack/dsptrf.f, revision 1.8

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

CVSweb interface <joel.bertrand@systella.fr>