001:       SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INFO, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       COMPLEX            AP( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  CHPTRF computes the factorization of a complex Hermitian packed
020: *  matrix A using the Bunch-Kaufman diagonal pivoting method:
021: *
022: *     A = U*D*U**H  or  A = L*D*L**H
023: *
024: *  where U (or L) is a product of permutation and unit upper (lower)
025: *  triangular matrices, and D is Hermitian and block diagonal with
026: *  1-by-1 and 2-by-2 diagonal blocks.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  UPLO    (input) CHARACTER*1
032: *          = 'U':  Upper triangle of A is stored;
033: *          = 'L':  Lower triangle of A is stored.
034: *
035: *  N       (input) INTEGER
036: *          The order of the matrix A.  N >= 0.
037: *
038: *  AP      (input/output) COMPLEX array, dimension (N*(N+1)/2)
039: *          On entry, the upper or lower triangle of the Hermitian matrix
040: *          A, packed columnwise in a linear array.  The j-th column of A
041: *          is stored in the array AP as follows:
042: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
043: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
044: *
045: *          On exit, the block diagonal matrix D and the multipliers used
046: *          to obtain the factor U or L, stored as a packed triangular
047: *          matrix overwriting A (see below for further details).
048: *
049: *  IPIV    (output) INTEGER array, dimension (N)
050: *          Details of the interchanges and the block structure of D.
051: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
052: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
053: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
054: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
055: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
056: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
057: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
058: *
059: *  INFO    (output) INTEGER
060: *          = 0: successful exit
061: *          < 0: if INFO = -i, the i-th argument had an illegal value
062: *          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
063: *               has been completed, but the block diagonal matrix D is
064: *               exactly singular, and division by zero will occur if it
065: *               is used to solve a system of equations.
066: *
067: *  Further Details
068: *  ===============
069: *
070: *  5-96 - Based on modifications by J. Lewis, Boeing Computer Services
071: *         Company
072: *
073: *  If UPLO = 'U', then A = U*D*U', where
074: *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
075: *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
076: *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
077: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
078: *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
079: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
080: *
081: *             (   I    v    0   )   k-s
082: *     U(k) =  (   0    I    0   )   s
083: *             (   0    0    I   )   n-k
084: *                k-s   s   n-k
085: *
086: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
087: *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
088: *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
089: *
090: *  If UPLO = 'L', then A = L*D*L', where
091: *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
092: *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
093: *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
094: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
095: *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
096: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
097: *
098: *             (   I    0     0   )  k-1
099: *     L(k) =  (   0    I     0   )  s
100: *             (   0    v     I   )  n-k-s+1
101: *                k-1   s  n-k-s+1
102: *
103: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
104: *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
105: *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
106: *
107: *  =====================================================================
108: *
109: *     .. Parameters ..
110:       REAL               ZERO, ONE
111:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
112:       REAL               EIGHT, SEVTEN
113:       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
114: *     ..
115: *     .. Local Scalars ..
116:       LOGICAL            UPPER
117:       INTEGER            I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
118:      $                   KSTEP, KX, NPP
119:       REAL               ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
120:      $                   TT
121:       COMPLEX            D12, D21, T, WK, WKM1, WKP1, ZDUM
122: *     ..
123: *     .. External Functions ..
124:       LOGICAL            LSAME
125:       INTEGER            ICAMAX
126:       REAL               SLAPY2
127:       EXTERNAL           LSAME, ICAMAX, SLAPY2
128: *     ..
129: *     .. External Subroutines ..
130:       EXTERNAL           CHPR, CSSCAL, CSWAP, XERBLA
131: *     ..
132: *     .. Intrinsic Functions ..
133:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
134: *     ..
135: *     .. Statement Functions ..
136:       REAL               CABS1
137: *     ..
138: *     .. Statement Function definitions ..
139:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
140: *     ..
141: *     .. Executable Statements ..
142: *
143: *     Test the input parameters.
144: *
145:       INFO = 0
146:       UPPER = LSAME( UPLO, 'U' )
147:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
148:          INFO = -1
149:       ELSE IF( N.LT.0 ) THEN
150:          INFO = -2
151:       END IF
152:       IF( INFO.NE.0 ) THEN
153:          CALL XERBLA( 'CHPTRF', -INFO )
154:          RETURN
155:       END IF
156: *
157: *     Initialize ALPHA for use in choosing pivot block size.
158: *
159:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
160: *
161:       IF( UPPER ) THEN
162: *
163: *        Factorize A as U*D*U' using the upper triangle of A
164: *
165: *        K is the main loop index, decreasing from N to 1 in steps of
166: *        1 or 2
167: *
168:          K = N
169:          KC = ( N-1 )*N / 2 + 1
170:    10    CONTINUE
171:          KNC = KC
172: *
173: *        If K < 1, exit from loop
174: *
175:          IF( K.LT.1 )
176:      $      GO TO 110
177:          KSTEP = 1
178: *
179: *        Determine rows and columns to be interchanged and whether
180: *        a 1-by-1 or 2-by-2 pivot block will be used
181: *
182:          ABSAKK = ABS( REAL( AP( KC+K-1 ) ) )
183: *
184: *        IMAX is the row-index of the largest off-diagonal element in
185: *        column K, and COLMAX is its absolute value
186: *
187:          IF( K.GT.1 ) THEN
188:             IMAX = ICAMAX( K-1, AP( KC ), 1 )
189:             COLMAX = CABS1( AP( KC+IMAX-1 ) )
190:          ELSE
191:             COLMAX = ZERO
192:          END IF
193: *
194:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
195: *
196: *           Column K is zero: set INFO and continue
197: *
198:             IF( INFO.EQ.0 )
199:      $         INFO = K
200:             KP = K
201:             AP( KC+K-1 ) = REAL( AP( KC+K-1 ) )
202:          ELSE
203:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
204: *
205: *              no interchange, use 1-by-1 pivot block
206: *
207:                KP = K
208:             ELSE
209: *
210: *              JMAX is the column-index of the largest off-diagonal
211: *              element in row IMAX, and ROWMAX is its absolute value
212: *
213:                ROWMAX = ZERO
214:                JMAX = IMAX
215:                KX = IMAX*( IMAX+1 ) / 2 + IMAX
216:                DO 20 J = IMAX + 1, K
217:                   IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
218:                      ROWMAX = CABS1( AP( KX ) )
219:                      JMAX = J
220:                   END IF
221:                   KX = KX + J
222:    20          CONTINUE
223:                KPC = ( IMAX-1 )*IMAX / 2 + 1
224:                IF( IMAX.GT.1 ) THEN
225:                   JMAX = ICAMAX( IMAX-1, AP( KPC ), 1 )
226:                   ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-1 ) ) )
227:                END IF
228: *
229:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
230: *
231: *                 no interchange, use 1-by-1 pivot block
232: *
233:                   KP = K
234:                ELSE IF( ABS( REAL( AP( KPC+IMAX-1 ) ) ).GE.ALPHA*
235:      $                  ROWMAX ) THEN
236: *
237: *                 interchange rows and columns K and IMAX, use 1-by-1
238: *                 pivot block
239: *
240:                   KP = IMAX
241:                ELSE
242: *
243: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
244: *                 pivot block
245: *
246:                   KP = IMAX
247:                   KSTEP = 2
248:                END IF
249:             END IF
250: *
251:             KK = K - KSTEP + 1
252:             IF( KSTEP.EQ.2 )
253:      $         KNC = KNC - K + 1
254:             IF( KP.NE.KK ) THEN
255: *
256: *              Interchange rows and columns KK and KP in the leading
257: *              submatrix A(1:k,1:k)
258: *
259:                CALL CSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
260:                KX = KPC + KP - 1
261:                DO 30 J = KP + 1, KK - 1
262:                   KX = KX + J - 1
263:                   T = CONJG( AP( KNC+J-1 ) )
264:                   AP( KNC+J-1 ) = CONJG( AP( KX ) )
265:                   AP( KX ) = T
266:    30          CONTINUE
267:                AP( KX+KK-1 ) = CONJG( AP( KX+KK-1 ) )
268:                R1 = REAL( AP( KNC+KK-1 ) )
269:                AP( KNC+KK-1 ) = REAL( AP( KPC+KP-1 ) )
270:                AP( KPC+KP-1 ) = R1
271:                IF( KSTEP.EQ.2 ) THEN
272:                   AP( KC+K-1 ) = REAL( AP( KC+K-1 ) )
273:                   T = AP( KC+K-2 )
274:                   AP( KC+K-2 ) = AP( KC+KP-1 )
275:                   AP( KC+KP-1 ) = T
276:                END IF
277:             ELSE
278:                AP( KC+K-1 ) = REAL( AP( KC+K-1 ) )
279:                IF( KSTEP.EQ.2 )
280:      $            AP( KC-1 ) = REAL( AP( KC-1 ) )
281:             END IF
282: *
283: *           Update the leading submatrix
284: *
285:             IF( KSTEP.EQ.1 ) THEN
286: *
287: *              1-by-1 pivot block D(k): column k now holds
288: *
289: *              W(k) = U(k)*D(k)
290: *
291: *              where U(k) is the k-th column of U
292: *
293: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
294: *
295: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
296: *
297:                R1 = ONE / REAL( AP( KC+K-1 ) )
298:                CALL CHPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
299: *
300: *              Store U(k) in column k
301: *
302:                CALL CSSCAL( K-1, R1, AP( KC ), 1 )
303:             ELSE
304: *
305: *              2-by-2 pivot block D(k): columns k and k-1 now hold
306: *
307: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
308: *
309: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
310: *              of U
311: *
312: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
313: *
314: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
315: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
316: *
317:                IF( K.GT.2 ) THEN
318: *
319:                   D = SLAPY2( REAL( AP( K-1+( K-1 )*K / 2 ) ),
320:      $                AIMAG( AP( K-1+( K-1 )*K / 2 ) ) )
321:                   D22 = REAL( AP( K-1+( K-2 )*( K-1 ) / 2 ) ) / D
322:                   D11 = REAL( AP( K+( K-1 )*K / 2 ) ) / D
323:                   TT = ONE / ( D11*D22-ONE )
324:                   D12 = AP( K-1+( K-1 )*K / 2 ) / D
325:                   D = TT / D
326: *
327:                   DO 50 J = K - 2, 1, -1
328:                      WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
329:      $                      CONJG( D12 )*AP( J+( K-1 )*K / 2 ) )
330:                      WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12*
331:      $                    AP( J+( K-2 )*( K-1 ) / 2 ) )
332:                      DO 40 I = J, 1, -1
333:                         AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
334:      $                     AP( I+( K-1 )*K / 2 )*CONJG( WK ) -
335:      $                     AP( I+( K-2 )*( K-1 ) / 2 )*CONJG( WKM1 )
336:    40                CONTINUE
337:                      AP( J+( K-1 )*K / 2 ) = WK
338:                      AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
339:                      AP( J+( J-1 )*J / 2 ) = CMPLX( REAL( AP( J+( J-1 )*
340:      $                                       J / 2 ) ), 0.0E+0 )
341:    50             CONTINUE
342: *
343:                END IF
344: *
345:             END IF
346:          END IF
347: *
348: *        Store details of the interchanges in IPIV
349: *
350:          IF( KSTEP.EQ.1 ) THEN
351:             IPIV( K ) = KP
352:          ELSE
353:             IPIV( K ) = -KP
354:             IPIV( K-1 ) = -KP
355:          END IF
356: *
357: *        Decrease K and return to the start of the main loop
358: *
359:          K = K - KSTEP
360:          KC = KNC - K
361:          GO TO 10
362: *
363:       ELSE
364: *
365: *        Factorize A as L*D*L' using the lower triangle of A
366: *
367: *        K is the main loop index, increasing from 1 to N in steps of
368: *        1 or 2
369: *
370:          K = 1
371:          KC = 1
372:          NPP = N*( N+1 ) / 2
373:    60    CONTINUE
374:          KNC = KC
375: *
376: *        If K > N, exit from loop
377: *
378:          IF( K.GT.N )
379:      $      GO TO 110
380:          KSTEP = 1
381: *
382: *        Determine rows and columns to be interchanged and whether
383: *        a 1-by-1 or 2-by-2 pivot block will be used
384: *
385:          ABSAKK = ABS( REAL( AP( KC ) ) )
386: *
387: *        IMAX is the row-index of the largest off-diagonal element in
388: *        column K, and COLMAX is its absolute value
389: *
390:          IF( K.LT.N ) THEN
391:             IMAX = K + ICAMAX( N-K, AP( KC+1 ), 1 )
392:             COLMAX = CABS1( AP( KC+IMAX-K ) )
393:          ELSE
394:             COLMAX = ZERO
395:          END IF
396: *
397:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
398: *
399: *           Column K is zero: set INFO and continue
400: *
401:             IF( INFO.EQ.0 )
402:      $         INFO = K
403:             KP = K
404:             AP( KC ) = REAL( AP( KC ) )
405:          ELSE
406:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
407: *
408: *              no interchange, use 1-by-1 pivot block
409: *
410:                KP = K
411:             ELSE
412: *
413: *              JMAX is the column-index of the largest off-diagonal
414: *              element in row IMAX, and ROWMAX is its absolute value
415: *
416:                ROWMAX = ZERO
417:                KX = KC + IMAX - K
418:                DO 70 J = K, IMAX - 1
419:                   IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
420:                      ROWMAX = CABS1( AP( KX ) )
421:                      JMAX = J
422:                   END IF
423:                   KX = KX + N - J
424:    70          CONTINUE
425:                KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
426:                IF( IMAX.LT.N ) THEN
427:                   JMAX = IMAX + ICAMAX( N-IMAX, AP( KPC+1 ), 1 )
428:                   ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-IMAX ) ) )
429:                END IF
430: *
431:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
432: *
433: *                 no interchange, use 1-by-1 pivot block
434: *
435:                   KP = K
436:                ELSE IF( ABS( REAL( AP( KPC ) ) ).GE.ALPHA*ROWMAX ) THEN
437: *
438: *                 interchange rows and columns K and IMAX, use 1-by-1
439: *                 pivot block
440: *
441:                   KP = IMAX
442:                ELSE
443: *
444: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
445: *                 pivot block
446: *
447:                   KP = IMAX
448:                   KSTEP = 2
449:                END IF
450:             END IF
451: *
452:             KK = K + KSTEP - 1
453:             IF( KSTEP.EQ.2 )
454:      $         KNC = KNC + N - K + 1
455:             IF( KP.NE.KK ) THEN
456: *
457: *              Interchange rows and columns KK and KP in the trailing
458: *              submatrix A(k:n,k:n)
459: *
460:                IF( KP.LT.N )
461:      $            CALL CSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
462:      $                        1 )
463:                KX = KNC + KP - KK
464:                DO 80 J = KK + 1, KP - 1
465:                   KX = KX + N - J + 1
466:                   T = CONJG( AP( KNC+J-KK ) )
467:                   AP( KNC+J-KK ) = CONJG( AP( KX ) )
468:                   AP( KX ) = T
469:    80          CONTINUE
470:                AP( KNC+KP-KK ) = CONJG( AP( KNC+KP-KK ) )
471:                R1 = REAL( AP( KNC ) )
472:                AP( KNC ) = REAL( AP( KPC ) )
473:                AP( KPC ) = R1
474:                IF( KSTEP.EQ.2 ) THEN
475:                   AP( KC ) = REAL( AP( KC ) )
476:                   T = AP( KC+1 )
477:                   AP( KC+1 ) = AP( KC+KP-K )
478:                   AP( KC+KP-K ) = T
479:                END IF
480:             ELSE
481:                AP( KC ) = REAL( AP( KC ) )
482:                IF( KSTEP.EQ.2 )
483:      $            AP( KNC ) = REAL( AP( KNC ) )
484:             END IF
485: *
486: *           Update the trailing submatrix
487: *
488:             IF( KSTEP.EQ.1 ) THEN
489: *
490: *              1-by-1 pivot block D(k): column k now holds
491: *
492: *              W(k) = L(k)*D(k)
493: *
494: *              where L(k) is the k-th column of L
495: *
496:                IF( K.LT.N ) THEN
497: *
498: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
499: *
500: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
501: *
502:                   R1 = ONE / REAL( AP( KC ) )
503:                   CALL CHPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
504:      $                       AP( KC+N-K+1 ) )
505: *
506: *                 Store L(k) in column K
507: *
508:                   CALL CSSCAL( N-K, R1, AP( KC+1 ), 1 )
509:                END IF
510:             ELSE
511: *
512: *              2-by-2 pivot block D(k): columns K and K+1 now hold
513: *
514: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
515: *
516: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
517: *              of L
518: *
519:                IF( K.LT.N-1 ) THEN
520: *
521: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
522: *
523: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
524: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
525: *
526: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
527: *                 columns of L
528: *
529:                   D = SLAPY2( REAL( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ),
530:      $                AIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) )
531:                   D11 = REAL( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D
532:                   D22 = REAL( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D
533:                   TT = ONE / ( D11*D22-ONE )
534:                   D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D
535:                   D = TT / D
536: *
537:                   DO 100 J = K + 2, N
538:                      WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21*
539:      $                    AP( J+K*( 2*N-K-1 ) / 2 ) )
540:                      WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
541:      $                      CONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 2 ) )
542:                      DO 90 I = J, N
543:                         AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
544:      $                     ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
545:      $                     2 )*CONJG( WK ) - AP( I+K*( 2*N-K-1 ) / 2 )*
546:      $                     CONJG( WKP1 )
547:    90                CONTINUE
548:                      AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
549:                      AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
550:                      AP( J+( J-1 )*( 2*N-J ) / 2 )
551:      $                  = CMPLX( REAL( AP( J+( J-1 )*( 2*N-J ) / 2 ) ),
552:      $                  0.0E+0 )
553:   100             CONTINUE
554:                END IF
555:             END IF
556:          END IF
557: *
558: *        Store details of the interchanges in IPIV
559: *
560:          IF( KSTEP.EQ.1 ) THEN
561:             IPIV( K ) = KP
562:          ELSE
563:             IPIV( K ) = -KP
564:             IPIV( K+1 ) = -KP
565:          END IF
566: *
567: *        Increase K and return to the start of the main loop
568: *
569:          K = K + KSTEP
570:          KC = KNC + N - K + 2
571:          GO TO 60
572: *
573:       END IF
574: *
575:   110 CONTINUE
576:       RETURN
577: *
578: *     End of CHPTRF
579: *
580:       END
581: