LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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)
Definition cblat2.f:3285
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 dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: