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

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