1*> \brief \b CPFTRI
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CPFTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpftri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpftri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpftri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CPFTRI( TRANSR, UPLO, N, A, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          TRANSR, UPLO
25*       INTEGER            INFO, N
26*       .. Array Arguments ..
27*       COMPLEX            A( 0: * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> CPFTRI computes the inverse of a complex Hermitian positive definite
37*> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
38*> computed by CPFTRF.
39*> \endverbatim
40*
41*  Arguments:
42*  ==========
43*
44*> \param[in] TRANSR
45*> \verbatim
46*>          TRANSR is CHARACTER*1
47*>          = 'N':  The Normal TRANSR of RFP A is stored;
48*>          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
49*> \endverbatim
50*>
51*> \param[in] UPLO
52*> \verbatim
53*>          UPLO is CHARACTER*1
54*>          = 'U':  Upper triangle of A is stored;
55*>          = 'L':  Lower triangle of A is stored.
56*> \endverbatim
57*>
58*> \param[in] N
59*> \verbatim
60*>          N is INTEGER
61*>          The order of the matrix A.  N >= 0.
62*> \endverbatim
63*>
64*> \param[in,out] A
65*> \verbatim
66*>          A is COMPLEX array, dimension ( N*(N+1)/2 );
67*>          On entry, the Hermitian matrix A in RFP format. RFP format is
68*>          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
69*>          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
70*>          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
71*>          the Conjugate-transpose of RFP A as defined when
72*>          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
73*>          follows: If UPLO = 'U' the RFP A contains the nt elements of
74*>          upper packed A. If UPLO = 'L' the RFP A contains the elements
75*>          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
76*>          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
77*>          is odd. See the Note below for more details.
78*>
79*>          On exit, the Hermitian inverse of the original matrix, in the
80*>          same storage format.
81*> \endverbatim
82*>
83*> \param[out] INFO
84*> \verbatim
85*>          INFO is INTEGER
86*>          = 0:  successful exit
87*>          < 0:  if INFO = -i, the i-th argument had an illegal value
88*>          > 0:  if INFO = i, the (i,i) element of the factor U or L is
89*>                zero, and the inverse could not be computed.
90*> \endverbatim
91*
92*  Authors:
93*  ========
94*
95*> \author Univ. of Tennessee
96*> \author Univ. of California Berkeley
97*> \author Univ. of Colorado Denver
98*> \author NAG Ltd.
99*
100*> \ingroup complexOTHERcomputational
101*
102*> \par Further Details:
103*  =====================
104*>
105*> \verbatim
106*>
107*>  We first consider Standard Packed Format when N is even.
108*>  We give an example where N = 6.
109*>
110*>      AP is Upper             AP is Lower
111*>
112*>   00 01 02 03 04 05       00
113*>      11 12 13 14 15       10 11
114*>         22 23 24 25       20 21 22
115*>            33 34 35       30 31 32 33
116*>               44 45       40 41 42 43 44
117*>                  55       50 51 52 53 54 55
118*>
119*>
120*>  Let TRANSR = 'N'. RFP holds AP as follows:
121*>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
122*>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
123*>  conjugate-transpose of the first three columns of AP upper.
124*>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
125*>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
126*>  conjugate-transpose of the last three columns of AP lower.
127*>  To denote conjugate we place -- above the element. This covers the
128*>  case N even and TRANSR = 'N'.
129*>
130*>         RFP A                   RFP A
131*>
132*>                                -- -- --
133*>        03 04 05                33 43 53
134*>                                   -- --
135*>        13 14 15                00 44 54
136*>                                      --
137*>        23 24 25                10 11 55
138*>
139*>        33 34 35                20 21 22
140*>        --
141*>        00 44 45                30 31 32
142*>        -- --
143*>        01 11 55                40 41 42
144*>        -- -- --
145*>        02 12 22                50 51 52
146*>
147*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
148*>  transpose of RFP A above. One therefore gets:
149*>
150*>
151*>           RFP A                   RFP A
152*>
153*>     -- -- -- --                -- -- -- -- -- --
154*>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
155*>     -- -- -- -- --                -- -- -- -- --
156*>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
157*>     -- -- -- -- -- --                -- -- -- --
158*>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
159*>
160*>
161*>  We next  consider Standard Packed Format when N is odd.
162*>  We give an example where N = 5.
163*>
164*>     AP is Upper                 AP is Lower
165*>
166*>   00 01 02 03 04              00
167*>      11 12 13 14              10 11
168*>         22 23 24              20 21 22
169*>            33 34              30 31 32 33
170*>               44              40 41 42 43 44
171*>
172*>
173*>  Let TRANSR = 'N'. RFP holds AP as follows:
174*>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
175*>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
176*>  conjugate-transpose of the first two   columns of AP upper.
177*>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
178*>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
179*>  conjugate-transpose of the last two   columns of AP lower.
180*>  To denote conjugate we place -- above the element. This covers the
181*>  case N odd  and TRANSR = 'N'.
182*>
183*>         RFP A                   RFP A
184*>
185*>                                   -- --
186*>        02 03 04                00 33 43
187*>                                      --
188*>        12 13 14                10 11 44
189*>
190*>        22 23 24                20 21 22
191*>        --
192*>        00 33 34                30 31 32
193*>        -- --
194*>        01 11 44                40 41 42
195*>
196*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
197*>  transpose of RFP A above. One therefore gets:
198*>
199*>
200*>           RFP A                   RFP A
201*>
202*>     -- -- --                   -- -- -- -- -- --
203*>     02 12 22 00 01             00 10 20 30 40 50
204*>     -- -- -- --                   -- -- -- -- --
205*>     03 13 23 33 11             33 11 21 31 41 51
206*>     -- -- -- -- --                   -- -- -- --
207*>     04 14 24 34 44             43 44 22 32 42 52
208*> \endverbatim
209*>
210*  =====================================================================
211      SUBROUTINE CPFTRI( TRANSR, UPLO, N, A, INFO )
212*
213*  -- LAPACK computational routine --
214*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
215*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217*     .. Scalar Arguments ..
218      CHARACTER          TRANSR, UPLO
219      INTEGER            INFO, N
220*     .. Array Arguments ..
221      COMPLEX            A( 0: * )
222*     ..
223*
224*  =====================================================================
225*
226*     .. Parameters ..
227      REAL               ONE
228      COMPLEX            CONE
229      PARAMETER          ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) )
230*     ..
231*     .. Local Scalars ..
232      LOGICAL            LOWER, NISODD, NORMALTRANSR
233      INTEGER            N1, N2, K
234*     ..
235*     .. External Functions ..
236      LOGICAL            LSAME
237      EXTERNAL           LSAME
238*     ..
239*     .. External Subroutines ..
240      EXTERNAL           XERBLA, CTFTRI, CLAUUM, CTRMM, CHERK
241*     ..
242*     .. Intrinsic Functions ..
243      INTRINSIC          MOD
244*     ..
245*     .. Executable Statements ..
246*
247*     Test the input parameters.
248*
249      INFO = 0
250      NORMALTRANSR = LSAME( TRANSR, 'N' )
251      LOWER = LSAME( UPLO, 'L' )
252      IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
253         INFO = -1
254      ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
255         INFO = -2
256      ELSE IF( N.LT.0 ) THEN
257         INFO = -3
258      END IF
259      IF( INFO.NE.0 ) THEN
260         CALL XERBLA( 'CPFTRI', -INFO )
261         RETURN
262      END IF
263*
264*     Quick return if possible
265*
266      IF( N.EQ.0 )
267     $   RETURN
268*
269*     Invert the triangular Cholesky factor U or L.
270*
271      CALL CTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
272      IF( INFO.GT.0 )
273     $   RETURN
274*
275*     If N is odd, set NISODD = .TRUE.
276*     If N is even, set K = N/2 and NISODD = .FALSE.
277*
278      IF( MOD( N, 2 ).EQ.0 ) THEN
279         K = N / 2
280         NISODD = .FALSE.
281      ELSE
282         NISODD = .TRUE.
283      END IF
284*
285*     Set N1 and N2 depending on LOWER
286*
287      IF( LOWER ) THEN
288         N2 = N / 2
289         N1 = N - N2
290      ELSE
291         N1 = N / 2
292         N2 = N - N1
293      END IF
294*
295*     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
296*     inv(L)^C*inv(L). There are eight cases.
297*
298      IF( NISODD ) THEN
299*
300*        N is odd
301*
302         IF( NORMALTRANSR ) THEN
303*
304*           N is odd and TRANSR = 'N'
305*
306            IF( LOWER ) THEN
307*
308*              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
309*              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
310*              T1 -> a(0), T2 -> a(n), S -> a(N1)
311*
312               CALL CLAUUM( 'L', N1, A( 0 ), N, INFO )
313               CALL CHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE,
314     $                     A( 0 ), N )
315               CALL CTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N,
316     $                     A( N1 ), N )
317               CALL CLAUUM( 'U', N2, A( N ), N, INFO )
318*
319            ELSE
320*
321*              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
322*              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
323*              T1 -> a(N2), T2 -> a(N1), S -> a(0)
324*
325               CALL CLAUUM( 'L', N1, A( N2 ), N, INFO )
326               CALL CHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
327     $                     A( N2 ), N )
328               CALL CTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N,
329     $                     A( 0 ), N )
330               CALL CLAUUM( 'U', N2, A( N1 ), N, INFO )
331*
332            END IF
333*
334         ELSE
335*
336*           N is odd and TRANSR = 'C'
337*
338            IF( LOWER ) THEN
339*
340*              SRPA for LOWER, TRANSPOSE, and N is odd
341*              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
342*
343               CALL CLAUUM( 'U', N1, A( 0 ), N1, INFO )
344               CALL CHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
345     $                     A( 0 ), N1 )
346               CALL CTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1,
347     $                     A( N1*N1 ), N1 )
348               CALL CLAUUM( 'L', N2, A( 1 ), N1, INFO )
349*
350            ELSE
351*
352*              SRPA for UPPER, TRANSPOSE, and N is odd
353*              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
354*
355               CALL CLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
356               CALL CHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE,
357     $                     A( N2*N2 ), N2 )
358               CALL CTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ),
359     $                     N2, A( 0 ), N2 )
360               CALL CLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
361*
362            END IF
363*
364         END IF
365*
366      ELSE
367*
368*        N is even
369*
370         IF( NORMALTRANSR ) THEN
371*
372*           N is even and TRANSR = 'N'
373*
374            IF( LOWER ) THEN
375*
376*              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
377*              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
378*              T1 -> a(1), T2 -> a(0), S -> a(k+1)
379*
380               CALL CLAUUM( 'L', K, A( 1 ), N+1, INFO )
381               CALL CHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE,
382     $                     A( 1 ), N+1 )
383               CALL CTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1,
384     $                     A( K+1 ), N+1 )
385               CALL CLAUUM( 'U', K, A( 0 ), N+1, INFO )
386*
387            ELSE
388*
389*              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
390*              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
391*              T1 -> a(k+1), T2 -> a(k), S -> a(0)
392*
393               CALL CLAUUM( 'L', K, A( K+1 ), N+1, INFO )
394               CALL CHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
395     $                     A( K+1 ), N+1 )
396               CALL CTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1,
397     $                     A( 0 ), N+1 )
398               CALL CLAUUM( 'U', K, A( K ), N+1, INFO )
399*
400            END IF
401*
402         ELSE
403*
404*           N is even and TRANSR = 'C'
405*
406            IF( LOWER ) THEN
407*
408*              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
409*              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
410*              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
411*
412               CALL CLAUUM( 'U', K, A( K ), K, INFO )
413               CALL CHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
414     $                     A( K ), K )
415               CALL CTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K,
416     $                     A( K*( K+1 ) ), K )
417               CALL CLAUUM( 'L', K, A( 0 ), K, INFO )
418*
419            ELSE
420*
421*              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
422*              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
423*              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
424*
425               CALL CLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
426               CALL CHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE,
427     $                     A( K*( K+1 ) ), K )
428               CALL CTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K,
429     $                     A( 0 ), K )
430               CALL CLAUUM( 'L', K, A( K*K ), K, INFO )
431*
432            END IF
433*
434         END IF
435*
436      END IF
437*
438      RETURN
439*
440*     End of CPFTRI
441*
442      END
443