1*> \brief \b ZSPTRI
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZSPTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsptri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsptri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsptri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, N
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            IPIV( * )
29*       COMPLEX*16         AP( * ), WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> ZSPTRI computes the inverse of a complex symmetric indefinite matrix
39*> A in packed storage using the factorization A = U*D*U**T or
40*> A = L*D*L**T computed by ZSPTRF.
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] AP
62*> \verbatim
63*>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
64*>          On entry, the block diagonal matrix D and the multipliers
65*>          used to obtain the factor U or L as computed by ZSPTRF,
66*>          stored as a packed triangular matrix.
67*>
68*>          On exit, if INFO = 0, the (symmetric) inverse of the original
69*>          matrix, stored as a packed triangular matrix. The j-th column
70*>          of inv(A) is stored in the array AP as follows:
71*>          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
72*>          if UPLO = 'L',
73*>             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
74*> \endverbatim
75*>
76*> \param[in] IPIV
77*> \verbatim
78*>          IPIV is INTEGER array, dimension (N)
79*>          Details of the interchanges and the block structure of D
80*>          as determined by ZSPTRF.
81*> \endverbatim
82*>
83*> \param[out] WORK
84*> \verbatim
85*>          WORK is COMPLEX*16 array, dimension (N)
86*> \endverbatim
87*>
88*> \param[out] INFO
89*> \verbatim
90*>          INFO is INTEGER
91*>          = 0: successful exit
92*>          < 0: if INFO = -i, the i-th argument had an illegal value
93*>          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
94*>               inverse could not be computed.
95*> \endverbatim
96*
97*  Authors:
98*  ========
99*
100*> \author Univ. of Tennessee
101*> \author Univ. of California Berkeley
102*> \author Univ. of Colorado Denver
103*> \author NAG Ltd.
104*
105*> \date November 2011
106*
107*> \ingroup complex16OTHERcomputational
108*
109*  =====================================================================
110      SUBROUTINE ZSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
111*
112*  -- LAPACK computational routine (version 3.4.0) --
113*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
114*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115*     November 2011
116*
117*     .. Scalar Arguments ..
118      CHARACTER          UPLO
119      INTEGER            INFO, N
120*     ..
121*     .. Array Arguments ..
122      INTEGER            IPIV( * )
123      COMPLEX*16         AP( * ), WORK( * )
124*     ..
125*
126*  =====================================================================
127*
128*     .. Parameters ..
129      COMPLEX*16         ONE, ZERO
130      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
131     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
132*     ..
133*     .. Local Scalars ..
134      LOGICAL            UPPER
135      INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
136      COMPLEX*16         AK, AKKP1, AKP1, D, T, TEMP
137*     ..
138*     .. External Functions ..
139      LOGICAL            LSAME
140      COMPLEX*16         ZDOTU
141      EXTERNAL           LSAME, ZDOTU
142*     ..
143*     .. External Subroutines ..
144      EXTERNAL           XERBLA, ZCOPY, ZSPMV, ZSWAP
145*     ..
146*     .. Intrinsic Functions ..
147      INTRINSIC          ABS
148*     ..
149*     .. Executable Statements ..
150*
151*     Test the input parameters.
152*
153      INFO = 0
154      UPPER = LSAME( UPLO, 'U' )
155      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
156         INFO = -1
157      ELSE IF( N.LT.0 ) THEN
158         INFO = -2
159      END IF
160      IF( INFO.NE.0 ) THEN
161         CALL XERBLA( 'ZSPTRI', -INFO )
162         RETURN
163      END IF
164*
165*     Quick return if possible
166*
167      IF( N.EQ.0 )
168     $   RETURN
169*
170*     Check that the diagonal matrix D is nonsingular.
171*
172      IF( UPPER ) THEN
173*
174*        Upper triangular storage: examine D from bottom to top
175*
176         KP = N*( N+1 ) / 2
177         DO 10 INFO = N, 1, -1
178            IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
179     $         RETURN
180            KP = KP - INFO
181   10    CONTINUE
182      ELSE
183*
184*        Lower triangular storage: examine D from top to bottom.
185*
186         KP = 1
187         DO 20 INFO = 1, N
188            IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
189     $         RETURN
190            KP = KP + N - INFO + 1
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         KC = 1
204   30    CONTINUE
205*
206*        If K > N, exit from loop.
207*
208         IF( K.GT.N )
209     $      GO TO 50
210*
211         KCNEXT = KC + K
212         IF( IPIV( K ).GT.0 ) THEN
213*
214*           1 x 1 diagonal block
215*
216*           Invert the diagonal block.
217*
218            AP( KC+K-1 ) = ONE / AP( KC+K-1 )
219*
220*           Compute column K of the inverse.
221*
222            IF( K.GT.1 ) THEN
223               CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
224               CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
225     $                     1 )
226               AP( KC+K-1 ) = AP( KC+K-1 ) -
227     $                        ZDOTU( K-1, WORK, 1, AP( KC ), 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 = AP( KCNEXT+K-1 )
237            AK = AP( KC+K-1 ) / T
238            AKP1 = AP( KCNEXT+K ) / T
239            AKKP1 = AP( KCNEXT+K-1 ) / T
240            D = T*( AK*AKP1-ONE )
241            AP( KC+K-1 ) = AKP1 / D
242            AP( KCNEXT+K ) = AK / D
243            AP( KCNEXT+K-1 ) = -AKKP1 / D
244*
245*           Compute columns K and K+1 of the inverse.
246*
247            IF( K.GT.1 ) THEN
248               CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
249               CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
250     $                     1 )
251               AP( KC+K-1 ) = AP( KC+K-1 ) -
252     $                        ZDOTU( K-1, WORK, 1, AP( KC ), 1 )
253               AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
254     $                            ZDOTU( K-1, AP( KC ), 1, AP( KCNEXT ),
255     $                            1 )
256               CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
257               CALL ZSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
258     $                     AP( KCNEXT ), 1 )
259               AP( KCNEXT+K ) = AP( KCNEXT+K ) -
260     $                          ZDOTU( K-1, WORK, 1, AP( KCNEXT ), 1 )
261            END IF
262            KSTEP = 2
263            KCNEXT = KCNEXT + K + 1
264         END IF
265*
266         KP = ABS( IPIV( K ) )
267         IF( KP.NE.K ) THEN
268*
269*           Interchange rows and columns K and KP in the leading
270*           submatrix A(1:k+1,1:k+1)
271*
272            KPC = ( KP-1 )*KP / 2 + 1
273            CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
274            KX = KPC + KP - 1
275            DO 40 J = KP + 1, K - 1
276               KX = KX + J - 1
277               TEMP = AP( KC+J-1 )
278               AP( KC+J-1 ) = AP( KX )
279               AP( KX ) = TEMP
280   40       CONTINUE
281            TEMP = AP( KC+K-1 )
282            AP( KC+K-1 ) = AP( KPC+KP-1 )
283            AP( KPC+KP-1 ) = TEMP
284            IF( KSTEP.EQ.2 ) THEN
285               TEMP = AP( KC+K+K-1 )
286               AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
287               AP( KC+K+KP-1 ) = TEMP
288            END IF
289         END IF
290*
291         K = K + KSTEP
292         KC = KCNEXT
293         GO TO 30
294   50    CONTINUE
295*
296      ELSE
297*
298*        Compute inv(A) from the factorization A = L*D*L**T.
299*
300*        K is the main loop index, increasing from 1 to N in steps of
301*        1 or 2, depending on the size of the diagonal blocks.
302*
303         NPP = N*( N+1 ) / 2
304         K = N
305         KC = NPP
306   60    CONTINUE
307*
308*        If K < 1, exit from loop.
309*
310         IF( K.LT.1 )
311     $      GO TO 80
312*
313         KCNEXT = KC - ( N-K+2 )
314         IF( IPIV( K ).GT.0 ) THEN
315*
316*           1 x 1 diagonal block
317*
318*           Invert the diagonal block.
319*
320            AP( KC ) = ONE / AP( KC )
321*
322*           Compute column K of the inverse.
323*
324            IF( K.LT.N ) THEN
325               CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
326               CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
327     $                     ZERO, AP( KC+1 ), 1 )
328               AP( KC ) = AP( KC ) - ZDOTU( N-K, WORK, 1, AP( KC+1 ),
329     $                    1 )
330            END IF
331            KSTEP = 1
332         ELSE
333*
334*           2 x 2 diagonal block
335*
336*           Invert the diagonal block.
337*
338            T = AP( KCNEXT+1 )
339            AK = AP( KCNEXT ) / T
340            AKP1 = AP( KC ) / T
341            AKKP1 = AP( KCNEXT+1 ) / T
342            D = T*( AK*AKP1-ONE )
343            AP( KCNEXT ) = AKP1 / D
344            AP( KC ) = AK / D
345            AP( KCNEXT+1 ) = -AKKP1 / D
346*
347*           Compute columns K-1 and K of the inverse.
348*
349            IF( K.LT.N ) THEN
350               CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
351               CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
352     $                     ZERO, AP( KC+1 ), 1 )
353               AP( KC ) = AP( KC ) - ZDOTU( N-K, WORK, 1, AP( KC+1 ),
354     $                    1 )
355               AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
356     $                          ZDOTU( N-K, AP( KC+1 ), 1,
357     $                          AP( KCNEXT+2 ), 1 )
358               CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
359               CALL ZSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
360     $                     ZERO, AP( KCNEXT+2 ), 1 )
361               AP( KCNEXT ) = AP( KCNEXT ) -
362     $                        ZDOTU( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
363            END IF
364            KSTEP = 2
365            KCNEXT = KCNEXT - ( N-K+3 )
366         END IF
367*
368         KP = ABS( IPIV( K ) )
369         IF( KP.NE.K ) THEN
370*
371*           Interchange rows and columns K and KP in the trailing
372*           submatrix A(k-1:n,k-1:n)
373*
374            KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
375            IF( KP.LT.N )
376     $         CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
377            KX = KC + KP - K
378            DO 70 J = K + 1, KP - 1
379               KX = KX + N - J + 1
380               TEMP = AP( KC+J-K )
381               AP( KC+J-K ) = AP( KX )
382               AP( KX ) = TEMP
383   70       CONTINUE
384            TEMP = AP( KC )
385            AP( KC ) = AP( KPC )
386            AP( KPC ) = TEMP
387            IF( KSTEP.EQ.2 ) THEN
388               TEMP = AP( KC-N+K-1 )
389               AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
390               AP( KC-N+KP-1 ) = TEMP
391            END IF
392         END IF
393*
394         K = K - KSTEP
395         KC = KCNEXT
396         GO TO 60
397   80    CONTINUE
398      END IF
399*
400      RETURN
401*
402*     End of ZSPTRI
403*
404      END
405