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