1 SUBROUTINE ZSYTRI( 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 COMPLEX*16 A( LDA, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZSYTRI computes the inverse of a complex symmetric indefinite matrix
21 * A using the factorization A = U*D*U**T or A = L*D*L**T computed by
22 * ZSYTRF.
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) COMPLEX*16 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 ZSYTRF.
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 ZSYTRF.
53 *
54 * WORK (workspace) COMPLEX*16 array, dimension (2*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 COMPLEX*16 ONE, ZERO
66 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
67 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 LOGICAL UPPER
71 INTEGER K, KP, KSTEP
72 COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 COMPLEX*16 ZDOTU
77 EXTERNAL LSAME, ZDOTU
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA, ZCOPY, ZSWAP, ZSYMV
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, MAX
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
96 INFO = -4
97 END IF
98 IF( INFO.NE.0 ) THEN
99 CALL XERBLA( 'ZSYTRI', -INFO )
100 RETURN
101 END IF
102 *
103 * Quick return if possible
104 *
105 IF( N.EQ.0 )
106 $ RETURN
107 *
108 * Check that the diagonal matrix D is nonsingular.
109 *
110 IF( UPPER ) THEN
111 *
112 * Upper triangular storage: examine D from bottom to top
113 *
114 DO 10 INFO = N, 1, -1
115 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
116 $ RETURN
117 10 CONTINUE
118 ELSE
119 *
120 * Lower triangular storage: examine D from top to bottom.
121 *
122 DO 20 INFO = 1, N
123 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
124 $ RETURN
125 20 CONTINUE
126 END IF
127 INFO = 0
128 *
129 IF( UPPER ) THEN
130 *
131 * Compute inv(A) from the factorization A = U*D*U**T.
132 *
133 * K is the main loop index, increasing from 1 to N in steps of
134 * 1 or 2, depending on the size of the diagonal blocks.
135 *
136 K = 1
137 30 CONTINUE
138 *
139 * If K > N, exit from loop.
140 *
141 IF( K.GT.N )
142 $ GO TO 40
143 *
144 IF( IPIV( K ).GT.0 ) THEN
145 *
146 * 1 x 1 diagonal block
147 *
148 * Invert the diagonal block.
149 *
150 A( K, K ) = ONE / A( K, K )
151 *
152 * Compute column K of the inverse.
153 *
154 IF( K.GT.1 ) THEN
155 CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
156 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
157 $ A( 1, K ), 1 )
158 A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
159 $ 1 )
160 END IF
161 KSTEP = 1
162 ELSE
163 *
164 * 2 x 2 diagonal block
165 *
166 * Invert the diagonal block.
167 *
168 T = A( K, K+1 )
169 AK = A( K, K ) / T
170 AKP1 = A( K+1, K+1 ) / T
171 AKKP1 = A( K, K+1 ) / T
172 D = T*( AK*AKP1-ONE )
173 A( K, K ) = AKP1 / D
174 A( K+1, K+1 ) = AK / D
175 A( K, K+1 ) = -AKKP1 / D
176 *
177 * Compute columns K and K+1 of the inverse.
178 *
179 IF( K.GT.1 ) THEN
180 CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
181 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
182 $ A( 1, K ), 1 )
183 A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
184 $ 1 )
185 A( K, K+1 ) = A( K, K+1 ) -
186 $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
187 CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
188 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
189 $ A( 1, K+1 ), 1 )
190 A( K+1, K+1 ) = A( K+1, K+1 ) -
191 $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
192 END IF
193 KSTEP = 2
194 END IF
195 *
196 KP = ABS( IPIV( K ) )
197 IF( KP.NE.K ) THEN
198 *
199 * Interchange rows and columns K and KP in the leading
200 * submatrix A(1:k+1,1:k+1)
201 *
202 CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
203 CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
204 TEMP = A( K, K )
205 A( K, K ) = A( KP, KP )
206 A( KP, KP ) = TEMP
207 IF( KSTEP.EQ.2 ) THEN
208 TEMP = A( K, K+1 )
209 A( K, K+1 ) = A( KP, K+1 )
210 A( KP, K+1 ) = TEMP
211 END IF
212 END IF
213 *
214 K = K + KSTEP
215 GO TO 30
216 40 CONTINUE
217 *
218 ELSE
219 *
220 * Compute inv(A) from the factorization A = L*D*L**T.
221 *
222 * K is the main loop index, increasing from 1 to N in steps of
223 * 1 or 2, depending on the size of the diagonal blocks.
224 *
225 K = N
226 50 CONTINUE
227 *
228 * If K < 1, exit from loop.
229 *
230 IF( K.LT.1 )
231 $ GO TO 60
232 *
233 IF( IPIV( K ).GT.0 ) THEN
234 *
235 * 1 x 1 diagonal block
236 *
237 * Invert the diagonal block.
238 *
239 A( K, K ) = ONE / A( K, K )
240 *
241 * Compute column K of the inverse.
242 *
243 IF( K.LT.N ) THEN
244 CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
245 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
246 $ ZERO, A( K+1, K ), 1 )
247 A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
248 $ 1 )
249 END IF
250 KSTEP = 1
251 ELSE
252 *
253 * 2 x 2 diagonal block
254 *
255 * Invert the diagonal block.
256 *
257 T = A( K, K-1 )
258 AK = A( K-1, K-1 ) / T
259 AKP1 = A( K, K ) / T
260 AKKP1 = A( K, K-1 ) / T
261 D = T*( AK*AKP1-ONE )
262 A( K-1, K-1 ) = AKP1 / D
263 A( K, K ) = AK / D
264 A( K, K-1 ) = -AKKP1 / D
265 *
266 * Compute columns K-1 and K of the inverse.
267 *
268 IF( K.LT.N ) THEN
269 CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
270 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
271 $ ZERO, A( K+1, K ), 1 )
272 A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
273 $ 1 )
274 A( K, K-1 ) = A( K, K-1 ) -
275 $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
276 $ 1 )
277 CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
278 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
279 $ ZERO, A( K+1, K-1 ), 1 )
280 A( K-1, K-1 ) = A( K-1, K-1 ) -
281 $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
282 END IF
283 KSTEP = 2
284 END IF
285 *
286 KP = ABS( IPIV( K ) )
287 IF( KP.NE.K ) THEN
288 *
289 * Interchange rows and columns K and KP in the trailing
290 * submatrix A(k-1:n,k-1:n)
291 *
292 IF( KP.LT.N )
293 $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
294 CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
295 TEMP = A( K, K )
296 A( K, K ) = A( KP, KP )
297 A( KP, KP ) = TEMP
298 IF( KSTEP.EQ.2 ) THEN
299 TEMP = A( K, K-1 )
300 A( K, K-1 ) = A( KP, K-1 )
301 A( KP, K-1 ) = TEMP
302 END IF
303 END IF
304 *
305 K = K - KSTEP
306 GO TO 50
307 60 CONTINUE
308 END IF
309 *
310 RETURN
311 *
312 * End of ZSYTRI
313 *
314 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 COMPLEX*16 A( LDA, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZSYTRI computes the inverse of a complex symmetric indefinite matrix
21 * A using the factorization A = U*D*U**T or A = L*D*L**T computed by
22 * ZSYTRF.
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) COMPLEX*16 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 ZSYTRF.
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 ZSYTRF.
53 *
54 * WORK (workspace) COMPLEX*16 array, dimension (2*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 COMPLEX*16 ONE, ZERO
66 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
67 $ ZERO = ( 0.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 LOGICAL UPPER
71 INTEGER K, KP, KSTEP
72 COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
73 * ..
74 * .. External Functions ..
75 LOGICAL LSAME
76 COMPLEX*16 ZDOTU
77 EXTERNAL LSAME, ZDOTU
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA, ZCOPY, ZSWAP, ZSYMV
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC ABS, MAX
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
96 INFO = -4
97 END IF
98 IF( INFO.NE.0 ) THEN
99 CALL XERBLA( 'ZSYTRI', -INFO )
100 RETURN
101 END IF
102 *
103 * Quick return if possible
104 *
105 IF( N.EQ.0 )
106 $ RETURN
107 *
108 * Check that the diagonal matrix D is nonsingular.
109 *
110 IF( UPPER ) THEN
111 *
112 * Upper triangular storage: examine D from bottom to top
113 *
114 DO 10 INFO = N, 1, -1
115 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
116 $ RETURN
117 10 CONTINUE
118 ELSE
119 *
120 * Lower triangular storage: examine D from top to bottom.
121 *
122 DO 20 INFO = 1, N
123 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
124 $ RETURN
125 20 CONTINUE
126 END IF
127 INFO = 0
128 *
129 IF( UPPER ) THEN
130 *
131 * Compute inv(A) from the factorization A = U*D*U**T.
132 *
133 * K is the main loop index, increasing from 1 to N in steps of
134 * 1 or 2, depending on the size of the diagonal blocks.
135 *
136 K = 1
137 30 CONTINUE
138 *
139 * If K > N, exit from loop.
140 *
141 IF( K.GT.N )
142 $ GO TO 40
143 *
144 IF( IPIV( K ).GT.0 ) THEN
145 *
146 * 1 x 1 diagonal block
147 *
148 * Invert the diagonal block.
149 *
150 A( K, K ) = ONE / A( K, K )
151 *
152 * Compute column K of the inverse.
153 *
154 IF( K.GT.1 ) THEN
155 CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
156 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
157 $ A( 1, K ), 1 )
158 A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
159 $ 1 )
160 END IF
161 KSTEP = 1
162 ELSE
163 *
164 * 2 x 2 diagonal block
165 *
166 * Invert the diagonal block.
167 *
168 T = A( K, K+1 )
169 AK = A( K, K ) / T
170 AKP1 = A( K+1, K+1 ) / T
171 AKKP1 = A( K, K+1 ) / T
172 D = T*( AK*AKP1-ONE )
173 A( K, K ) = AKP1 / D
174 A( K+1, K+1 ) = AK / D
175 A( K, K+1 ) = -AKKP1 / D
176 *
177 * Compute columns K and K+1 of the inverse.
178 *
179 IF( K.GT.1 ) THEN
180 CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
181 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
182 $ A( 1, K ), 1 )
183 A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
184 $ 1 )
185 A( K, K+1 ) = A( K, K+1 ) -
186 $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
187 CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
188 CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
189 $ A( 1, K+1 ), 1 )
190 A( K+1, K+1 ) = A( K+1, K+1 ) -
191 $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
192 END IF
193 KSTEP = 2
194 END IF
195 *
196 KP = ABS( IPIV( K ) )
197 IF( KP.NE.K ) THEN
198 *
199 * Interchange rows and columns K and KP in the leading
200 * submatrix A(1:k+1,1:k+1)
201 *
202 CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
203 CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
204 TEMP = A( K, K )
205 A( K, K ) = A( KP, KP )
206 A( KP, KP ) = TEMP
207 IF( KSTEP.EQ.2 ) THEN
208 TEMP = A( K, K+1 )
209 A( K, K+1 ) = A( KP, K+1 )
210 A( KP, K+1 ) = TEMP
211 END IF
212 END IF
213 *
214 K = K + KSTEP
215 GO TO 30
216 40 CONTINUE
217 *
218 ELSE
219 *
220 * Compute inv(A) from the factorization A = L*D*L**T.
221 *
222 * K is the main loop index, increasing from 1 to N in steps of
223 * 1 or 2, depending on the size of the diagonal blocks.
224 *
225 K = N
226 50 CONTINUE
227 *
228 * If K < 1, exit from loop.
229 *
230 IF( K.LT.1 )
231 $ GO TO 60
232 *
233 IF( IPIV( K ).GT.0 ) THEN
234 *
235 * 1 x 1 diagonal block
236 *
237 * Invert the diagonal block.
238 *
239 A( K, K ) = ONE / A( K, K )
240 *
241 * Compute column K of the inverse.
242 *
243 IF( K.LT.N ) THEN
244 CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
245 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
246 $ ZERO, A( K+1, K ), 1 )
247 A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
248 $ 1 )
249 END IF
250 KSTEP = 1
251 ELSE
252 *
253 * 2 x 2 diagonal block
254 *
255 * Invert the diagonal block.
256 *
257 T = A( K, K-1 )
258 AK = A( K-1, K-1 ) / T
259 AKP1 = A( K, K ) / T
260 AKKP1 = A( K, K-1 ) / T
261 D = T*( AK*AKP1-ONE )
262 A( K-1, K-1 ) = AKP1 / D
263 A( K, K ) = AK / D
264 A( K, K-1 ) = -AKKP1 / D
265 *
266 * Compute columns K-1 and K of the inverse.
267 *
268 IF( K.LT.N ) THEN
269 CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
270 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
271 $ ZERO, A( K+1, K ), 1 )
272 A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
273 $ 1 )
274 A( K, K-1 ) = A( K, K-1 ) -
275 $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
276 $ 1 )
277 CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
278 CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
279 $ ZERO, A( K+1, K-1 ), 1 )
280 A( K-1, K-1 ) = A( K-1, K-1 ) -
281 $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
282 END IF
283 KSTEP = 2
284 END IF
285 *
286 KP = ABS( IPIV( K ) )
287 IF( KP.NE.K ) THEN
288 *
289 * Interchange rows and columns K and KP in the trailing
290 * submatrix A(k-1:n,k-1:n)
291 *
292 IF( KP.LT.N )
293 $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
294 CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
295 TEMP = A( K, K )
296 A( K, K ) = A( KP, KP )
297 A( KP, KP ) = TEMP
298 IF( KSTEP.EQ.2 ) THEN
299 TEMP = A( K, K-1 )
300 A( K, K-1 ) = A( KP, K-1 )
301 A( KP, K-1 ) = TEMP
302 END IF
303 END IF
304 *
305 K = K - KSTEP
306 GO TO 50
307 60 CONTINUE
308 END IF
309 *
310 RETURN
311 *
312 * End of ZSYTRI
313 *
314 END