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