Annotation of rpl/lapack/lapack/zhptrf.f, revision 1.1
1.1 ! bertrand 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>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhptrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack