LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
April 2012
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 131 of file dsytri_rook.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: