Annotation of rpl/lapack/lapack/zsytri.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZSYTRI
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZSYTRI + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZSYTRI( 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: *> ZSYTRI computes the inverse of a complex symmetric indefinite matrix
! 39: *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
! 40: *> ZSYTRF.
! 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**T;
! 52: *> = 'L': Lower triangular, form is A = L*D*L**T.
! 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 ZSYTRF.
! 66: *>
! 67: *> On exit, if INFO = 0, the (symmetric) 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 ZSYTRF.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[out] WORK
! 89: *> \verbatim
! 90: *> WORK is COMPLEX*16 array, dimension (2*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 complex16SYcomputational
! 113: *
! 114: * =====================================================================
1.1 bertrand 115: SUBROUTINE ZSYTRI( 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: COMPLEX*16 ONE, ZERO
135: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
136: $ ZERO = ( 0.0D+0, 0.0D+0 ) )
137: * ..
138: * .. Local Scalars ..
139: LOGICAL UPPER
140: INTEGER K, KP, KSTEP
141: COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
142: * ..
143: * .. External Functions ..
144: LOGICAL LSAME
145: COMPLEX*16 ZDOTU
146: EXTERNAL LSAME, ZDOTU
147: * ..
148: * .. External Subroutines ..
149: EXTERNAL XERBLA, ZCOPY, ZSWAP, ZSYMV
150: * ..
151: * .. Intrinsic Functions ..
152: INTRINSIC ABS, MAX
153: * ..
154: * .. Executable Statements ..
155: *
156: * Test the input parameters.
157: *
158: INFO = 0
159: UPPER = LSAME( UPLO, 'U' )
160: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
161: INFO = -1
162: ELSE IF( N.LT.0 ) THEN
163: INFO = -2
164: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
165: INFO = -4
166: END IF
167: IF( INFO.NE.0 ) THEN
168: CALL XERBLA( 'ZSYTRI', -INFO )
169: RETURN
170: END IF
171: *
172: * Quick return if possible
173: *
174: IF( N.EQ.0 )
175: $ RETURN
176: *
177: * Check that the diagonal matrix D is nonsingular.
178: *
179: IF( UPPER ) THEN
180: *
181: * Upper triangular storage: examine D from bottom to top
182: *
183: DO 10 INFO = N, 1, -1
184: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
185: $ RETURN
186: 10 CONTINUE
187: ELSE
188: *
189: * Lower triangular storage: examine D from top to bottom.
190: *
191: DO 20 INFO = 1, N
192: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
193: $ RETURN
194: 20 CONTINUE
195: END IF
196: INFO = 0
197: *
198: IF( UPPER ) THEN
199: *
1.8 bertrand 200: * Compute inv(A) from the factorization A = U*D*U**T.
1.1 bertrand 201: *
202: * K is the main loop index, increasing from 1 to N in steps of
203: * 1 or 2, depending on the size of the diagonal blocks.
204: *
205: K = 1
206: 30 CONTINUE
207: *
208: * If K > N, exit from loop.
209: *
210: IF( K.GT.N )
211: $ GO TO 40
212: *
213: IF( IPIV( K ).GT.0 ) THEN
214: *
215: * 1 x 1 diagonal block
216: *
217: * Invert the diagonal block.
218: *
219: A( K, K ) = ONE / A( K, K )
220: *
221: * Compute column K of the inverse.
222: *
223: IF( K.GT.1 ) THEN
224: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
225: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
226: $ A( 1, K ), 1 )
227: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
228: $ 1 )
229: END IF
230: KSTEP = 1
231: ELSE
232: *
233: * 2 x 2 diagonal block
234: *
235: * Invert the diagonal block.
236: *
237: T = A( K, K+1 )
238: AK = A( K, K ) / T
239: AKP1 = A( K+1, K+1 ) / T
240: AKKP1 = A( K, K+1 ) / T
241: D = T*( AK*AKP1-ONE )
242: A( K, K ) = AKP1 / D
243: A( K+1, K+1 ) = AK / D
244: A( K, K+1 ) = -AKKP1 / D
245: *
246: * Compute columns K and K+1 of the inverse.
247: *
248: IF( K.GT.1 ) THEN
249: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
250: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
251: $ A( 1, K ), 1 )
252: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
253: $ 1 )
254: A( K, K+1 ) = A( K, K+1 ) -
255: $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
256: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
257: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
258: $ A( 1, K+1 ), 1 )
259: A( K+1, K+1 ) = A( K+1, K+1 ) -
260: $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 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: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
273: TEMP = A( K, K )
274: A( K, K ) = A( KP, KP )
275: A( KP, KP ) = TEMP
276: IF( KSTEP.EQ.2 ) THEN
277: TEMP = A( K, K+1 )
278: A( K, K+1 ) = A( KP, K+1 )
279: A( KP, K+1 ) = TEMP
280: END IF
281: END IF
282: *
283: K = K + KSTEP
284: GO TO 30
285: 40 CONTINUE
286: *
287: ELSE
288: *
1.8 bertrand 289: * Compute inv(A) from the factorization A = L*D*L**T.
1.1 bertrand 290: *
291: * K is the main loop index, increasing from 1 to N in steps of
292: * 1 or 2, depending on the size of the diagonal blocks.
293: *
294: K = N
295: 50 CONTINUE
296: *
297: * If K < 1, exit from loop.
298: *
299: IF( K.LT.1 )
300: $ GO TO 60
301: *
302: IF( IPIV( K ).GT.0 ) THEN
303: *
304: * 1 x 1 diagonal block
305: *
306: * Invert the diagonal block.
307: *
308: A( K, K ) = ONE / A( K, K )
309: *
310: * Compute column K of the inverse.
311: *
312: IF( K.LT.N ) THEN
313: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
314: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
315: $ ZERO, A( K+1, K ), 1 )
316: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
317: $ 1 )
318: END IF
319: KSTEP = 1
320: ELSE
321: *
322: * 2 x 2 diagonal block
323: *
324: * Invert the diagonal block.
325: *
326: T = A( K, K-1 )
327: AK = A( K-1, K-1 ) / T
328: AKP1 = A( K, K ) / T
329: AKKP1 = A( K, K-1 ) / T
330: D = T*( AK*AKP1-ONE )
331: A( K-1, K-1 ) = AKP1 / D
332: A( K, K ) = AK / D
333: A( K, K-1 ) = -AKKP1 / D
334: *
335: * Compute columns K-1 and K of the inverse.
336: *
337: IF( K.LT.N ) THEN
338: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
339: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
340: $ ZERO, A( K+1, K ), 1 )
341: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
342: $ 1 )
343: A( K, K-1 ) = A( K, K-1 ) -
344: $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
345: $ 1 )
346: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
347: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
348: $ ZERO, A( K+1, K-1 ), 1 )
349: A( K-1, K-1 ) = A( K-1, K-1 ) -
350: $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
351: END IF
352: KSTEP = 2
353: END IF
354: *
355: KP = ABS( IPIV( K ) )
356: IF( KP.NE.K ) THEN
357: *
358: * Interchange rows and columns K and KP in the trailing
359: * submatrix A(k-1:n,k-1:n)
360: *
361: IF( KP.LT.N )
362: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
363: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
364: TEMP = A( K, K )
365: A( K, K ) = A( KP, KP )
366: A( KP, KP ) = TEMP
367: IF( KSTEP.EQ.2 ) THEN
368: TEMP = A( K, K-1 )
369: A( K, K-1 ) = A( KP, K-1 )
370: A( KP, K-1 ) = TEMP
371: END IF
372: END IF
373: *
374: K = K - KSTEP
375: GO TO 50
376: 60 CONTINUE
377: END IF
378: *
379: RETURN
380: *
381: * End of ZSYTRI
382: *
383: 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