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