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

SSYTRI

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

Purpose: 
 SSYTRI computes the inverse of a real symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 SSYTRF.
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.

          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.
[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
November 2011

Definition at line 116 of file ssytri.f.

116 *
117 * -- LAPACK computational routine (version 3.4.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2011
121 *
122 * .. Scalar Arguments ..
123  CHARACTER uplo
124  INTEGER info, lda, n
125 * ..
126 * .. Array Arguments ..
127  INTEGER ipiv( * )
128  REAL a( lda, * ), work( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  REAL one, zero
135  parameter ( one = 1.0e+0, zero = 0.0e+0 )
136 * ..
137 * .. Local Scalars ..
138  LOGICAL upper
139  INTEGER k, kp, kstep
140  REAL ak, akkp1, akp1, d, t, temp
141 * ..
142 * .. External Functions ..
143  LOGICAL lsame
144  REAL sdot
145  EXTERNAL lsame, sdot
146 * ..
147 * .. External Subroutines ..
148  EXTERNAL scopy, sswap, ssymv, xerbla
149 * ..
150 * .. Intrinsic Functions ..
151  INTRINSIC abs, max
152 * ..
153 * .. Executable Statements ..
154 *
155 * Test the input parameters.
156 *
157  info = 0
158  upper = lsame( uplo, 'U' )
159  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
160  info = -1
161  ELSE IF( n.LT.0 ) THEN
162  info = -2
163  ELSE IF( lda.LT.max( 1, n ) ) THEN
164  info = -4
165  END IF
166  IF( info.NE.0 ) THEN
167  CALL xerbla( 'SSYTRI', -info )
168  RETURN
169  END IF
170 *
171 * Quick return if possible
172 *
173  IF( n.EQ.0 )
174  $ RETURN
175 *
176 * Check that the diagonal matrix D is nonsingular.
177 *
178  IF( upper ) THEN
179 *
180 * Upper triangular storage: examine D from bottom to top
181 *
182  DO 10 info = n, 1, -1
183  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
184  $ RETURN
185  10 CONTINUE
186  ELSE
187 *
188 * Lower triangular storage: examine D from top to bottom.
189 *
190  DO 20 info = 1, n
191  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
192  $ RETURN
193  20 CONTINUE
194  END IF
195  info = 0
196 *
197  IF( upper ) THEN
198 *
199 * Compute inv(A) from the factorization A = U*D*U**T.
200 *
201 * K is the main loop index, increasing from 1 to N in steps of
202 * 1 or 2, depending on the size of the diagonal blocks.
203 *
204  k = 1
205  30 CONTINUE
206 *
207 * If K > N, exit from loop.
208 *
209  IF( k.GT.n )
210  $ GO TO 40
211 *
212  IF( ipiv( k ).GT.0 ) THEN
213 *
214 * 1 x 1 diagonal block
215 *
216 * Invert the diagonal block.
217 *
218  a( k, k ) = one / a( k, k )
219 *
220 * Compute column K of the inverse.
221 *
222  IF( k.GT.1 ) THEN
223  CALL scopy( k-1, a( 1, k ), 1, work, 1 )
224  CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
225  $ a( 1, k ), 1 )
226  a( k, k ) = a( k, k ) - sdot( k-1, work, 1, a( 1, k ),
227  $ 1 )
228  END IF
229  kstep = 1
230  ELSE
231 *
232 * 2 x 2 diagonal block
233 *
234 * Invert the diagonal block.
235 *
236  t = abs( a( k, k+1 ) )
237  ak = a( k, k ) / t
238  akp1 = a( k+1, k+1 ) / t
239  akkp1 = a( k, k+1 ) / t
240  d = t*( ak*akp1-one )
241  a( k, k ) = akp1 / d
242  a( k+1, k+1 ) = ak / d
243  a( k, k+1 ) = -akkp1 / d
244 *
245 * Compute columns K and K+1 of the inverse.
246 *
247  IF( k.GT.1 ) THEN
248  CALL scopy( k-1, a( 1, k ), 1, work, 1 )
249  CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
250  $ a( 1, k ), 1 )
251  a( k, k ) = a( k, k ) - sdot( k-1, work, 1, a( 1, k ),
252  $ 1 )
253  a( k, k+1 ) = a( k, k+1 ) -
254  $ sdot( k-1, a( 1, k ), 1, a( 1, k+1 ), 1 )
255  CALL scopy( k-1, a( 1, k+1 ), 1, work, 1 )
256  CALL ssymv( uplo, k-1, -one, a, lda, work, 1, zero,
257  $ a( 1, k+1 ), 1 )
258  a( k+1, k+1 ) = a( k+1, k+1 ) -
259  $ sdot( k-1, work, 1, a( 1, k+1 ), 1 )
260  END IF
261  kstep = 2
262  END IF
263 *
264  kp = abs( ipiv( k ) )
265  IF( kp.NE.k ) THEN
266 *
267 * Interchange rows and columns K and KP in the leading
268 * submatrix A(1:k+1,1:k+1)
269 *
270  CALL sswap( kp-1, a( 1, k ), 1, a( 1, kp ), 1 )
271  CALL sswap( k-kp-1, a( kp+1, k ), 1, a( kp, kp+1 ), lda )
272  temp = a( k, k )
273  a( k, k ) = a( kp, kp )
274  a( kp, kp ) = temp
275  IF( kstep.EQ.2 ) THEN
276  temp = a( k, k+1 )
277  a( k, k+1 ) = a( kp, k+1 )
278  a( kp, k+1 ) = temp
279  END IF
280  END IF
281 *
282  k = k + kstep
283  GO TO 30
284  40 CONTINUE
285 *
286  ELSE
287 *
288 * Compute inv(A) from the factorization A = L*D*L**T.
289 *
290 * K is the main loop index, increasing from 1 to N in steps of
291 * 1 or 2, depending on the size of the diagonal blocks.
292 *
293  k = n
294  50 CONTINUE
295 *
296 * If K < 1, exit from loop.
297 *
298  IF( k.LT.1 )
299  $ GO TO 60
300 *
301  IF( ipiv( k ).GT.0 ) THEN
302 *
303 * 1 x 1 diagonal block
304 *
305 * Invert the diagonal block.
306 *
307  a( k, k ) = one / a( k, k )
308 *
309 * Compute column K of the inverse.
310 *
311  IF( k.LT.n ) THEN
312  CALL scopy( n-k, a( k+1, k ), 1, work, 1 )
313  CALL ssymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
314  $ zero, a( k+1, k ), 1 )
315  a( k, k ) = a( k, k ) - sdot( n-k, work, 1, a( k+1, k ),
316  $ 1 )
317  END IF
318  kstep = 1
319  ELSE
320 *
321 * 2 x 2 diagonal block
322 *
323 * Invert the diagonal block.
324 *
325  t = abs( a( k, k-1 ) )
326  ak = a( k-1, k-1 ) / t
327  akp1 = a( k, k ) / t
328  akkp1 = a( k, k-1 ) / t
329  d = t*( ak*akp1-one )
330  a( k-1, k-1 ) = akp1 / d
331  a( k, k ) = ak / d
332  a( k, k-1 ) = -akkp1 / d
333 *
334 * Compute columns K-1 and K of the inverse.
335 *
336  IF( k.LT.n ) THEN
337  CALL scopy( n-k, a( k+1, k ), 1, work, 1 )
338  CALL ssymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
339  $ zero, a( k+1, k ), 1 )
340  a( k, k ) = a( k, k ) - sdot( n-k, work, 1, a( k+1, k ),
341  $ 1 )
342  a( k, k-1 ) = a( k, k-1 ) -
343  $ sdot( n-k, a( k+1, k ), 1, a( k+1, k-1 ),
344  $ 1 )
345  CALL scopy( n-k, a( k+1, k-1 ), 1, work, 1 )
346  CALL ssymv( uplo, n-k, -one, a( k+1, k+1 ), lda, work, 1,
347  $ zero, a( k+1, k-1 ), 1 )
348  a( k-1, k-1 ) = a( k-1, k-1 ) -
349  $ sdot( n-k, work, 1, a( k+1, k-1 ), 1 )
350  END IF
351  kstep = 2
352  END IF
353 *
354  kp = abs( ipiv( k ) )
355  IF( kp.NE.k ) THEN
356 *
357 * Interchange rows and columns K and KP in the trailing
358 * submatrix A(k-1:n,k-1:n)
359 *
360  IF( kp.LT.n )
361  $ CALL sswap( n-kp, a( kp+1, k ), 1, a( kp+1, kp ), 1 )
362  CALL sswap( kp-k-1, a( k+1, k ), 1, a( kp, k+1 ), lda )
363  temp = a( k, k )
364  a( k, k ) = a( kp, kp )
365  a( kp, kp ) = temp
366  IF( kstep.EQ.2 ) THEN
367  temp = a( k, k-1 )
368  a( k, k-1 ) = a( kp, k-1 )
369  a( kp, k-1 ) = temp
370  END IF
371  END IF
372 *
373  k = k - kstep
374  GO TO 50
375  60 CONTINUE
376  END IF
377 *
378  RETURN
379 *
380 * End of SSYTRI
381 *
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: