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