Annotation of rpl/lapack/lapack/zhetri.f, revision 1.9
1.9 ! bertrand 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 November 2011
! 111: *
! 112: *> \ingroup complex16HEcomputational
! 113: *
! 114: * =====================================================================
1.1 bertrand 115: SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
116: *
1.9 ! bertrand 117: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 118: * -- LAPACK is a software package provided by Univ. of Tennessee, --
119: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 120: * November 2011
1.1 bertrand 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: *
1.8 bertrand 202: * Compute inv(A) from the factorization A = U*D*U**H.
1.1 bertrand 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: *
1.8 bertrand 297: * Compute inv(A) from the factorization A = L*D*L**H.
1.1 bertrand 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