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