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