File:  [local] / rpl / lapack / lapack / zhptrf.f
Revision 1.5: download - view: text, annotated - select for diffs - revision graph
Sat Aug 7 13:22:35 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       CHARACTER          UPLO
   10:       INTEGER            INFO, N
   11: *     ..
   12: *     .. Array Arguments ..
   13:       INTEGER            IPIV( * )
   14:       COMPLEX*16         AP( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  ZHPTRF computes the factorization of a complex Hermitian packed
   21: *  matrix A using the Bunch-Kaufman diagonal pivoting method:
   22: *
   23: *     A = U*D*U**H  or  A = L*D*L**H
   24: *
   25: *  where U (or L) is a product of permutation and unit upper (lower)
   26: *  triangular matrices, and D is Hermitian 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) COMPLEX*16 array, dimension (N*(N+1)/2)
   40: *          On entry, the upper or lower triangle of the Hermitian 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: *
   74: *  If UPLO = 'U', then A = U*D*U', where
   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: *
   91: *  If UPLO = 'L', then A = L*D*L', where
   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, D, D11, D22, R1, ROWMAX,
  121:      $                   TT
  122:       COMPLEX*16         D12, D21, T, WK, WKM1, WKP1, ZDUM
  123: *     ..
  124: *     .. External Functions ..
  125:       LOGICAL            LSAME
  126:       INTEGER            IZAMAX
  127:       DOUBLE PRECISION   DLAPY2
  128:       EXTERNAL           LSAME, IZAMAX, DLAPY2
  129: *     ..
  130: *     .. External Subroutines ..
  131:       EXTERNAL           XERBLA, ZDSCAL, ZHPR, ZSWAP
  132: *     ..
  133: *     .. Intrinsic Functions ..
  134:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
  135: *     ..
  136: *     .. Statement Functions ..
  137:       DOUBLE PRECISION   CABS1
  138: *     ..
  139: *     .. Statement Function definitions ..
  140:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  141: *     ..
  142: *     .. Executable Statements ..
  143: *
  144: *     Test the input parameters.
  145: *
  146:       INFO = 0
  147:       UPPER = LSAME( UPLO, 'U' )
  148:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  149:          INFO = -1
  150:       ELSE IF( N.LT.0 ) THEN
  151:          INFO = -2
  152:       END IF
  153:       IF( INFO.NE.0 ) THEN
  154:          CALL XERBLA( 'ZHPTRF', -INFO )
  155:          RETURN
  156:       END IF
  157: *
  158: *     Initialize ALPHA for use in choosing pivot block size.
  159: *
  160:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
  161: *
  162:       IF( UPPER ) THEN
  163: *
  164: *        Factorize A as U*D*U' using the upper triangle of A
  165: *
  166: *        K is the main loop index, decreasing from N to 1 in steps of
  167: *        1 or 2
  168: *
  169:          K = N
  170:          KC = ( N-1 )*N / 2 + 1
  171:    10    CONTINUE
  172:          KNC = KC
  173: *
  174: *        If K < 1, exit from loop
  175: *
  176:          IF( K.LT.1 )
  177:      $      GO TO 110
  178:          KSTEP = 1
  179: *
  180: *        Determine rows and columns to be interchanged and whether
  181: *        a 1-by-1 or 2-by-2 pivot block will be used
  182: *
  183:          ABSAKK = ABS( DBLE( AP( KC+K-1 ) ) )
  184: *
  185: *        IMAX is the row-index of the largest off-diagonal element in
  186: *        column K, and COLMAX is its absolute value
  187: *
  188:          IF( K.GT.1 ) THEN
  189:             IMAX = IZAMAX( K-1, AP( KC ), 1 )
  190:             COLMAX = CABS1( AP( KC+IMAX-1 ) )
  191:          ELSE
  192:             COLMAX = ZERO
  193:          END IF
  194: *
  195:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  196: *
  197: *           Column K is zero: set INFO and continue
  198: *
  199:             IF( INFO.EQ.0 )
  200:      $         INFO = K
  201:             KP = K
  202:             AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
  203:          ELSE
  204:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  205: *
  206: *              no interchange, use 1-by-1 pivot block
  207: *
  208:                KP = K
  209:             ELSE
  210: *
  211: *              JMAX is the column-index of the largest off-diagonal
  212: *              element in row IMAX, and ROWMAX is its absolute value
  213: *
  214:                ROWMAX = ZERO
  215:                JMAX = IMAX
  216:                KX = IMAX*( IMAX+1 ) / 2 + IMAX
  217:                DO 20 J = IMAX + 1, K
  218:                   IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
  219:                      ROWMAX = CABS1( AP( KX ) )
  220:                      JMAX = J
  221:                   END IF
  222:                   KX = KX + J
  223:    20          CONTINUE
  224:                KPC = ( IMAX-1 )*IMAX / 2 + 1
  225:                IF( IMAX.GT.1 ) THEN
  226:                   JMAX = IZAMAX( IMAX-1, AP( KPC ), 1 )
  227:                   ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-1 ) ) )
  228:                END IF
  229: *
  230:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  231: *
  232: *                 no interchange, use 1-by-1 pivot block
  233: *
  234:                   KP = K
  235:                ELSE IF( ABS( DBLE( AP( KPC+IMAX-1 ) ) ).GE.ALPHA*
  236:      $                  ROWMAX ) THEN
  237: *
  238: *                 interchange rows and columns K and IMAX, use 1-by-1
  239: *                 pivot block
  240: *
  241:                   KP = IMAX
  242:                ELSE
  243: *
  244: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
  245: *                 pivot block
  246: *
  247:                   KP = IMAX
  248:                   KSTEP = 2
  249:                END IF
  250:             END IF
  251: *
  252:             KK = K - KSTEP + 1
  253:             IF( KSTEP.EQ.2 )
  254:      $         KNC = KNC - K + 1
  255:             IF( KP.NE.KK ) THEN
  256: *
  257: *              Interchange rows and columns KK and KP in the leading
  258: *              submatrix A(1:k,1:k)
  259: *
  260:                CALL ZSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
  261:                KX = KPC + KP - 1
  262:                DO 30 J = KP + 1, KK - 1
  263:                   KX = KX + J - 1
  264:                   T = DCONJG( AP( KNC+J-1 ) )
  265:                   AP( KNC+J-1 ) = DCONJG( AP( KX ) )
  266:                   AP( KX ) = T
  267:    30          CONTINUE
  268:                AP( KX+KK-1 ) = DCONJG( AP( KX+KK-1 ) )
  269:                R1 = DBLE( AP( KNC+KK-1 ) )
  270:                AP( KNC+KK-1 ) = DBLE( AP( KPC+KP-1 ) )
  271:                AP( KPC+KP-1 ) = R1
  272:                IF( KSTEP.EQ.2 ) THEN
  273:                   AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
  274:                   T = AP( KC+K-2 )
  275:                   AP( KC+K-2 ) = AP( KC+KP-1 )
  276:                   AP( KC+KP-1 ) = T
  277:                END IF
  278:             ELSE
  279:                AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
  280:                IF( KSTEP.EQ.2 )
  281:      $            AP( KC-1 ) = DBLE( AP( KC-1 ) )
  282:             END IF
  283: *
  284: *           Update the leading submatrix
  285: *
  286:             IF( KSTEP.EQ.1 ) THEN
  287: *
  288: *              1-by-1 pivot block D(k): column k now holds
  289: *
  290: *              W(k) = U(k)*D(k)
  291: *
  292: *              where U(k) is the k-th column of U
  293: *
  294: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
  295: *
  296: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
  297: *
  298:                R1 = ONE / DBLE( AP( KC+K-1 ) )
  299:                CALL ZHPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
  300: *
  301: *              Store U(k) in column k
  302: *
  303:                CALL ZDSCAL( K-1, R1, AP( KC ), 1 )
  304:             ELSE
  305: *
  306: *              2-by-2 pivot block D(k): columns k and k-1 now hold
  307: *
  308: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
  309: *
  310: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
  311: *              of U
  312: *
  313: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
  314: *
  315: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
  316: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
  317: *
  318:                IF( K.GT.2 ) THEN
  319: *
  320:                   D = DLAPY2( DBLE( AP( K-1+( K-1 )*K / 2 ) ),
  321:      $                DIMAG( AP( K-1+( K-1 )*K / 2 ) ) )
  322:                   D22 = DBLE( AP( K-1+( K-2 )*( K-1 ) / 2 ) ) / D
  323:                   D11 = DBLE( AP( K+( K-1 )*K / 2 ) ) / D
  324:                   TT = ONE / ( D11*D22-ONE )
  325:                   D12 = AP( K-1+( K-1 )*K / 2 ) / D
  326:                   D = TT / D
  327: *
  328:                   DO 50 J = K - 2, 1, -1
  329:                      WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
  330:      $                      DCONJG( D12 )*AP( J+( K-1 )*K / 2 ) )
  331:                      WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12*
  332:      $                    AP( J+( K-2 )*( K-1 ) / 2 ) )
  333:                      DO 40 I = J, 1, -1
  334:                         AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
  335:      $                     AP( I+( K-1 )*K / 2 )*DCONJG( WK ) -
  336:      $                     AP( I+( K-2 )*( K-1 ) / 2 )*DCONJG( WKM1 )
  337:    40                CONTINUE
  338:                      AP( J+( K-1 )*K / 2 ) = WK
  339:                      AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
  340:                      AP( J+( J-1 )*J / 2 ) = DCMPLX( DBLE( AP( J+( J-
  341:      $                                       1 )*J / 2 ) ), 0.0D+0 )
  342:    50             CONTINUE
  343: *
  344:                END IF
  345: *
  346:             END IF
  347:          END IF
  348: *
  349: *        Store details of the interchanges in IPIV
  350: *
  351:          IF( KSTEP.EQ.1 ) THEN
  352:             IPIV( K ) = KP
  353:          ELSE
  354:             IPIV( K ) = -KP
  355:             IPIV( K-1 ) = -KP
  356:          END IF
  357: *
  358: *        Decrease K and return to the start of the main loop
  359: *
  360:          K = K - KSTEP
  361:          KC = KNC - K
  362:          GO TO 10
  363: *
  364:       ELSE
  365: *
  366: *        Factorize A as L*D*L' using the lower triangle of A
  367: *
  368: *        K is the main loop index, increasing from 1 to N in steps of
  369: *        1 or 2
  370: *
  371:          K = 1
  372:          KC = 1
  373:          NPP = N*( N+1 ) / 2
  374:    60    CONTINUE
  375:          KNC = KC
  376: *
  377: *        If K > N, exit from loop
  378: *
  379:          IF( K.GT.N )
  380:      $      GO TO 110
  381:          KSTEP = 1
  382: *
  383: *        Determine rows and columns to be interchanged and whether
  384: *        a 1-by-1 or 2-by-2 pivot block will be used
  385: *
  386:          ABSAKK = ABS( DBLE( AP( KC ) ) )
  387: *
  388: *        IMAX is the row-index of the largest off-diagonal element in
  389: *        column K, and COLMAX is its absolute value
  390: *
  391:          IF( K.LT.N ) THEN
  392:             IMAX = K + IZAMAX( N-K, AP( KC+1 ), 1 )
  393:             COLMAX = CABS1( AP( KC+IMAX-K ) )
  394:          ELSE
  395:             COLMAX = ZERO
  396:          END IF
  397: *
  398:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  399: *
  400: *           Column K is zero: set INFO and continue
  401: *
  402:             IF( INFO.EQ.0 )
  403:      $         INFO = K
  404:             KP = K
  405:             AP( KC ) = DBLE( AP( KC ) )
  406:          ELSE
  407:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  408: *
  409: *              no interchange, use 1-by-1 pivot block
  410: *
  411:                KP = K
  412:             ELSE
  413: *
  414: *              JMAX is the column-index of the largest off-diagonal
  415: *              element in row IMAX, and ROWMAX is its absolute value
  416: *
  417:                ROWMAX = ZERO
  418:                KX = KC + IMAX - K
  419:                DO 70 J = K, IMAX - 1
  420:                   IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
  421:                      ROWMAX = CABS1( AP( KX ) )
  422:                      JMAX = J
  423:                   END IF
  424:                   KX = KX + N - J
  425:    70          CONTINUE
  426:                KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
  427:                IF( IMAX.LT.N ) THEN
  428:                   JMAX = IMAX + IZAMAX( N-IMAX, AP( KPC+1 ), 1 )
  429:                   ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-IMAX ) ) )
  430:                END IF
  431: *
  432:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  433: *
  434: *                 no interchange, use 1-by-1 pivot block
  435: *
  436:                   KP = K
  437:                ELSE IF( ABS( DBLE( AP( KPC ) ) ).GE.ALPHA*ROWMAX ) THEN
  438: *
  439: *                 interchange rows and columns K and IMAX, use 1-by-1
  440: *                 pivot block
  441: *
  442:                   KP = IMAX
  443:                ELSE
  444: *
  445: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
  446: *                 pivot block
  447: *
  448:                   KP = IMAX
  449:                   KSTEP = 2
  450:                END IF
  451:             END IF
  452: *
  453:             KK = K + KSTEP - 1
  454:             IF( KSTEP.EQ.2 )
  455:      $         KNC = KNC + N - K + 1
  456:             IF( KP.NE.KK ) THEN
  457: *
  458: *              Interchange rows and columns KK and KP in the trailing
  459: *              submatrix A(k:n,k:n)
  460: *
  461:                IF( KP.LT.N )
  462:      $            CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
  463:      $                        1 )
  464:                KX = KNC + KP - KK
  465:                DO 80 J = KK + 1, KP - 1
  466:                   KX = KX + N - J + 1
  467:                   T = DCONJG( AP( KNC+J-KK ) )
  468:                   AP( KNC+J-KK ) = DCONJG( AP( KX ) )
  469:                   AP( KX ) = T
  470:    80          CONTINUE
  471:                AP( KNC+KP-KK ) = DCONJG( AP( KNC+KP-KK ) )
  472:                R1 = DBLE( AP( KNC ) )
  473:                AP( KNC ) = DBLE( AP( KPC ) )
  474:                AP( KPC ) = R1
  475:                IF( KSTEP.EQ.2 ) THEN
  476:                   AP( KC ) = DBLE( AP( KC ) )
  477:                   T = AP( KC+1 )
  478:                   AP( KC+1 ) = AP( KC+KP-K )
  479:                   AP( KC+KP-K ) = T
  480:                END IF
  481:             ELSE
  482:                AP( KC ) = DBLE( AP( KC ) )
  483:                IF( KSTEP.EQ.2 )
  484:      $            AP( KNC ) = DBLE( AP( KNC ) )
  485:             END IF
  486: *
  487: *           Update the trailing submatrix
  488: *
  489:             IF( KSTEP.EQ.1 ) THEN
  490: *
  491: *              1-by-1 pivot block D(k): column k now holds
  492: *
  493: *              W(k) = L(k)*D(k)
  494: *
  495: *              where L(k) is the k-th column of L
  496: *
  497:                IF( K.LT.N ) THEN
  498: *
  499: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
  500: *
  501: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
  502: *
  503:                   R1 = ONE / DBLE( AP( KC ) )
  504:                   CALL ZHPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
  505:      $                       AP( KC+N-K+1 ) )
  506: *
  507: *                 Store L(k) in column K
  508: *
  509:                   CALL ZDSCAL( N-K, R1, AP( KC+1 ), 1 )
  510:                END IF
  511:             ELSE
  512: *
  513: *              2-by-2 pivot block D(k): columns K and K+1 now hold
  514: *
  515: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
  516: *
  517: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
  518: *              of L
  519: *
  520:                IF( K.LT.N-1 ) THEN
  521: *
  522: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
  523: *
  524: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
  525: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
  526: *
  527: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
  528: *                 columns of L
  529: *
  530:                   D = DLAPY2( DBLE( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ),
  531:      $                DIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) )
  532:                   D11 = DBLE( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D
  533:                   D22 = DBLE( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D
  534:                   TT = ONE / ( D11*D22-ONE )
  535:                   D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D
  536:                   D = TT / D
  537: *
  538:                   DO 100 J = K + 2, N
  539:                      WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21*
  540:      $                    AP( J+K*( 2*N-K-1 ) / 2 ) )
  541:                      WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
  542:      $                      DCONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) /
  543:      $                      2 ) )
  544:                      DO 90 I = J, N
  545:                         AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
  546:      $                     ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
  547:      $                     2 )*DCONJG( WK ) - AP( I+K*( 2*N-K-1 ) / 2 )*
  548:      $                     DCONJG( WKP1 )
  549:    90                CONTINUE
  550:                      AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
  551:                      AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
  552:                      AP( J+( J-1 )*( 2*N-J ) / 2 )
  553:      $                  = DCMPLX( DBLE( AP( J+( J-1 )*( 2*N-J ) / 2 ) ),
  554:      $                  0.0D+0 )
  555:   100             CONTINUE
  556:                END IF
  557:             END IF
  558:          END IF
  559: *
  560: *        Store details of the interchanges in IPIV
  561: *
  562:          IF( KSTEP.EQ.1 ) THEN
  563:             IPIV( K ) = KP
  564:          ELSE
  565:             IPIV( K ) = -KP
  566:             IPIV( K+1 ) = -KP
  567:          END IF
  568: *
  569: *        Increase K and return to the start of the main loop
  570: *
  571:          K = K + KSTEP
  572:          KC = KNC + N - K + 2
  573:          GO TO 60
  574: *
  575:       END IF
  576: *
  577:   110 CONTINUE
  578:       RETURN
  579: *
  580: *     End of ZHPTRF
  581: *
  582:       END

CVSweb interface <joel.bertrand@systella.fr>