1*> \brief \b ZSYTRI
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZSYTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZSYTRI( 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*       COMPLEX*16         A( LDA, * ), WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> ZSYTRI computes the inverse of a complex symmetric indefinite matrix
39*> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40*> ZSYTRF.
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 COMPLEX*16 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 ZSYTRF.
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 ZSYTRF.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*>          WORK is COMPLEX*16 array, dimension (2*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 complex16SYcomputational
111*
112*  =====================================================================
113      SUBROUTINE ZSYTRI( 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      COMPLEX*16         A( LDA, * ), WORK( * )
126*     ..
127*
128*  =====================================================================
129*
130*     .. Parameters ..
131      COMPLEX*16         ONE, ZERO
132      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
133     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
134*     ..
135*     .. Local Scalars ..
136      LOGICAL            UPPER
137      INTEGER            K, KP, KSTEP
138      COMPLEX*16         AK, AKKP1, AKP1, D, T, TEMP
139*     ..
140*     .. External Functions ..
141      LOGICAL            LSAME
142      COMPLEX*16         ZDOTU
143      EXTERNAL           LSAME, ZDOTU
144*     ..
145*     .. External Subroutines ..
146      EXTERNAL           XERBLA, ZCOPY, ZSWAP, ZSYMV
147*     ..
148*     .. Intrinsic Functions ..
149      INTRINSIC          ABS, MAX
150*     ..
151*     .. Executable Statements ..
152*
153*     Test the input parameters.
154*
155      INFO = 0
156      UPPER = LSAME( UPLO, 'U' )
157      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
158         INFO = -1
159      ELSE IF( N.LT.0 ) THEN
160         INFO = -2
161      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162         INFO = -4
163      END IF
164      IF( INFO.NE.0 ) THEN
165         CALL XERBLA( 'ZSYTRI', -INFO )
166         RETURN
167      END IF
168*
169*     Quick return if possible
170*
171      IF( N.EQ.0 )
172     $   RETURN
173*
174*     Check that the diagonal matrix D is nonsingular.
175*
176      IF( UPPER ) THEN
177*
178*        Upper triangular storage: examine D from bottom to top
179*
180         DO 10 INFO = N, 1, -1
181            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
182     $         RETURN
183   10    CONTINUE
184      ELSE
185*
186*        Lower triangular storage: examine D from top to bottom.
187*
188         DO 20 INFO = 1, N
189            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
190     $         RETURN
191   20    CONTINUE
192      END IF
193      INFO = 0
194*
195      IF( UPPER ) THEN
196*
197*        Compute inv(A) from the factorization A = U*D*U**T.
198*
199*        K is the main loop index, increasing from 1 to N in steps of
200*        1 or 2, depending on the size of the diagonal blocks.
201*
202         K = 1
203   30    CONTINUE
204*
205*        If K > N, exit from loop.
206*
207         IF( K.GT.N )
208     $      GO TO 40
209*
210         IF( IPIV( K ).GT.0 ) THEN
211*
212*           1 x 1 diagonal block
213*
214*           Invert the diagonal block.
215*
216            A( K, K ) = ONE / A( K, K )
217*
218*           Compute column K of the inverse.
219*
220            IF( K.GT.1 ) THEN
221               CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
222               CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
223     $                     A( 1, K ), 1 )
224               A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
225     $                     1 )
226            END IF
227            KSTEP = 1
228         ELSE
229*
230*           2 x 2 diagonal block
231*
232*           Invert the diagonal block.
233*
234            T = A( K, K+1 )
235            AK = A( K, K ) / T
236            AKP1 = A( K+1, K+1 ) / T
237            AKKP1 = A( K, K+1 ) / T
238            D = T*( AK*AKP1-ONE )
239            A( K, K ) = AKP1 / D
240            A( K+1, K+1 ) = AK / D
241            A( K, K+1 ) = -AKKP1 / D
242*
243*           Compute columns K and K+1 of the inverse.
244*
245            IF( K.GT.1 ) THEN
246               CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
247               CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
248     $                     A( 1, K ), 1 )
249               A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
250     $                     1 )
251               A( K, K+1 ) = A( K, K+1 ) -
252     $                       ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
253               CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
254               CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
255     $                     A( 1, K+1 ), 1 )
256               A( K+1, K+1 ) = A( K+1, K+1 ) -
257     $                         ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
258            END IF
259            KSTEP = 2
260         END IF
261*
262         KP = ABS( IPIV( K ) )
263         IF( KP.NE.K ) THEN
264*
265*           Interchange rows and columns K and KP in the leading
266*           submatrix A(1:k+1,1:k+1)
267*
268            CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
269            CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
270            TEMP = A( K, K )
271            A( K, K ) = A( KP, KP )
272            A( KP, KP ) = TEMP
273            IF( KSTEP.EQ.2 ) THEN
274               TEMP = A( K, K+1 )
275               A( K, K+1 ) = A( KP, K+1 )
276               A( KP, K+1 ) = TEMP
277            END IF
278         END IF
279*
280         K = K + KSTEP
281         GO TO 30
282   40    CONTINUE
283*
284      ELSE
285*
286*        Compute inv(A) from the factorization A = L*D*L**T.
287*
288*        K is the main loop index, increasing from 1 to N in steps of
289*        1 or 2, depending on the size of the diagonal blocks.
290*
291         K = N
292   50    CONTINUE
293*
294*        If K < 1, exit from loop.
295*
296         IF( K.LT.1 )
297     $      GO TO 60
298*
299         IF( IPIV( K ).GT.0 ) THEN
300*
301*           1 x 1 diagonal block
302*
303*           Invert the diagonal block.
304*
305            A( K, K ) = ONE / A( K, K )
306*
307*           Compute column K of the inverse.
308*
309            IF( K.LT.N ) THEN
310               CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
311               CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
312     $                     ZERO, A( K+1, K ), 1 )
313               A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
314     $                     1 )
315            END IF
316            KSTEP = 1
317         ELSE
318*
319*           2 x 2 diagonal block
320*
321*           Invert the diagonal block.
322*
323            T = A( K, K-1 )
324            AK = A( K-1, K-1 ) / T
325            AKP1 = A( K, K ) / T
326            AKKP1 = A( K, K-1 ) / T
327            D = T*( AK*AKP1-ONE )
328            A( K-1, K-1 ) = AKP1 / D
329            A( K, K ) = AK / D
330            A( K, K-1 ) = -AKKP1 / D
331*
332*           Compute columns K-1 and K of the inverse.
333*
334            IF( K.LT.N ) THEN
335               CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
336               CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
337     $                     ZERO, A( K+1, K ), 1 )
338               A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
339     $                     1 )
340               A( K, K-1 ) = A( K, K-1 ) -
341     $                       ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
342     $                       1 )
343               CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
344               CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
345     $                     ZERO, A( K+1, K-1 ), 1 )
346               A( K-1, K-1 ) = A( K-1, K-1 ) -
347     $                         ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
348            END IF
349            KSTEP = 2
350         END IF
351*
352         KP = ABS( IPIV( K ) )
353         IF( KP.NE.K ) THEN
354*
355*           Interchange rows and columns K and KP in the trailing
356*           submatrix A(k-1:n,k-1:n)
357*
358            IF( KP.LT.N )
359     $         CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
360            CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
361            TEMP = A( K, K )
362            A( K, K ) = A( KP, KP )
363            A( KP, KP ) = TEMP
364            IF( KSTEP.EQ.2 ) THEN
365               TEMP = A( K, K-1 )
366               A( K, K-1 ) = A( KP, K-1 )
367               A( KP, K-1 ) = TEMP
368            END IF
369         END IF
370*
371         K = K - KSTEP
372         GO TO 50
373   60    CONTINUE
374      END IF
375*
376      RETURN
377*
378*     End of ZSYTRI
379*
380      END
381