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
380