LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ssytri()

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.

Definition at line 113 of file ssytri.f.

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