![[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|
Return to zhptrf.f CVS log | Up to [local] / rpl / lapack / lapack |
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