1*> \brief \b CSYTRI_ROOK
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CSYTRI_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CSYTRI_ROOK( 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            A( LDA, * ), WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> CSYTRI_ROOK computes the inverse of a complex symmetric
39*> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
40*> computed by CSYTRF_ROOK.
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 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 CSYTRF_ROOK.
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 CSYTRF_ROOK.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*>          WORK is COMPLEX 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*> \date November 2015
111*
112*> \ingroup complexSYcomputational
113*
114*> \par Contributors:
115*  ==================
116*>
117*> \verbatim
118*>
119*>   November 2015, Igor Kozachenko,
120*>                  Computer Science Division,
121*>                  University of California, Berkeley
122*>
123*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
124*>                  School of Mathematics,
125*>                  University of Manchester
126*>
127*> \endverbatim
128*
129*  =====================================================================
130      SUBROUTINE CSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
131*
132*  -- LAPACK computational routine (version 3.6.0) --
133*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
134*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*     November 2015
136*
137*     .. Scalar Arguments ..
138      CHARACTER          UPLO
139      INTEGER            INFO, LDA, N
140*     ..
141*     .. Array Arguments ..
142      INTEGER            IPIV( * )
143      COMPLEX            A( LDA, * ), WORK( * )
144*     ..
145*
146*  =====================================================================
147*
148*     .. Parameters ..
149      COMPLEX            CONE, CZERO
150      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ),
151     $                   CZERO = ( 0.0E+0, 0.0E+0 ) )
152*     ..
153*     .. Local Scalars ..
154      LOGICAL            UPPER
155      INTEGER            K, KP, KSTEP
156      COMPLEX            AK, AKKP1, AKP1, D, T, TEMP
157*     ..
158*     .. External Functions ..
159      LOGICAL            LSAME
160      COMPLEX            CDOTU
161      EXTERNAL           LSAME, CDOTU
162*     ..
163*     .. External Subroutines ..
164      EXTERNAL           CCOPY, CSWAP, CSYMV, XERBLA
165*     ..
166*     .. Intrinsic Functions ..
167      INTRINSIC          MAX
168*     ..
169*     .. Executable Statements ..
170*
171*     Test the input parameters.
172*
173      INFO = 0
174      UPPER = LSAME( UPLO, 'U' )
175      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
176         INFO = -1
177      ELSE IF( N.LT.0 ) THEN
178         INFO = -2
179      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
180         INFO = -4
181      END IF
182      IF( INFO.NE.0 ) THEN
183         CALL XERBLA( 'CSYTRI_ROOK', -INFO )
184         RETURN
185      END IF
186*
187*     Quick return if possible
188*
189      IF( N.EQ.0 )
190     $   RETURN
191*
192*     Check that the diagonal matrix D is nonsingular.
193*
194      IF( UPPER ) THEN
195*
196*        Upper triangular storage: examine D from bottom to top
197*
198         DO 10 INFO = N, 1, -1
199            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
200     $         RETURN
201   10    CONTINUE
202      ELSE
203*
204*        Lower triangular storage: examine D from top to bottom.
205*
206         DO 20 INFO = 1, N
207            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
208     $         RETURN
209   20    CONTINUE
210      END IF
211      INFO = 0
212*
213      IF( UPPER ) THEN
214*
215*        Compute inv(A) from the factorization A = U*D*U**T.
216*
217*        K is the main loop index, increasing from 1 to N in steps of
218*        1 or 2, depending on the size of the diagonal blocks.
219*
220         K = 1
221   30    CONTINUE
222*
223*        If K > N, exit from loop.
224*
225         IF( K.GT.N )
226     $      GO TO 40
227*
228         IF( IPIV( K ).GT.0 ) THEN
229*
230*           1 x 1 diagonal block
231*
232*           Invert the diagonal block.
233*
234            A( K, K ) = CONE / A( K, K )
235*
236*           Compute column K of the inverse.
237*
238            IF( K.GT.1 ) THEN
239               CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
240               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
241     $                     A( 1, K ), 1 )
242               A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ),
243     $                     1 )
244            END IF
245            KSTEP = 1
246         ELSE
247*
248*           2 x 2 diagonal block
249*
250*           Invert the diagonal block.
251*
252            T = A( K, K+1 )
253            AK = A( K, K ) / T
254            AKP1 = A( K+1, K+1 ) / T
255            AKKP1 = A( K, K+1 ) / T
256            D = T*( AK*AKP1-CONE )
257            A( K, K ) = AKP1 / D
258            A( K+1, K+1 ) = AK / D
259            A( K, K+1 ) = -AKKP1 / D
260*
261*           Compute columns K and K+1 of the inverse.
262*
263            IF( K.GT.1 ) THEN
264               CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
265               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
266     $                     A( 1, K ), 1 )
267               A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ),
268     $                     1 )
269               A( K, K+1 ) = A( K, K+1 ) -
270     $                       CDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
271               CALL CCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
272               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
273     $                     A( 1, K+1 ), 1 )
274               A( K+1, K+1 ) = A( K+1, K+1 ) -
275     $                         CDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
276            END IF
277            KSTEP = 2
278         END IF
279*
280         IF( KSTEP.EQ.1 ) THEN
281*
282*           Interchange rows and columns K and IPIV(K) in the leading
283*           submatrix A(1:k+1,1:k+1)
284*
285            KP = IPIV( K )
286            IF( KP.NE.K ) THEN
287               IF( KP.GT.1 )
288     $             CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
289               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
290               TEMP = A( K, K )
291               A( K, K ) = A( KP, KP )
292               A( KP, KP ) = TEMP
293            END IF
294         ELSE
295*
296*           Interchange rows and columns K and K+1 with -IPIV(K) and
297*           -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
298*
299            KP = -IPIV( K )
300            IF( KP.NE.K ) THEN
301               IF( KP.GT.1 )
302     $            CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
303               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
304*
305               TEMP = A( K, K )
306               A( K, K ) = A( KP, KP )
307               A( KP, KP ) = TEMP
308               TEMP = A( K, K+1 )
309               A( K, K+1 ) = A( KP, K+1 )
310               A( KP, K+1 ) = TEMP
311            END IF
312*
313            K = K + 1
314            KP = -IPIV( K )
315            IF( KP.NE.K ) THEN
316               IF( KP.GT.1 )
317     $            CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
318               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
319               TEMP = A( K, K )
320               A( K, K ) = A( KP, KP )
321               A( KP, KP ) = TEMP
322            END IF
323         END IF
324*
325         K = K + 1
326         GO TO 30
327   40    CONTINUE
328*
329      ELSE
330*
331*        Compute inv(A) from the factorization A = L*D*L**T.
332*
333*        K is the main loop index, increasing from 1 to N in steps of
334*        1 or 2, depending on the size of the diagonal blocks.
335*
336         K = N
337   50    CONTINUE
338*
339*        If K < 1, exit from loop.
340*
341         IF( K.LT.1 )
342     $      GO TO 60
343*
344         IF( IPIV( K ).GT.0 ) THEN
345*
346*           1 x 1 diagonal block
347*
348*           Invert the diagonal block.
349*
350            A( K, K ) = CONE / A( K, K )
351*
352*           Compute column K of the inverse.
353*
354            IF( K.LT.N ) THEN
355               CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
356               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
357     $                     CZERO, A( K+1, K ), 1 )
358               A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ),
359     $                     1 )
360            END IF
361            KSTEP = 1
362         ELSE
363*
364*           2 x 2 diagonal block
365*
366*           Invert the diagonal block.
367*
368            T = A( K, K-1 )
369            AK = A( K-1, K-1 ) / T
370            AKP1 = A( K, K ) / T
371            AKKP1 = A( K, K-1 ) / T
372            D = T*( AK*AKP1-CONE )
373            A( K-1, K-1 ) = AKP1 / D
374            A( K, K ) = AK / D
375            A( K, K-1 ) = -AKKP1 / D
376*
377*           Compute columns K-1 and K of the inverse.
378*
379            IF( K.LT.N ) THEN
380               CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
381               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
382     $                     CZERO, A( K+1, K ), 1 )
383               A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ),
384     $                     1 )
385               A( K, K-1 ) = A( K, K-1 ) -
386     $                       CDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
387     $                       1 )
388               CALL CCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
389               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
390     $                     CZERO, A( K+1, K-1 ), 1 )
391               A( K-1, K-1 ) = A( K-1, K-1 ) -
392     $                         CDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
393            END IF
394            KSTEP = 2
395         END IF
396*
397         IF( KSTEP.EQ.1 ) THEN
398*
399*           Interchange rows and columns K and IPIV(K) in the trailing
400*           submatrix A(k-1:n,k-1:n)
401*
402            KP = IPIV( K )
403            IF( KP.NE.K ) THEN
404               IF( KP.LT.N )
405     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
406               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
407               TEMP = A( K, K )
408               A( K, K ) = A( KP, KP )
409               A( KP, KP ) = TEMP
410            END IF
411         ELSE
412*
413*           Interchange rows and columns K and K-1 with -IPIV(K) and
414*           -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
415*
416            KP = -IPIV( K )
417            IF( KP.NE.K ) THEN
418               IF( KP.LT.N )
419     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
420               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
421*
422               TEMP = A( K, K )
423               A( K, K ) = A( KP, KP )
424               A( KP, KP ) = TEMP
425               TEMP = A( K, K-1 )
426               A( K, K-1 ) = A( KP, K-1 )
427               A( KP, K-1 ) = TEMP
428            END IF
429*
430            K = K - 1
431            KP = -IPIV( K )
432            IF( KP.NE.K ) THEN
433               IF( KP.LT.N )
434     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
435               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
436               TEMP = A( K, K )
437               A( K, K ) = A( KP, KP )
438               A( KP, KP ) = TEMP
439            END IF
440         END IF
441*
442         K = K - 1
443         GO TO 50
444   60    CONTINUE
445      END IF
446*
447      RETURN
448*
449*     End of CSYTRI_ROOK
450*
451      END
452