LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dsytri()

 subroutine dsytri ( character UPLO, integer N, double precision, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, double precision, dimension( * ) WORK, integer INFO )

DSYTRI

Purpose:
``` DSYTRI computes the inverse of a real symmetric indefinite matrix
A using the factorization A = U*D*U**T or A = L*D*L**T computed by
DSYTRF.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the details of the factorization are stored as an upper or lower triangular matrix. = 'U': Upper triangular, form is A = U*D*U**T; = 'L': Lower triangular, form is A = L*D*L**T.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is DOUBLE PRECISION array, dimension (LDA,N) On entry, the block diagonal matrix D and the multipliers used to obtain the factor U or L as computed by DSYTRF. On exit, if INFO = 0, the (symmetric) inverse of the original matrix. If UPLO = 'U', the upper triangular part of the inverse is formed and the part of A below the diagonal is not referenced; if UPLO = 'L' the lower triangular part of the inverse is formed and the part of A above the diagonal is not referenced.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in] IPIV ``` IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D as determined by DSYTRF.``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its inverse could not be computed.```

Definition at line 113 of file dsytri.f.

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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:152
Here is the call graph for this function:
Here is the caller graph for this function: