Annotation of rpl/lapack/lapack/zsptrf.f, revision 1.16
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: *
1.15 bertrand 111: *> \date December 2016
1.9 bertrand 112: *
113: *> \ingroup complex16OTHERcomputational
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: * =====================================================================
1.1 bertrand 159: SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO )
160: *
1.15 bertrand 161: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 162: * -- LAPACK is a software package provided by Univ. of Tennessee, --
163: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 bertrand 164: * December 2016
1.1 bertrand 165: *
166: * .. Scalar Arguments ..
167: CHARACTER UPLO
168: INTEGER INFO, N
169: * ..
170: * .. Array Arguments ..
171: INTEGER IPIV( * )
172: COMPLEX*16 AP( * )
173: * ..
174: *
175: * =====================================================================
176: *
177: * .. Parameters ..
178: DOUBLE PRECISION ZERO, ONE
179: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
180: DOUBLE PRECISION EIGHT, SEVTEN
181: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
182: COMPLEX*16 CONE
183: PARAMETER ( CONE = ( 1.0D+0, 0.0D+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: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
190: COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM
191: * ..
192: * .. External Functions ..
193: LOGICAL LSAME
194: INTEGER IZAMAX
195: EXTERNAL LSAME, IZAMAX
196: * ..
197: * .. External Subroutines ..
198: EXTERNAL XERBLA, ZSCAL, ZSPR, ZSWAP
199: * ..
200: * .. Intrinsic Functions ..
201: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
202: * ..
203: * .. Statement Functions ..
204: DOUBLE PRECISION CABS1
205: * ..
206: * .. Statement Function definitions ..
207: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( 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( 'ZSPTRF', -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: *
1.8 bertrand 231: * Factorize A as U*D*U**T using the upper triangle of A
1.1 bertrand 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 = IZAMAX( 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 = IZAMAX( 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 ZSWAP( 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: *
1.8 bertrand 352: * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
1.1 bertrand 353: *
354: R1 = CONE / AP( KC+K-1 )
355: CALL ZSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
356: *
357: * Store U(k) in column k
358: *
359: CALL ZSCAL( 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: *
1.8 bertrand 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
1.1 bertrand 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: *
1.8 bertrand 417: * Factorize A as L*D*L**T using the lower triangle of A
1.1 bertrand 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 + IZAMAX( 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 + IZAMAX( 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 ZSWAP( 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: *
1.8 bertrand 545: * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
1.1 bertrand 546: *
547: R1 = CONE / AP( KC )
548: CALL ZSPR( 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 ZSCAL( 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: *
1.8 bertrand 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
1.1 bertrand 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 ZSPTRF
618: *
619: 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