1: SUBROUTINE ZSYTRI( 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: 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'.
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'.
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytri.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack