LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
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

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

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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 126 of file dsytri_rook.f.

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