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*> \ingroup complex16OTHERcomputational 111* 112*> \par Further Details: 113* ===================== 114*> 115*> \verbatim 116*> 117*> If UPLO = 'U', then A = U*D*U**H, where 118*> U = P(n)*U(n)* ... *P(k)U(k)* ..., 119*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to 120*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 121*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 122*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 123*> that if the diagonal block D(k) is of order s (s = 1 or 2), then 124*> 125*> ( I v 0 ) k-s 126*> U(k) = ( 0 I 0 ) s 127*> ( 0 0 I ) n-k 128*> k-s s n-k 129*> 130*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 131*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 132*> and A(k,k), and v overwrites A(1:k-2,k-1:k). 133*> 134*> If UPLO = 'L', then A = L*D*L**H, where 135*> L = P(1)*L(1)* ... *P(k)*L(k)* ..., 136*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 137*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 138*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 139*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 140*> that if the diagonal block D(k) is of order s (s = 1 or 2), then 141*> 142*> ( I 0 0 ) k-1 143*> L(k) = ( 0 I 0 ) s 144*> ( 0 v I ) n-k-s+1 145*> k-1 s n-k-s+1 146*> 147*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 148*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 149*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 150*> \endverbatim 151* 152*> \par Contributors: 153* ================== 154*> 155*> J. Lewis, Boeing Computer Services Company 156* 157* ===================================================================== 158 SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) 159* 160* -- LAPACK computational routine -- 161* -- LAPACK is a software package provided by Univ. of Tennessee, -- 162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 163* 164* .. Scalar Arguments .. 165 CHARACTER UPLO 166 INTEGER INFO, N 167* .. 168* .. Array Arguments .. 169 INTEGER IPIV( * ) 170 COMPLEX*16 AP( * ) 171* .. 172* 173* ===================================================================== 174* 175* .. Parameters .. 176 DOUBLE PRECISION ZERO, ONE 177 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 178 DOUBLE PRECISION EIGHT, SEVTEN 179 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 180* .. 181* .. Local Scalars .. 182 LOGICAL UPPER 183 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC, 184 $ KSTEP, KX, NPP 185 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX, 186 $ TT 187 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM 188* .. 189* .. External Functions .. 190 LOGICAL LSAME 191 INTEGER IZAMAX 192 DOUBLE PRECISION DLAPY2 193 EXTERNAL LSAME, IZAMAX, DLAPY2 194* .. 195* .. External Subroutines .. 196 EXTERNAL XERBLA, ZDSCAL, ZHPR, ZSWAP 197* .. 198* .. Intrinsic Functions .. 199 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT 200* .. 201* .. Statement Functions .. 202 DOUBLE PRECISION CABS1 203* .. 204* .. Statement Function definitions .. 205 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 206* .. 207* .. Executable Statements .. 208* 209* Test the input parameters. 210* 211 INFO = 0 212 UPPER = LSAME( UPLO, 'U' ) 213 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 214 INFO = -1 215 ELSE IF( N.LT.0 ) THEN 216 INFO = -2 217 END IF 218 IF( INFO.NE.0 ) THEN 219 CALL XERBLA( 'ZHPTRF', -INFO ) 220 RETURN 221 END IF 222* 223* Initialize ALPHA for use in choosing pivot block size. 224* 225 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 226* 227 IF( UPPER ) THEN 228* 229* Factorize A as U*D*U**H using the upper triangle of A 230* 231* K is the main loop index, decreasing from N to 1 in steps of 232* 1 or 2 233* 234 K = N 235 KC = ( N-1 )*N / 2 + 1 236 10 CONTINUE 237 KNC = KC 238* 239* If K < 1, exit from loop 240* 241 IF( K.LT.1 ) 242 $ GO TO 110 243 KSTEP = 1 244* 245* Determine rows and columns to be interchanged and whether 246* a 1-by-1 or 2-by-2 pivot block will be used 247* 248 ABSAKK = ABS( DBLE( AP( KC+K-1 ) ) ) 249* 250* IMAX is the row-index of the largest off-diagonal element in 251* column K, and COLMAX is its absolute value 252* 253 IF( K.GT.1 ) THEN 254 IMAX = IZAMAX( K-1, AP( KC ), 1 ) 255 COLMAX = CABS1( AP( KC+IMAX-1 ) ) 256 ELSE 257 COLMAX = ZERO 258 END IF 259* 260 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 261* 262* Column K is zero: set INFO and continue 263* 264 IF( INFO.EQ.0 ) 265 $ INFO = K 266 KP = K 267 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) ) 268 ELSE 269 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 270* 271* no interchange, use 1-by-1 pivot block 272* 273 KP = K 274 ELSE 275* 276* JMAX is the column-index of the largest off-diagonal 277* element in row IMAX, and ROWMAX is its absolute value 278* 279 ROWMAX = ZERO 280 JMAX = IMAX 281 KX = IMAX*( IMAX+1 ) / 2 + IMAX 282 DO 20 J = IMAX + 1, K 283 IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN 284 ROWMAX = CABS1( AP( KX ) ) 285 JMAX = J 286 END IF 287 KX = KX + J 288 20 CONTINUE 289 KPC = ( IMAX-1 )*IMAX / 2 + 1 290 IF( IMAX.GT.1 ) THEN 291 JMAX = IZAMAX( IMAX-1, AP( KPC ), 1 ) 292 ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-1 ) ) ) 293 END IF 294* 295 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 296* 297* no interchange, use 1-by-1 pivot block 298* 299 KP = K 300 ELSE IF( ABS( DBLE( AP( KPC+IMAX-1 ) ) ).GE.ALPHA* 301 $ ROWMAX ) THEN 302* 303* interchange rows and columns K and IMAX, use 1-by-1 304* pivot block 305* 306 KP = IMAX 307 ELSE 308* 309* interchange rows and columns K-1 and IMAX, use 2-by-2 310* pivot block 311* 312 KP = IMAX 313 KSTEP = 2 314 END IF 315 END IF 316* 317 KK = K - KSTEP + 1 318 IF( KSTEP.EQ.2 ) 319 $ KNC = KNC - K + 1 320 IF( KP.NE.KK ) THEN 321* 322* Interchange rows and columns KK and KP in the leading 323* submatrix A(1:k,1:k) 324* 325 CALL ZSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 ) 326 KX = KPC + KP - 1 327 DO 30 J = KP + 1, KK - 1 328 KX = KX + J - 1 329 T = DCONJG( AP( KNC+J-1 ) ) 330 AP( KNC+J-1 ) = DCONJG( AP( KX ) ) 331 AP( KX ) = T 332 30 CONTINUE 333 AP( KX+KK-1 ) = DCONJG( AP( KX+KK-1 ) ) 334 R1 = DBLE( AP( KNC+KK-1 ) ) 335 AP( KNC+KK-1 ) = DBLE( AP( KPC+KP-1 ) ) 336 AP( KPC+KP-1 ) = R1 337 IF( KSTEP.EQ.2 ) THEN 338 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) ) 339 T = AP( KC+K-2 ) 340 AP( KC+K-2 ) = AP( KC+KP-1 ) 341 AP( KC+KP-1 ) = T 342 END IF 343 ELSE 344 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) ) 345 IF( KSTEP.EQ.2 ) 346 $ AP( KC-1 ) = DBLE( AP( KC-1 ) ) 347 END IF 348* 349* Update the leading submatrix 350* 351 IF( KSTEP.EQ.1 ) THEN 352* 353* 1-by-1 pivot block D(k): column k now holds 354* 355* W(k) = U(k)*D(k) 356* 357* where U(k) is the k-th column of U 358* 359* Perform a rank-1 update of A(1:k-1,1:k-1) as 360* 361* A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H 362* 363 R1 = ONE / DBLE( AP( KC+K-1 ) ) 364 CALL ZHPR( UPLO, K-1, -R1, AP( KC ), 1, AP ) 365* 366* Store U(k) in column k 367* 368 CALL ZDSCAL( K-1, R1, AP( KC ), 1 ) 369 ELSE 370* 371* 2-by-2 pivot block D(k): columns k and k-1 now hold 372* 373* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 374* 375* where U(k) and U(k-1) are the k-th and (k-1)-th columns 376* of U 377* 378* Perform a rank-2 update of A(1:k-2,1:k-2) as 379* 380* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H 381* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H 382* 383 IF( K.GT.2 ) THEN 384* 385 D = DLAPY2( DBLE( AP( K-1+( K-1 )*K / 2 ) ), 386 $ DIMAG( AP( K-1+( K-1 )*K / 2 ) ) ) 387 D22 = DBLE( AP( K-1+( K-2 )*( K-1 ) / 2 ) ) / D 388 D11 = DBLE( AP( K+( K-1 )*K / 2 ) ) / D 389 TT = ONE / ( D11*D22-ONE ) 390 D12 = AP( K-1+( K-1 )*K / 2 ) / D 391 D = TT / D 392* 393 DO 50 J = K - 2, 1, -1 394 WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- 395 $ DCONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) 396 WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12* 397 $ AP( J+( K-2 )*( K-1 ) / 2 ) ) 398 DO 40 I = J, 1, -1 399 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - 400 $ AP( I+( K-1 )*K / 2 )*DCONJG( WK ) - 401 $ AP( I+( K-2 )*( K-1 ) / 2 )*DCONJG( WKM1 ) 402 40 CONTINUE 403 AP( J+( K-1 )*K / 2 ) = WK 404 AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1 405 AP( J+( J-1 )*J / 2 ) = DCMPLX( DBLE( AP( J+( J- 406 $ 1 )*J / 2 ) ), 0.0D+0 ) 407 50 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 KC = KNC - K 427 GO TO 10 428* 429 ELSE 430* 431* Factorize A as L*D*L**H using the lower triangle of A 432* 433* K is the main loop index, increasing from 1 to N in steps of 434* 1 or 2 435* 436 K = 1 437 KC = 1 438 NPP = N*( N+1 ) / 2 439 60 CONTINUE 440 KNC = KC 441* 442* If K > N, exit from loop 443* 444 IF( K.GT.N ) 445 $ GO TO 110 446 KSTEP = 1 447* 448* Determine rows and columns to be interchanged and whether 449* a 1-by-1 or 2-by-2 pivot block will be used 450* 451 ABSAKK = ABS( DBLE( AP( KC ) ) ) 452* 453* IMAX is the row-index of the largest off-diagonal element in 454* column K, and COLMAX is its absolute value 455* 456 IF( K.LT.N ) THEN 457 IMAX = K + IZAMAX( N-K, AP( KC+1 ), 1 ) 458 COLMAX = CABS1( AP( KC+IMAX-K ) ) 459 ELSE 460 COLMAX = ZERO 461 END IF 462* 463 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 464* 465* Column K is zero: set INFO and continue 466* 467 IF( INFO.EQ.0 ) 468 $ INFO = K 469 KP = K 470 AP( KC ) = DBLE( AP( KC ) ) 471 ELSE 472 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 473* 474* no interchange, use 1-by-1 pivot block 475* 476 KP = K 477 ELSE 478* 479* JMAX is the column-index of the largest off-diagonal 480* element in row IMAX, and ROWMAX is its absolute value 481* 482 ROWMAX = ZERO 483 KX = KC + IMAX - K 484 DO 70 J = K, IMAX - 1 485 IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN 486 ROWMAX = CABS1( AP( KX ) ) 487 JMAX = J 488 END IF 489 KX = KX + N - J 490 70 CONTINUE 491 KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1 492 IF( IMAX.LT.N ) THEN 493 JMAX = IMAX + IZAMAX( N-IMAX, AP( KPC+1 ), 1 ) 494 ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-IMAX ) ) ) 495 END IF 496* 497 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 498* 499* no interchange, use 1-by-1 pivot block 500* 501 KP = K 502 ELSE IF( ABS( DBLE( AP( KPC ) ) ).GE.ALPHA*ROWMAX ) THEN 503* 504* interchange rows and columns K and IMAX, use 1-by-1 505* pivot block 506* 507 KP = IMAX 508 ELSE 509* 510* interchange rows and columns K+1 and IMAX, use 2-by-2 511* pivot block 512* 513 KP = IMAX 514 KSTEP = 2 515 END IF 516 END IF 517* 518 KK = K + KSTEP - 1 519 IF( KSTEP.EQ.2 ) 520 $ KNC = KNC + N - K + 1 521 IF( KP.NE.KK ) THEN 522* 523* Interchange rows and columns KK and KP in the trailing 524* submatrix A(k:n,k:n) 525* 526 IF( KP.LT.N ) 527 $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), 528 $ 1 ) 529 KX = KNC + KP - KK 530 DO 80 J = KK + 1, KP - 1 531 KX = KX + N - J + 1 532 T = DCONJG( AP( KNC+J-KK ) ) 533 AP( KNC+J-KK ) = DCONJG( AP( KX ) ) 534 AP( KX ) = T 535 80 CONTINUE 536 AP( KNC+KP-KK ) = DCONJG( AP( KNC+KP-KK ) ) 537 R1 = DBLE( AP( KNC ) ) 538 AP( KNC ) = DBLE( AP( KPC ) ) 539 AP( KPC ) = R1 540 IF( KSTEP.EQ.2 ) THEN 541 AP( KC ) = DBLE( AP( KC ) ) 542 T = AP( KC+1 ) 543 AP( KC+1 ) = AP( KC+KP-K ) 544 AP( KC+KP-K ) = T 545 END IF 546 ELSE 547 AP( KC ) = DBLE( AP( KC ) ) 548 IF( KSTEP.EQ.2 ) 549 $ AP( KNC ) = DBLE( AP( KNC ) ) 550 END IF 551* 552* Update the trailing submatrix 553* 554 IF( KSTEP.EQ.1 ) THEN 555* 556* 1-by-1 pivot block D(k): column k now holds 557* 558* W(k) = L(k)*D(k) 559* 560* where L(k) is the k-th column of L 561* 562 IF( K.LT.N ) THEN 563* 564* Perform a rank-1 update of A(k+1:n,k+1:n) as 565* 566* A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H 567* 568 R1 = ONE / DBLE( AP( KC ) ) 569 CALL ZHPR( UPLO, N-K, -R1, AP( KC+1 ), 1, 570 $ AP( KC+N-K+1 ) ) 571* 572* Store L(k) in column K 573* 574 CALL ZDSCAL( N-K, R1, AP( KC+1 ), 1 ) 575 END IF 576 ELSE 577* 578* 2-by-2 pivot block D(k): columns K and K+1 now hold 579* 580* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 581* 582* where L(k) and L(k+1) are the k-th and (k+1)-th columns 583* of L 584* 585 IF( K.LT.N-1 ) THEN 586* 587* Perform a rank-2 update of A(k+2:n,k+2:n) as 588* 589* A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H 590* = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H 591* 592* where L(k) and L(k+1) are the k-th and (k+1)-th 593* columns of L 594* 595 D = DLAPY2( DBLE( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ), 596 $ DIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) ) 597 D11 = DBLE( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D 598 D22 = DBLE( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D 599 TT = ONE / ( D11*D22-ONE ) 600 D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D 601 D = TT / D 602* 603 DO 100 J = K + 2, N 604 WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21* 605 $ AP( J+K*( 2*N-K-1 ) / 2 ) ) 606 WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- 607 $ DCONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 608 $ 2 ) ) 609 DO 90 I = J, N 610 AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* 611 $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / 612 $ 2 )*DCONJG( WK ) - AP( I+K*( 2*N-K-1 ) / 2 )* 613 $ DCONJG( WKP1 ) 614 90 CONTINUE 615 AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK 616 AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1 617 AP( J+( J-1 )*( 2*N-J ) / 2 ) 618 $ = DCMPLX( DBLE( AP( J+( J-1 )*( 2*N-J ) / 2 ) ), 619 $ 0.0D+0 ) 620 100 CONTINUE 621 END IF 622 END IF 623 END IF 624* 625* Store details of the interchanges in IPIV 626* 627 IF( KSTEP.EQ.1 ) THEN 628 IPIV( K ) = KP 629 ELSE 630 IPIV( K ) = -KP 631 IPIV( K+1 ) = -KP 632 END IF 633* 634* Increase K and return to the start of the main loop 635* 636 K = K + KSTEP 637 KC = KNC + N - K + 2 638 GO TO 60 639* 640 END IF 641* 642 110 CONTINUE 643 RETURN 644* 645* End of ZHPTRF 646* 647 END 648