LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dsytri ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK,
integer  INFO 
)

DSYTRI

Download DSYTRI + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 116 of file dsytri.f.

116 *
117 * -- LAPACK computational routine (version 3.4.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2011
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 *
199 * Compute inv(A) from the factorization A = U*D*U**T.
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 *
288 * Compute inv(A) from the factorization A = L*D*L**T.
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 *
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dsymv
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:154

Here is the call graph for this function:

Here is the caller graph for this function: