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