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: *> \date December 2016
111: *
112: *> \ingroup complex16HEcomputational
113: *
114: * =====================================================================
115: SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
116: *
117: * -- LAPACK computational routine (version 3.7.0) --
118: * -- LAPACK is a software package provided by Univ. of Tennessee, --
119: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120: * December 2016
121: *
122: * .. Scalar Arguments ..
123: CHARACTER UPLO
124: INTEGER INFO, LDA, N
125: * ..
126: * .. Array Arguments ..
127: INTEGER IPIV( * )
128: COMPLEX*16 A( LDA, * ), WORK( * )
129: * ..
130: *
131: * =====================================================================
132: *
133: * .. Parameters ..
134: DOUBLE PRECISION ONE
135: COMPLEX*16 CONE, ZERO
136: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
137: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
138: * ..
139: * .. Local Scalars ..
140: LOGICAL UPPER
141: INTEGER J, K, KP, KSTEP
142: DOUBLE PRECISION AK, AKP1, D, T
143: COMPLEX*16 AKKP1, TEMP
144: * ..
145: * .. External Functions ..
146: LOGICAL LSAME
147: COMPLEX*16 ZDOTC
148: EXTERNAL LSAME, ZDOTC
149: * ..
150: * .. External Subroutines ..
151: EXTERNAL XERBLA, ZCOPY, ZHEMV, ZSWAP
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC ABS, DBLE, DCONJG, MAX
155: * ..
156: * .. Executable Statements ..
157: *
158: * Test the input parameters.
159: *
160: INFO = 0
161: UPPER = LSAME( UPLO, 'U' )
162: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
163: INFO = -1
164: ELSE IF( N.LT.0 ) THEN
165: INFO = -2
166: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
167: INFO = -4
168: END IF
169: IF( INFO.NE.0 ) THEN
170: CALL XERBLA( 'ZHETRI', -INFO )
171: RETURN
172: END IF
173: *
174: * Quick return if possible
175: *
176: IF( N.EQ.0 )
177: $ RETURN
178: *
179: * Check that the diagonal matrix D is nonsingular.
180: *
181: IF( UPPER ) THEN
182: *
183: * Upper triangular storage: examine D from bottom to top
184: *
185: DO 10 INFO = N, 1, -1
186: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
187: $ RETURN
188: 10 CONTINUE
189: ELSE
190: *
191: * Lower triangular storage: examine D from top to bottom.
192: *
193: DO 20 INFO = 1, N
194: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
195: $ RETURN
196: 20 CONTINUE
197: END IF
198: INFO = 0
199: *
200: IF( UPPER ) THEN
201: *
202: * Compute inv(A) from the factorization A = U*D*U**H.
203: *
204: * K is the main loop index, increasing from 1 to N in steps of
205: * 1 or 2, depending on the size of the diagonal blocks.
206: *
207: K = 1
208: 30 CONTINUE
209: *
210: * If K > N, exit from loop.
211: *
212: IF( K.GT.N )
213: $ GO TO 50
214: *
215: IF( IPIV( K ).GT.0 ) THEN
216: *
217: * 1 x 1 diagonal block
218: *
219: * Invert the diagonal block.
220: *
221: A( K, K ) = ONE / DBLE( A( K, K ) )
222: *
223: * Compute column K of the inverse.
224: *
225: IF( K.GT.1 ) THEN
226: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
227: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
228: $ A( 1, K ), 1 )
229: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
230: $ K ), 1 ) )
231: END IF
232: KSTEP = 1
233: ELSE
234: *
235: * 2 x 2 diagonal block
236: *
237: * Invert the diagonal block.
238: *
239: T = ABS( A( K, K+1 ) )
240: AK = DBLE( A( K, K ) ) / T
241: AKP1 = DBLE( A( K+1, K+1 ) ) / T
242: AKKP1 = A( K, K+1 ) / T
243: D = T*( AK*AKP1-ONE )
244: A( K, K ) = AKP1 / D
245: A( K+1, K+1 ) = AK / D
246: A( K, K+1 ) = -AKKP1 / D
247: *
248: * Compute columns K and K+1 of the inverse.
249: *
250: IF( K.GT.1 ) THEN
251: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
252: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
253: $ A( 1, K ), 1 )
254: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
255: $ K ), 1 ) )
256: A( K, K+1 ) = A( K, K+1 ) -
257: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
258: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
259: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
260: $ A( 1, K+1 ), 1 )
261: A( K+1, K+1 ) = A( K+1, K+1 ) -
262: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
263: $ 1 ) )
264: END IF
265: KSTEP = 2
266: END IF
267: *
268: KP = ABS( IPIV( K ) )
269: IF( KP.NE.K ) THEN
270: *
271: * Interchange rows and columns K and KP in the leading
272: * submatrix A(1:k+1,1:k+1)
273: *
274: CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
275: DO 40 J = KP + 1, K - 1
276: TEMP = DCONJG( A( J, K ) )
277: A( J, K ) = DCONJG( A( KP, J ) )
278: A( KP, J ) = TEMP
279: 40 CONTINUE
280: A( KP, K ) = DCONJG( A( KP, K ) )
281: TEMP = A( K, K )
282: A( K, K ) = A( KP, KP )
283: A( KP, KP ) = TEMP
284: IF( KSTEP.EQ.2 ) THEN
285: TEMP = A( K, K+1 )
286: A( K, K+1 ) = A( KP, K+1 )
287: A( KP, K+1 ) = TEMP
288: END IF
289: END IF
290: *
291: K = K + KSTEP
292: GO TO 30
293: 50 CONTINUE
294: *
295: ELSE
296: *
297: * Compute inv(A) from the factorization A = L*D*L**H.
298: *
299: * K is the main loop index, increasing from 1 to N in steps of
300: * 1 or 2, depending on the size of the diagonal blocks.
301: *
302: K = N
303: 60 CONTINUE
304: *
305: * If K < 1, exit from loop.
306: *
307: IF( K.LT.1 )
308: $ GO TO 80
309: *
310: IF( IPIV( K ).GT.0 ) THEN
311: *
312: * 1 x 1 diagonal block
313: *
314: * Invert the diagonal block.
315: *
316: A( K, K ) = ONE / DBLE( A( K, K ) )
317: *
318: * Compute column K of the inverse.
319: *
320: IF( K.LT.N ) THEN
321: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
322: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
323: $ 1, ZERO, A( K+1, K ), 1 )
324: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
325: $ A( K+1, K ), 1 ) )
326: END IF
327: KSTEP = 1
328: ELSE
329: *
330: * 2 x 2 diagonal block
331: *
332: * Invert the diagonal block.
333: *
334: T = ABS( A( K, K-1 ) )
335: AK = DBLE( A( K-1, K-1 ) ) / T
336: AKP1 = DBLE( A( K, K ) ) / T
337: AKKP1 = A( K, K-1 ) / T
338: D = T*( AK*AKP1-ONE )
339: A( K-1, K-1 ) = AKP1 / D
340: A( K, K ) = AK / D
341: A( K, K-1 ) = -AKKP1 / D
342: *
343: * Compute columns K-1 and K of the inverse.
344: *
345: IF( K.LT.N ) THEN
346: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
347: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
348: $ 1, ZERO, A( K+1, K ), 1 )
349: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
350: $ A( K+1, K ), 1 ) )
351: A( K, K-1 ) = A( K, K-1 ) -
352: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
353: $ 1 )
354: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
355: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
356: $ 1, ZERO, A( K+1, K-1 ), 1 )
357: A( K-1, K-1 ) = A( K-1, K-1 ) -
358: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
359: $ 1 ) )
360: END IF
361: KSTEP = 2
362: END IF
363: *
364: KP = ABS( IPIV( K ) )
365: IF( KP.NE.K ) THEN
366: *
367: * Interchange rows and columns K and KP in the trailing
368: * submatrix A(k-1:n,k-1:n)
369: *
370: IF( KP.LT.N )
371: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
372: DO 70 J = K + 1, KP - 1
373: TEMP = DCONJG( A( J, K ) )
374: A( J, K ) = DCONJG( A( KP, J ) )
375: A( KP, J ) = TEMP
376: 70 CONTINUE
377: A( KP, K ) = DCONJG( A( KP, K ) )
378: TEMP = A( K, K )
379: A( K, K ) = A( KP, KP )
380: A( KP, KP ) = TEMP
381: IF( KSTEP.EQ.2 ) THEN
382: TEMP = A( K, K-1 )
383: A( K, K-1 ) = A( KP, K-1 )
384: A( KP, K-1 ) = TEMP
385: END IF
386: END IF
387: *
388: K = K - KSTEP
389: GO TO 60
390: 80 CONTINUE
391: END IF
392: *
393: RETURN
394: *
395: * End of ZHETRI
396: *
397: 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