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