LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ssytri_rook ( character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
real, dimension( * )  WORK,
integer  INFO 
)

SSYTRI_ROOK

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

Purpose: 
 SSYTRI_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 SSYTRF_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 REAL 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 SSYTRF_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 SSYTRF_ROOK.
[out]WORK
          WORK is REAL 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 ssytri_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  REAL a( lda, * ), work( * )
144 * ..
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  REAL one, zero
150  parameter ( one = 1.0e+0, zero = 0.0e+0 )
151 * ..
152 * .. Local Scalars ..
153  LOGICAL upper
154  INTEGER k, kp, kstep
155  REAL ak, akkp1, akp1, d, t, temp
156 * ..
157 * .. External Functions ..
158  LOGICAL lsame
159  REAL sdot
160  EXTERNAL lsame, sdot
161 * ..
162 * .. External Subroutines ..
163  EXTERNAL scopy, sswap, ssymv, 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( 'SSYTRI_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 scopy( k-1, a( 1, k ), 1, work, 1 )
239  CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
240  $ a( 1, k ), 1 )
241  a( k, k ) = a( k, k ) - sdot( 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 scopy( k-1, a( 1, k ), 1, work, 1 )
264  CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
265  $ a( 1, k ), 1 )
266  a( k, k ) = a( k, k ) - sdot( k-1, work, 1, a( 1, k ),
267  $ 1 )
268  a( k, k+1 ) = a( k, k+1 ) -
269  $ sdot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
270  CALL scopy( k-1, a( 1, k+1 ), 1, work, 1 )
271  CALL ssymv( 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  $ sdot( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
288  CALL sswap( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
302  CALL sswap( 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 sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
317  CALL sswap( 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 scopy( n-k, a( k+1, k ), 1, work, 1 )
355  CALL ssymv( 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 ) - sdot( 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 scopy( n-k, a( k+1, k ), 1, work, 1 )
380  CALL ssymv( 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 ) - sdot( n-k, work, 1, a( k+1, k ),
383  $ 1 )
384  a( k, k-1 ) = a( k, k-1 ) -
385  $ sdot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
386  $ 1 )
387  CALL scopy( n-k, a( k+1, k-1 ), 1, work, 1 )
388  CALL ssymv( 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  $ sdot( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
405  CALL sswap( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
419  CALL sswap( 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 sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
434  CALL sswap( 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 SSYTRI_ROOK
449 *
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:53
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
ssymv
subroutine ssymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SSYMV
Definition: ssymv.f:154
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: