![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsptri.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 dsptri.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, N
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: DOUBLE PRECISION AP( * ), WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DSPTRI computes the inverse of a real symmetric indefinite matrix
21: * A in packed storage using the factorization A = U*D*U**T or
22: * A = L*D*L**T computed by DSPTRF.
23: *
24: * Arguments
25: * =========
26: *
27: * UPLO (input) CHARACTER*1
28: * Specifies whether the details of the factorization are stored
29: * as an upper or lower triangular matrix.
30: * = 'U': Upper triangular, form is A = U*D*U**T;
31: * = 'L': Lower triangular, form is A = L*D*L**T.
32: *
33: * N (input) INTEGER
34: * The order of the matrix A. N >= 0.
35: *
36: * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
37: * On entry, the block diagonal matrix D and the multipliers
38: * used to obtain the factor U or L as computed by DSPTRF,
39: * stored as a packed triangular matrix.
40: *
41: * On exit, if INFO = 0, the (symmetric) inverse of the original
42: * matrix, stored as a packed triangular matrix. The j-th column
43: * of inv(A) is stored in the array AP as follows:
44: * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45: * if UPLO = 'L',
46: * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47: *
48: * IPIV (input) INTEGER array, dimension (N)
49: * Details of the interchanges and the block structure of D
50: * as determined by DSPTRF.
51: *
52: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
53: *
54: * INFO (output) INTEGER
55: * = 0: successful exit
56: * < 0: if INFO = -i, the i-th argument had an illegal value
57: * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58: * inverse could not be computed.
59: *
60: * =====================================================================
61: *
62: * .. Parameters ..
63: DOUBLE PRECISION ONE, ZERO
64: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
65: * ..
66: * .. Local Scalars ..
67: LOGICAL UPPER
68: INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
69: DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
70: * ..
71: * .. External Functions ..
72: LOGICAL LSAME
73: DOUBLE PRECISION DDOT
74: EXTERNAL LSAME, DDOT
75: * ..
76: * .. External Subroutines ..
77: EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA
78: * ..
79: * .. Intrinsic Functions ..
80: INTRINSIC ABS
81: * ..
82: * .. Executable Statements ..
83: *
84: * Test the input parameters.
85: *
86: INFO = 0
87: UPPER = LSAME( UPLO, 'U' )
88: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89: INFO = -1
90: ELSE IF( N.LT.0 ) THEN
91: INFO = -2
92: END IF
93: IF( INFO.NE.0 ) THEN
94: CALL XERBLA( 'DSPTRI', -INFO )
95: RETURN
96: END IF
97: *
98: * Quick return if possible
99: *
100: IF( N.EQ.0 )
101: $ RETURN
102: *
103: * Check that the diagonal matrix D is nonsingular.
104: *
105: IF( UPPER ) THEN
106: *
107: * Upper triangular storage: examine D from bottom to top
108: *
109: KP = N*( N+1 ) / 2
110: DO 10 INFO = N, 1, -1
111: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
112: $ RETURN
113: KP = KP - INFO
114: 10 CONTINUE
115: ELSE
116: *
117: * Lower triangular storage: examine D from top to bottom.
118: *
119: KP = 1
120: DO 20 INFO = 1, N
121: IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
122: $ RETURN
123: KP = KP + N - INFO + 1
124: 20 CONTINUE
125: END IF
126: INFO = 0
127: *
128: IF( UPPER ) THEN
129: *
130: * Compute inv(A) from the factorization A = U*D*U'.
131: *
132: * K is the main loop index, increasing from 1 to N in steps of
133: * 1 or 2, depending on the size of the diagonal blocks.
134: *
135: K = 1
136: KC = 1
137: 30 CONTINUE
138: *
139: * If K > N, exit from loop.
140: *
141: IF( K.GT.N )
142: $ GO TO 50
143: *
144: KCNEXT = KC + K
145: IF( IPIV( K ).GT.0 ) THEN
146: *
147: * 1 x 1 diagonal block
148: *
149: * Invert the diagonal block.
150: *
151: AP( KC+K-1 ) = ONE / AP( KC+K-1 )
152: *
153: * Compute column K of the inverse.
154: *
155: IF( K.GT.1 ) THEN
156: CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
157: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
158: $ 1 )
159: AP( KC+K-1 ) = AP( KC+K-1 ) -
160: $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
161: END IF
162: KSTEP = 1
163: ELSE
164: *
165: * 2 x 2 diagonal block
166: *
167: * Invert the diagonal block.
168: *
169: T = ABS( AP( KCNEXT+K-1 ) )
170: AK = AP( KC+K-1 ) / T
171: AKP1 = AP( KCNEXT+K ) / T
172: AKKP1 = AP( KCNEXT+K-1 ) / T
173: D = T*( AK*AKP1-ONE )
174: AP( KC+K-1 ) = AKP1 / D
175: AP( KCNEXT+K ) = AK / D
176: AP( KCNEXT+K-1 ) = -AKKP1 / D
177: *
178: * Compute columns K and K+1 of the inverse.
179: *
180: IF( K.GT.1 ) THEN
181: CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
182: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
183: $ 1 )
184: AP( KC+K-1 ) = AP( KC+K-1 ) -
185: $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
186: AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
187: $ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
188: $ 1 )
189: CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
190: CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
191: $ AP( KCNEXT ), 1 )
192: AP( KCNEXT+K ) = AP( KCNEXT+K ) -
193: $ DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 )
194: END IF
195: KSTEP = 2
196: KCNEXT = KCNEXT + K + 1
197: END IF
198: *
199: KP = ABS( IPIV( K ) )
200: IF( KP.NE.K ) THEN
201: *
202: * Interchange rows and columns K and KP in the leading
203: * submatrix A(1:k+1,1:k+1)
204: *
205: KPC = ( KP-1 )*KP / 2 + 1
206: CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
207: KX = KPC + KP - 1
208: DO 40 J = KP + 1, K - 1
209: KX = KX + J - 1
210: TEMP = AP( KC+J-1 )
211: AP( KC+J-1 ) = AP( KX )
212: AP( KX ) = TEMP
213: 40 CONTINUE
214: TEMP = AP( KC+K-1 )
215: AP( KC+K-1 ) = AP( KPC+KP-1 )
216: AP( KPC+KP-1 ) = TEMP
217: IF( KSTEP.EQ.2 ) THEN
218: TEMP = AP( KC+K+K-1 )
219: AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
220: AP( KC+K+KP-1 ) = TEMP
221: END IF
222: END IF
223: *
224: K = K + KSTEP
225: KC = KCNEXT
226: GO TO 30
227: 50 CONTINUE
228: *
229: ELSE
230: *
231: * Compute inv(A) from the factorization A = L*D*L'.
232: *
233: * K is the main loop index, increasing from 1 to N in steps of
234: * 1 or 2, depending on the size of the diagonal blocks.
235: *
236: NPP = N*( N+1 ) / 2
237: K = N
238: KC = NPP
239: 60 CONTINUE
240: *
241: * If K < 1, exit from loop.
242: *
243: IF( K.LT.1 )
244: $ GO TO 80
245: *
246: KCNEXT = KC - ( N-K+2 )
247: IF( IPIV( K ).GT.0 ) THEN
248: *
249: * 1 x 1 diagonal block
250: *
251: * Invert the diagonal block.
252: *
253: AP( KC ) = ONE / AP( KC )
254: *
255: * Compute column K of the inverse.
256: *
257: IF( K.LT.N ) THEN
258: CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
259: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
260: $ ZERO, AP( KC+1 ), 1 )
261: AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
262: END IF
263: KSTEP = 1
264: ELSE
265: *
266: * 2 x 2 diagonal block
267: *
268: * Invert the diagonal block.
269: *
270: T = ABS( AP( KCNEXT+1 ) )
271: AK = AP( KCNEXT ) / T
272: AKP1 = AP( KC ) / T
273: AKKP1 = AP( KCNEXT+1 ) / T
274: D = T*( AK*AKP1-ONE )
275: AP( KCNEXT ) = AKP1 / D
276: AP( KC ) = AK / D
277: AP( KCNEXT+1 ) = -AKKP1 / D
278: *
279: * Compute columns K-1 and K of the inverse.
280: *
281: IF( K.LT.N ) THEN
282: CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
283: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
284: $ ZERO, AP( KC+1 ), 1 )
285: AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
286: AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
287: $ DDOT( N-K, AP( KC+1 ), 1,
288: $ AP( KCNEXT+2 ), 1 )
289: CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
290: CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
291: $ ZERO, AP( KCNEXT+2 ), 1 )
292: AP( KCNEXT ) = AP( KCNEXT ) -
293: $ DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
294: END IF
295: KSTEP = 2
296: KCNEXT = KCNEXT - ( N-K+3 )
297: END IF
298: *
299: KP = ABS( IPIV( K ) )
300: IF( KP.NE.K ) THEN
301: *
302: * Interchange rows and columns K and KP in the trailing
303: * submatrix A(k-1:n,k-1:n)
304: *
305: KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
306: IF( KP.LT.N )
307: $ CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
308: KX = KC + KP - K
309: DO 70 J = K + 1, KP - 1
310: KX = KX + N - J + 1
311: TEMP = AP( KC+J-K )
312: AP( KC+J-K ) = AP( KX )
313: AP( KX ) = TEMP
314: 70 CONTINUE
315: TEMP = AP( KC )
316: AP( KC ) = AP( KPC )
317: AP( KPC ) = TEMP
318: IF( KSTEP.EQ.2 ) THEN
319: TEMP = AP( KC-N+K-1 )
320: AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
321: AP( KC-N+KP-1 ) = TEMP
322: END IF
323: END IF
324: *
325: K = K - KSTEP
326: KC = KCNEXT
327: GO TO 60
328: 80 CONTINUE
329: END IF
330: *
331: RETURN
332: *
333: * End of DSPTRI
334: *
335: END