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*> \date November 2011
101*
102*> \ingroup complexOTHERcomputational
103*
104*> \par Further Details:
105*  =====================
106*>
107*> \verbatim
108*>
109*>  We first consider Standard Packed Format when N is even.
110*>  We give an example where N = 6.
111*>
112*>      AP is Upper             AP is Lower
113*>
114*>   00 01 02 03 04 05       00
115*>      11 12 13 14 15       10 11
116*>         22 23 24 25       20 21 22
117*>            33 34 35       30 31 32 33
118*>               44 45       40 41 42 43 44
119*>                  55       50 51 52 53 54 55
120*>
121*>
122*>  Let TRANSR = 'N'. RFP holds AP as follows:
123*>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
124*>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
125*>  conjugate-transpose of the first three columns of AP upper.
126*>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
127*>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
128*>  conjugate-transpose of the last three columns of AP lower.
129*>  To denote conjugate we place -- above the element. This covers the
130*>  case N even and TRANSR = 'N'.
131*>
132*>         RFP A                   RFP A
133*>
134*>                                -- -- --
135*>        03 04 05                33 43 53
136*>                                   -- --
137*>        13 14 15                00 44 54
138*>                                      --
139*>        23 24 25                10 11 55
140*>
141*>        33 34 35                20 21 22
142*>        --
143*>        00 44 45                30 31 32
144*>        -- --
145*>        01 11 55                40 41 42
146*>        -- -- --
147*>        02 12 22                50 51 52
148*>
149*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
150*>  transpose of RFP A above. One therefore gets:
151*>
152*>
153*>           RFP A                   RFP A
154*>
155*>     -- -- -- --                -- -- -- -- -- --
156*>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
157*>     -- -- -- -- --                -- -- -- -- --
158*>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
159*>     -- -- -- -- -- --                -- -- -- --
160*>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
161*>
162*>
163*>  We next  consider Standard Packed Format when N is odd.
164*>  We give an example where N = 5.
165*>
166*>     AP is Upper                 AP is Lower
167*>
168*>   00 01 02 03 04              00
169*>      11 12 13 14              10 11
170*>         22 23 24              20 21 22
171*>            33 34              30 31 32 33
172*>               44              40 41 42 43 44
173*>
174*>
175*>  Let TRANSR = 'N'. RFP holds AP as follows:
176*>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
177*>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
178*>  conjugate-transpose of the first two   columns of AP upper.
179*>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
180*>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
181*>  conjugate-transpose of the last two   columns of AP lower.
182*>  To denote conjugate we place -- above the element. This covers the
183*>  case N odd  and TRANSR = 'N'.
184*>
185*>         RFP A                   RFP A
186*>
187*>                                   -- --
188*>        02 03 04                00 33 43
189*>                                      --
190*>        12 13 14                10 11 44
191*>
192*>        22 23 24                20 21 22
193*>        --
194*>        00 33 34                30 31 32
195*>        -- --
196*>        01 11 44                40 41 42
197*>
198*>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
199*>  transpose of RFP A above. One therefore gets:
200*>
201*>
202*>           RFP A                   RFP A
203*>
204*>     -- -- --                   -- -- -- -- -- --
205*>     02 12 22 00 01             00 10 20 30 40 50
206*>     -- -- -- --                   -- -- -- -- --
207*>     03 13 23 33 11             33 11 21 31 41 51
208*>     -- -- -- -- --                   -- -- -- --
209*>     04 14 24 34 44             43 44 22 32 42 52
210*> \endverbatim
211*>
212*  =====================================================================
213      SUBROUTINE CPFTRI( TRANSR, UPLO, N, A, INFO )
214*
215*  -- LAPACK computational routine (version 3.4.0) --
216*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
217*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218*     November 2011
219*
220*     .. Scalar Arguments ..
221      CHARACTER          TRANSR, UPLO
222      INTEGER            INFO, N
223*     .. Array Arguments ..
224      COMPLEX            A( 0: * )
225*     ..
226*
227*  =====================================================================
228*
229*     .. Parameters ..
230      REAL               ONE
231      COMPLEX            CONE
232      PARAMETER          ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) )
233*     ..
234*     .. Local Scalars ..
235      LOGICAL            LOWER, NISODD, NORMALTRANSR
236      INTEGER            N1, N2, K
237*     ..
238*     .. External Functions ..
239      LOGICAL            LSAME
240      EXTERNAL           LSAME
241*     ..
242*     .. External Subroutines ..
243      EXTERNAL           XERBLA, CTFTRI, CLAUUM, CTRMM, CHERK
244*     ..
245*     .. Intrinsic Functions ..
246      INTRINSIC          MOD
247*     ..
248*     .. Executable Statements ..
249*
250*     Test the input parameters.
251*
252      INFO = 0
253      NORMALTRANSR = LSAME( TRANSR, 'N' )
254      LOWER = LSAME( UPLO, 'L' )
255      IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
256         INFO = -1
257      ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
258         INFO = -2
259      ELSE IF( N.LT.0 ) THEN
260         INFO = -3
261      END IF
262      IF( INFO.NE.0 ) THEN
263         CALL XERBLA( 'CPFTRI', -INFO )
264         RETURN
265      END IF
266*
267*     Quick return if possible
268*
269      IF( N.EQ.0 )
270     $   RETURN
271*
272*     Invert the triangular Cholesky factor U or L.
273*
274      CALL CTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
275      IF( INFO.GT.0 )
276     $   RETURN
277*
278*     If N is odd, set NISODD = .TRUE.
279*     If N is even, set K = N/2 and NISODD = .FALSE.
280*
281      IF( MOD( N, 2 ).EQ.0 ) THEN
282         K = N / 2
283         NISODD = .FALSE.
284      ELSE
285         NISODD = .TRUE.
286      END IF
287*
288*     Set N1 and N2 depending on LOWER
289*
290      IF( LOWER ) THEN
291         N2 = N / 2
292         N1 = N - N2
293      ELSE
294         N1 = N / 2
295         N2 = N - N1
296      END IF
297*
298*     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
299*     inv(L)^C*inv(L). There are eight cases.
300*
301      IF( NISODD ) THEN
302*
303*        N is odd
304*
305         IF( NORMALTRANSR ) THEN
306*
307*           N is odd and TRANSR = 'N'
308*
309            IF( LOWER ) THEN
310*
311*              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
312*              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
313*              T1 -> a(0), T2 -> a(n), S -> a(N1)
314*
315               CALL CLAUUM( 'L', N1, A( 0 ), N, INFO )
316               CALL CHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE,
317     $                     A( 0 ), N )
318               CALL CTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N,
319     $                     A( N1 ), N )
320               CALL CLAUUM( 'U', N2, A( N ), N, INFO )
321*
322            ELSE
323*
324*              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
325*              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
326*              T1 -> a(N2), T2 -> a(N1), S -> a(0)
327*
328               CALL CLAUUM( 'L', N1, A( N2 ), N, INFO )
329               CALL CHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
330     $                     A( N2 ), N )
331               CALL CTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N,
332     $                     A( 0 ), N )
333               CALL CLAUUM( 'U', N2, A( N1 ), N, INFO )
334*
335            END IF
336*
337         ELSE
338*
339*           N is odd and TRANSR = 'C'
340*
341            IF( LOWER ) THEN
342*
343*              SRPA for LOWER, TRANSPOSE, and N is odd
344*              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
345*
346               CALL CLAUUM( 'U', N1, A( 0 ), N1, INFO )
347               CALL CHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
348     $                     A( 0 ), N1 )
349               CALL CTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1,
350     $                     A( N1*N1 ), N1 )
351               CALL CLAUUM( 'L', N2, A( 1 ), N1, INFO )
352*
353            ELSE
354*
355*              SRPA for UPPER, TRANSPOSE, and N is odd
356*              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
357*
358               CALL CLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
359               CALL CHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE,
360     $                     A( N2*N2 ), N2 )
361               CALL CTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ),
362     $                     N2, A( 0 ), N2 )
363               CALL CLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
364*
365            END IF
366*
367         END IF
368*
369      ELSE
370*
371*        N is even
372*
373         IF( NORMALTRANSR ) THEN
374*
375*           N is even and TRANSR = 'N'
376*
377            IF( LOWER ) THEN
378*
379*              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
380*              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
381*              T1 -> a(1), T2 -> a(0), S -> a(k+1)
382*
383               CALL CLAUUM( 'L', K, A( 1 ), N+1, INFO )
384               CALL CHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE,
385     $                     A( 1 ), N+1 )
386               CALL CTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1,
387     $                     A( K+1 ), N+1 )
388               CALL CLAUUM( 'U', K, A( 0 ), N+1, INFO )
389*
390            ELSE
391*
392*              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
393*              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
394*              T1 -> a(k+1), T2 -> a(k), S -> a(0)
395*
396               CALL CLAUUM( 'L', K, A( K+1 ), N+1, INFO )
397               CALL CHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
398     $                     A( K+1 ), N+1 )
399               CALL CTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1,
400     $                     A( 0 ), N+1 )
401               CALL CLAUUM( 'U', K, A( K ), N+1, INFO )
402*
403            END IF
404*
405         ELSE
406*
407*           N is even and TRANSR = 'C'
408*
409            IF( LOWER ) THEN
410*
411*              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
412*              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
413*              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
414*
415               CALL CLAUUM( 'U', K, A( K ), K, INFO )
416               CALL CHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
417     $                     A( K ), K )
418               CALL CTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K,
419     $                     A( K*( K+1 ) ), K )
420               CALL CLAUUM( 'L', K, A( 0 ), K, INFO )
421*
422            ELSE
423*
424*              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
425*              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
426*              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
427*
428               CALL CLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
429               CALL CHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE,
430     $                     A( K*( K+1 ) ), K )
431               CALL CTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K,
432     $                     A( 0 ), K )
433               CALL CLAUUM( 'L', K, A( K*K ), K, INFO )
434*
435            END IF
436*
437         END IF
438*
439      END IF
440*
441      RETURN
442*
443*     End of CPFTRI
444*
445      END
446