LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssytri.f
Go to the documentation of this file.
1 *> \brief \b SSYTRI
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSYTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * REAL A( LDA, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> SSYTRI computes the inverse of a real symmetric indefinite matrix
39 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40 *> SSYTRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**T;
52 *> = 'L': Lower triangular, form is A = L*D*L**T.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] A
62 *> \verbatim
63 *> A is REAL array, dimension (LDA,N)
64 *> On entry, the block diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by SSYTRF.
66 *>
67 *> On exit, if INFO = 0, the (symmetric) inverse of the original
68 *> matrix. If UPLO = 'U', the upper triangular part of the
69 *> inverse is formed and the part of A below the diagonal is not
70 *> referenced; if UPLO = 'L' the lower triangular part of the
71 *> inverse is formed and the part of A above the diagonal is
72 *> not referenced.
73 *> \endverbatim
74 *>
75 *> \param[in] LDA
76 *> \verbatim
77 *> LDA is INTEGER
78 *> The leading dimension of the array A. LDA >= max(1,N).
79 *> \endverbatim
80 *>
81 *> \param[in] IPIV
82 *> \verbatim
83 *> IPIV is INTEGER array, dimension (N)
84 *> Details of the interchanges and the block structure of D
85 *> as determined by SSYTRF.
86 *> \endverbatim
87 *>
88 *> \param[out] WORK
89 *> \verbatim
90 *> WORK is REAL array, dimension (N)
91 *> \endverbatim
92 *>
93 *> \param[out] INFO
94 *> \verbatim
95 *> INFO is INTEGER
96 *> = 0: successful exit
97 *> < 0: if INFO = -i, the i-th argument had an illegal value
98 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
99 *> inverse could not be computed.
100 *> \endverbatim
101 *
102 * Authors:
103 * ========
104 *
105 *> \author Univ. of Tennessee
106 *> \author Univ. of California Berkeley
107 *> \author Univ. of Colorado Denver
108 *> \author NAG Ltd.
109 *
110 *> \ingroup realSYcomputational
111 *
112 * =====================================================================
113  SUBROUTINE ssytri( UPLO, N, A, LDA, IPIV, WORK, INFO )
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 *
379  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssytri(UPLO, N, A, LDA, IPIV, WORK, INFO)
SSYTRI
Definition: ssytri.f:114
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine ssymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SSYMV
Definition: ssymv.f:152