1 SUBROUTINE DSYTRI( UPLO, N, A, LDA, IPIV, WORK, 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, 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**T.
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**T.
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
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, 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**T.
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**T.
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