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