LAPACK 3.12.0
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 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*
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: