 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dsytri_rook()

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

DSYTRI_ROOK

Purpose:
``` DSYTRI_ROOK computes the inverse of a real symmetric
matrix A using the factorization A = U*D*U**T or A = L*D*L**T
computed by DSYTRF_ROOK.```
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_ROOK. 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_ROOK.``` [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.```
Contributors:
```   April 2012, Igor Kozachenko,
Computer Science Division,
University of California, Berkeley

September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
School of Mathematics,
University of Manchester```

Definition at line 128 of file dsytri_rook.f.

129 *
130 * -- LAPACK computational routine --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 *
134 * .. Scalar Arguments ..
135  CHARACTER UPLO
136  INTEGER INFO, LDA, N
137 * ..
138 * .. Array Arguments ..
139  INTEGER IPIV( * )
140  DOUBLE PRECISION A( LDA, * ), WORK( * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION ONE, ZERO
147  parameter( one = 1.0d+0, zero = 0.0d+0 )
148 * ..
149 * .. Local Scalars ..
150  LOGICAL UPPER
151  INTEGER K, KP, KSTEP
152  DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
153 * ..
154 * .. External Functions ..
155  LOGICAL LSAME
156  DOUBLE PRECISION DDOT
157  EXTERNAL lsame, ddot
158 * ..
159 * .. External Subroutines ..
160  EXTERNAL dcopy, dswap, dsymv, xerbla
161 * ..
162 * .. Intrinsic Functions ..
163  INTRINSIC abs, max
164 * ..
165 * .. Executable Statements ..
166 *
167 * Test the input parameters.
168 *
169  info = 0
170  upper = lsame( uplo, 'U' )
171  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
172  info = -1
173  ELSE IF( n.LT.0 ) THEN
174  info = -2
175  ELSE IF( lda.LT.max( 1, n ) ) THEN
176  info = -4
177  END IF
178  IF( info.NE.0 ) THEN
179  CALL xerbla( 'DSYTRI_ROOK', -info )
180  RETURN
181  END IF
182 *
183 * Quick return if possible
184 *
185  IF( n.EQ.0 )
186  \$ RETURN
187 *
188 * Check that the diagonal matrix D is nonsingular.
189 *
190  IF( upper ) THEN
191 *
192 * Upper triangular storage: examine D from bottom to top
193 *
194  DO 10 info = n, 1, -1
195  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
196  \$ RETURN
197  10 CONTINUE
198  ELSE
199 *
200 * Lower triangular storage: examine D from top to bottom.
201 *
202  DO 20 info = 1, n
203  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
204  \$ RETURN
205  20 CONTINUE
206  END IF
207  info = 0
208 *
209  IF( upper ) THEN
210 *
211 * Compute inv(A) from the factorization A = U*D*U**T.
212 *
213 * K is the main loop index, increasing from 1 to N in steps of
214 * 1 or 2, depending on the size of the diagonal blocks.
215 *
216  k = 1
217  30 CONTINUE
218 *
219 * If K > N, exit from loop.
220 *
221  IF( k.GT.n )
222  \$ GO TO 40
223 *
224  IF( ipiv( k ).GT.0 ) THEN
225 *
226 * 1 x 1 diagonal block
227 *
228 * Invert the diagonal block.
229 *
230  a( k, k ) = one / a( k, k )
231 *
232 * Compute column K of the inverse.
233 *
234  IF( k.GT.1 ) THEN
235  CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
236  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
237  \$ a( 1, k ), 1 )
238  a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
239  \$ 1 )
240  END IF
241  kstep = 1
242  ELSE
243 *
244 * 2 x 2 diagonal block
245 *
246 * Invert the diagonal block.
247 *
248  t = abs( a( k, k+1 ) )
249  ak = a( k, k ) / t
250  akp1 = a( k+1, k+1 ) / t
251  akkp1 = a( k, k+1 ) / t
252  d = t*( ak*akp1-one )
253  a( k, k ) = akp1 / d
254  a( k+1, k+1 ) = ak / d
255  a( k, k+1 ) = -akkp1 / d
256 *
257 * Compute columns K and K+1 of the inverse.
258 *
259  IF( k.GT.1 ) THEN
260  CALL dcopy( k-1, a( 1, k ), 1, work, 1 )
261  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
262  \$ a( 1, k ), 1 )
263  a( k, k ) = a( k, k ) - ddot( k-1, work, 1, a( 1, k ),
264  \$ 1 )
265  a( k, k+1 ) = a( k, k+1 ) -
266  \$ ddot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
267  CALL dcopy( k-1, a( 1, k+1 ), 1, work, 1 )
268  CALL dsymv( uplo, k-1, -one, a, lda, work, 1, zero,
269  \$ a( 1, k+1 ), 1 )
270  a( k+1, k+1 ) = a( k+1, k+1 ) -
271  \$ ddot( k-1, work, 1, a( 1, k+1 ), 1 )
272  END IF
273  kstep = 2
274  END IF
275 *
276  IF( kstep.EQ.1 ) THEN
277 *
278 * Interchange rows and columns K and IPIV(K) in the leading
279 * submatrix A(1:k+1,1:k+1)
280 *
281  kp = ipiv( k )
282  IF( kp.NE.k ) THEN
283  IF( kp.GT.1 )
284  \$ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
285  CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
286  temp = a( k, k )
287  a( k, k ) = a( kp, kp )
288  a( kp, kp ) = temp
289  END IF
290  ELSE
291 *
292 * Interchange rows and columns K and K+1 with -IPIV(K) and
293 * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
294 *
295  kp = -ipiv( k )
296  IF( kp.NE.k ) THEN
297  IF( kp.GT.1 )
298  \$ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
299  CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
300 *
301  temp = a( k, k )
302  a( k, k ) = a( kp, kp )
303  a( kp, kp ) = temp
304  temp = a( k, k+1 )
305  a( k, k+1 ) = a( kp, k+1 )
306  a( kp, k+1 ) = temp
307  END IF
308 *
309  k = k + 1
310  kp = -ipiv( k )
311  IF( kp.NE.k ) THEN
312  IF( kp.GT.1 )
313  \$ CALL dswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
314  CALL dswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
315  temp = a( k, k )
316  a( k, k ) = a( kp, kp )
317  a( kp, kp ) = temp
318  END IF
319  END IF
320 *
321  k = k + 1
322  GO TO 30
323  40 CONTINUE
324 *
325  ELSE
326 *
327 * Compute inv(A) from the factorization A = L*D*L**T.
328 *
329 * K is the main loop index, increasing from 1 to N in steps of
330 * 1 or 2, depending on the size of the diagonal blocks.
331 *
332  k = n
333  50 CONTINUE
334 *
335 * If K < 1, exit from loop.
336 *
337  IF( k.LT.1 )
338  \$ GO TO 60
339 *
340  IF( ipiv( k ).GT.0 ) THEN
341 *
342 * 1 x 1 diagonal block
343 *
344 * Invert the diagonal block.
345 *
346  a( k, k ) = one / a( k, k )
347 *
348 * Compute column K of the inverse.
349 *
350  IF( k.LT.n ) THEN
351  CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
352  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
353  \$ zero, a( k+1, k ), 1 )
354  a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
355  \$ 1 )
356  END IF
357  kstep = 1
358  ELSE
359 *
360 * 2 x 2 diagonal block
361 *
362 * Invert the diagonal block.
363 *
364  t = abs( a( k, k-1 ) )
365  ak = a( k-1, k-1 ) / t
366  akp1 = a( k, k ) / t
367  akkp1 = a( k, k-1 ) / t
368  d = t*( ak*akp1-one )
369  a( k-1, k-1 ) = akp1 / d
370  a( k, k ) = ak / d
371  a( k, k-1 ) = -akkp1 / d
372 *
373 * Compute columns K-1 and K of the inverse.
374 *
375  IF( k.LT.n ) THEN
376  CALL dcopy( n-k, a( k+1, k ), 1, work, 1 )
377  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
378  \$ zero, a( k+1, k ), 1 )
379  a( k, k ) = a( k, k ) - ddot( n-k, work, 1, a( k+1, k ),
380  \$ 1 )
381  a( k, k-1 ) = a( k, k-1 ) -
382  \$ ddot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
383  \$ 1 )
384  CALL dcopy( n-k, a( k+1, k-1 ), 1, work, 1 )
385  CALL dsymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
386  \$ zero, a( k+1, k-1 ), 1 )
387  a( k-1, k-1 ) = a( k-1, k-1 ) -
388  \$ ddot( n-k, work, 1, a( k+1, k-1 ), 1 )
389  END IF
390  kstep = 2
391  END IF
392 *
393  IF( kstep.EQ.1 ) THEN
394 *
395 * Interchange rows and columns K and IPIV(K) in the trailing
396 * submatrix A(k-1:n,k-1:n)
397 *
398  kp = ipiv( k )
399  IF( kp.NE.k ) THEN
400  IF( kp.LT.n )
401  \$ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
402  CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
403  temp = a( k, k )
404  a( k, k ) = a( kp, kp )
405  a( kp, kp ) = temp
406  END IF
407  ELSE
408 *
409 * Interchange rows and columns K and K-1 with -IPIV(K) and
410 * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
411 *
412  kp = -ipiv( k )
413  IF( kp.NE.k ) THEN
414  IF( kp.LT.n )
415  \$ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
416  CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
417 *
418  temp = a( k, k )
419  a( k, k ) = a( kp, kp )
420  a( kp, kp ) = temp
421  temp = a( k, k-1 )
422  a( k, k-1 ) = a( kp, k-1 )
423  a( kp, k-1 ) = temp
424  END IF
425 *
426  k = k - 1
427  kp = -ipiv( k )
428  IF( kp.NE.k ) THEN
429  IF( kp.LT.n )
430  \$ CALL dswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
431  CALL dswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
432  temp = a( k, k )
433  a( k, k ) = a( kp, kp )
434  a( kp, kp ) = temp
435  END IF
436  END IF
437 *
438  k = k - 1
439  GO TO 50
440  60 CONTINUE
441  END IF
442 *
443  RETURN
444 *
445 * End of DSYTRI_ROOK
446 *
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: