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