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