1: *> \brief \b ZHETRI
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZHETRI + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHETRI computes the inverse of a complex Hermitian indefinite matrix
39: *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
40: *> ZHETRF.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] UPLO
47: *> \verbatim
48: *> UPLO is CHARACTER*1
49: *> Specifies whether the details of the factorization are stored
50: *> as an upper or lower triangular matrix.
51: *> = 'U': Upper triangular, form is A = U*D*U**H;
52: *> = 'L': Lower triangular, form is A = L*D*L**H.
53: *> \endverbatim
54: *>
55: *> \param[in] N
56: *> \verbatim
57: *> N is INTEGER
58: *> The order of the matrix A. N >= 0.
59: *> \endverbatim
60: *>
61: *> \param[in,out] A
62: *> \verbatim
63: *> A is COMPLEX*16 array, dimension (LDA,N)
64: *> On entry, the block diagonal matrix D and the multipliers
65: *> used to obtain the factor U or L as computed by ZHETRF.
66: *>
67: *> On exit, if INFO = 0, the (Hermitian) inverse of the original
68: *> matrix. If UPLO = 'U', the upper triangular part of the
69: *> inverse is formed and the part of A below the diagonal is not
70: *> referenced; if UPLO = 'L' the lower triangular part of the
71: *> inverse is formed and the part of A above the diagonal is
72: *> not referenced.
73: *> \endverbatim
74: *>
75: *> \param[in] LDA
76: *> \verbatim
77: *> LDA is INTEGER
78: *> The leading dimension of the array A. LDA >= max(1,N).
79: *> \endverbatim
80: *>
81: *> \param[in] IPIV
82: *> \verbatim
83: *> IPIV is INTEGER array, dimension (N)
84: *> Details of the interchanges and the block structure of D
85: *> as determined by ZHETRF.
86: *> \endverbatim
87: *>
88: *> \param[out] WORK
89: *> \verbatim
90: *> WORK is COMPLEX*16 array, dimension (N)
91: *> \endverbatim
92: *>
93: *> \param[out] INFO
94: *> \verbatim
95: *> INFO is INTEGER
96: *> = 0: successful exit
97: *> < 0: if INFO = -i, the i-th argument had an illegal value
98: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99: *> inverse could not be computed.
100: *> \endverbatim
101: *
102: * Authors:
103: * ========
104: *
105: *> \author Univ. of Tennessee
106: *> \author Univ. of California Berkeley
107: *> \author Univ. of Colorado Denver
108: *> \author NAG Ltd.
109: *
110: *> \ingroup complex16HEcomputational
111: *
112: * =====================================================================
113: SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
114: *
115: * -- LAPACK computational routine --
116: * -- LAPACK is a software package provided by Univ. of Tennessee, --
117: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118: *
119: * .. Scalar Arguments ..
120: CHARACTER UPLO
121: INTEGER INFO, LDA, N
122: * ..
123: * .. Array Arguments ..
124: INTEGER IPIV( * )
125: COMPLEX*16 A( LDA, * ), WORK( * )
126: * ..
127: *
128: * =====================================================================
129: *
130: * .. Parameters ..
131: DOUBLE PRECISION ONE
132: COMPLEX*16 CONE, ZERO
133: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
134: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
135: * ..
136: * .. Local Scalars ..
137: LOGICAL UPPER
138: INTEGER J, K, KP, KSTEP
139: DOUBLE PRECISION AK, AKP1, D, T
140: COMPLEX*16 AKKP1, TEMP
141: * ..
142: * .. External Functions ..
143: LOGICAL LSAME
144: COMPLEX*16 ZDOTC
145: EXTERNAL LSAME, ZDOTC
146: * ..
147: * .. External Subroutines ..
148: EXTERNAL XERBLA, ZCOPY, ZHEMV, ZSWAP
149: * ..
150: * .. Intrinsic Functions ..
151: INTRINSIC ABS, DBLE, DCONJG, MAX
152: * ..
153: * .. Executable Statements ..
154: *
155: * Test the input parameters.
156: *
157: INFO = 0
158: UPPER = LSAME( UPLO, 'U' )
159: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
160: INFO = -1
161: ELSE IF( N.LT.0 ) THEN
162: INFO = -2
163: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
164: INFO = -4
165: END IF
166: IF( INFO.NE.0 ) THEN
167: CALL XERBLA( 'ZHETRI', -INFO )
168: RETURN
169: END IF
170: *
171: * Quick return if possible
172: *
173: IF( N.EQ.0 )
174: $ RETURN
175: *
176: * Check that the diagonal matrix D is nonsingular.
177: *
178: IF( UPPER ) THEN
179: *
180: * Upper triangular storage: examine D from bottom to top
181: *
182: DO 10 INFO = N, 1, -1
183: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
184: $ RETURN
185: 10 CONTINUE
186: ELSE
187: *
188: * Lower triangular storage: examine D from top to bottom.
189: *
190: DO 20 INFO = 1, N
191: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
192: $ RETURN
193: 20 CONTINUE
194: END IF
195: INFO = 0
196: *
197: IF( UPPER ) THEN
198: *
199: * Compute inv(A) from the factorization A = U*D*U**H.
200: *
201: * K is the main loop index, increasing from 1 to N in steps of
202: * 1 or 2, depending on the size of the diagonal blocks.
203: *
204: K = 1
205: 30 CONTINUE
206: *
207: * If K > N, exit from loop.
208: *
209: IF( K.GT.N )
210: $ GO TO 50
211: *
212: IF( IPIV( K ).GT.0 ) THEN
213: *
214: * 1 x 1 diagonal block
215: *
216: * Invert the diagonal block.
217: *
218: A( K, K ) = ONE / DBLE( A( K, K ) )
219: *
220: * Compute column K of the inverse.
221: *
222: IF( K.GT.1 ) THEN
223: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
224: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
225: $ A( 1, K ), 1 )
226: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
227: $ K ), 1 ) )
228: END IF
229: KSTEP = 1
230: ELSE
231: *
232: * 2 x 2 diagonal block
233: *
234: * Invert the diagonal block.
235: *
236: T = ABS( A( K, K+1 ) )
237: AK = DBLE( A( K, K ) ) / T
238: AKP1 = DBLE( A( K+1, K+1 ) ) / T
239: AKKP1 = A( K, K+1 ) / T
240: D = T*( AK*AKP1-ONE )
241: A( K, K ) = AKP1 / D
242: A( K+1, K+1 ) = AK / D
243: A( K, K+1 ) = -AKKP1 / D
244: *
245: * Compute columns K and K+1 of the inverse.
246: *
247: IF( K.GT.1 ) THEN
248: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
249: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
250: $ A( 1, K ), 1 )
251: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
252: $ K ), 1 ) )
253: A( K, K+1 ) = A( K, K+1 ) -
254: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
255: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
256: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
257: $ A( 1, K+1 ), 1 )
258: A( K+1, K+1 ) = A( K+1, K+1 ) -
259: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
260: $ 1 ) )
261: END IF
262: KSTEP = 2
263: END IF
264: *
265: KP = ABS( IPIV( K ) )
266: IF( KP.NE.K ) THEN
267: *
268: * Interchange rows and columns K and KP in the leading
269: * submatrix A(1:k+1,1:k+1)
270: *
271: CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
272: DO 40 J = KP + 1, K - 1
273: TEMP = DCONJG( A( J, K ) )
274: A( J, K ) = DCONJG( A( KP, J ) )
275: A( KP, J ) = TEMP
276: 40 CONTINUE
277: A( KP, K ) = DCONJG( A( KP, K ) )
278: TEMP = A( K, K )
279: A( K, K ) = A( KP, KP )
280: A( KP, KP ) = TEMP
281: IF( KSTEP.EQ.2 ) THEN
282: TEMP = A( K, K+1 )
283: A( K, K+1 ) = A( KP, K+1 )
284: A( KP, K+1 ) = TEMP
285: END IF
286: END IF
287: *
288: K = K + KSTEP
289: GO TO 30
290: 50 CONTINUE
291: *
292: ELSE
293: *
294: * Compute inv(A) from the factorization A = L*D*L**H.
295: *
296: * K is the main loop index, increasing from 1 to N in steps of
297: * 1 or 2, depending on the size of the diagonal blocks.
298: *
299: K = N
300: 60 CONTINUE
301: *
302: * If K < 1, exit from loop.
303: *
304: IF( K.LT.1 )
305: $ GO TO 80
306: *
307: IF( IPIV( K ).GT.0 ) THEN
308: *
309: * 1 x 1 diagonal block
310: *
311: * Invert the diagonal block.
312: *
313: A( K, K ) = ONE / DBLE( A( K, K ) )
314: *
315: * Compute column K of the inverse.
316: *
317: IF( K.LT.N ) THEN
318: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
319: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
320: $ 1, ZERO, A( K+1, K ), 1 )
321: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
322: $ A( K+1, K ), 1 ) )
323: END IF
324: KSTEP = 1
325: ELSE
326: *
327: * 2 x 2 diagonal block
328: *
329: * Invert the diagonal block.
330: *
331: T = ABS( A( K, K-1 ) )
332: AK = DBLE( A( K-1, K-1 ) ) / T
333: AKP1 = DBLE( A( K, K ) ) / T
334: AKKP1 = A( K, K-1 ) / T
335: D = T*( AK*AKP1-ONE )
336: A( K-1, K-1 ) = AKP1 / D
337: A( K, K ) = AK / D
338: A( K, K-1 ) = -AKKP1 / D
339: *
340: * Compute columns K-1 and K of the inverse.
341: *
342: IF( K.LT.N ) THEN
343: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
344: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
345: $ 1, ZERO, A( K+1, K ), 1 )
346: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
347: $ A( K+1, K ), 1 ) )
348: A( K, K-1 ) = A( K, K-1 ) -
349: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
350: $ 1 )
351: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
352: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
353: $ 1, ZERO, A( K+1, K-1 ), 1 )
354: A( K-1, K-1 ) = A( K-1, K-1 ) -
355: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
356: $ 1 ) )
357: END IF
358: KSTEP = 2
359: END IF
360: *
361: KP = ABS( IPIV( K ) )
362: IF( KP.NE.K ) THEN
363: *
364: * Interchange rows and columns K and KP in the trailing
365: * submatrix A(k-1:n,k-1:n)
366: *
367: IF( KP.LT.N )
368: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
369: DO 70 J = K + 1, KP - 1
370: TEMP = DCONJG( A( J, K ) )
371: A( J, K ) = DCONJG( A( KP, J ) )
372: A( KP, J ) = TEMP
373: 70 CONTINUE
374: A( KP, K ) = DCONJG( A( KP, K ) )
375: TEMP = A( K, K )
376: A( K, K ) = A( KP, KP )
377: A( KP, KP ) = TEMP
378: IF( KSTEP.EQ.2 ) THEN
379: TEMP = A( K, K-1 )
380: A( K, K-1 ) = A( KP, K-1 )
381: A( KP, K-1 ) = TEMP
382: END IF
383: END IF
384: *
385: K = K - KSTEP
386: GO TO 60
387: 80 CONTINUE
388: END IF
389: *
390: RETURN
391: *
392: * End of ZHETRI
393: *
394: 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