001:       SUBROUTINE CPFTRF( TRANSR, UPLO, N, A, INFO )
002: *
003: *  -- LAPACK routine (version 3.2)                                    --
004: *
005: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
006: *  -- November 2008                                                   --
007: *
008: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
009: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
010: *
011: *     ..
012: *     .. Scalar Arguments ..
013:       CHARACTER          TRANSR, UPLO
014:       INTEGER            N, INFO
015: *     ..
016: *     .. Array Arguments ..
017:       COMPLEX            A( 0: * )
018: *
019: *  Purpose
020: *  =======
021: *
022: *  CPFTRF computes the Cholesky factorization of a complex Hermitian
023: *  positive definite matrix A.
024: *
025: *  The factorization has the form
026: *     A = U**H * U,  if UPLO = 'U', or
027: *     A = L  * L**H,  if UPLO = 'L',
028: *  where U is an upper triangular matrix and L is lower triangular.
029: *
030: *  This is the block version of the algorithm, calling Level 3 BLAS.
031: *
032: *  Arguments
033: *  =========
034: *
035: *  TRANSR    (input) CHARACTER
036: *          = 'N':  The Normal TRANSR of RFP A is stored;
037: *          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.
038: *
039: *  UPLO    (input) CHARACTER
040: *          = 'U':  Upper triangle of RFP A is stored;
041: *          = 'L':  Lower triangle of RFP A is stored.
042: *
043: *  N       (input) INTEGER
044: *          The order of the matrix A.  N >= 0.
045: *
046: *  A       (input/output) COMPLEX array, dimension ( N*(N+1)/2 );
047: *          On entry, the Hermitian matrix A in RFP format. RFP format is
048: *          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
049: *          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
050: *          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
051: *          the Conjugate-transpose of RFP A as defined when
052: *          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
053: *          follows: If UPLO = 'U' the RFP A contains the nt elements of
054: *          upper packed A. If UPLO = 'L' the RFP A contains the elements
055: *          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
056: *          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
057: *          is odd. See the Note below for more details.
058: *
059: *          On exit, if INFO = 0, the factor U or L from the Cholesky
060: *          factorization RFP A = U**H*U or RFP A = L*L**H.
061: *
062: *  INFO    (output) INTEGER
063: *          = 0:  successful exit
064: *          < 0:  if INFO = -i, the i-th argument had an illegal value
065: *          > 0:  if INFO = i, the leading minor of order i is not
066: *                positive definite, and the factorization could not be
067: *                completed.
068: *
069: *  Further Notes on RFP Format:
070: *  ============================
071: *
072: *
073: *  We first consider Standard Packed Format when N is even.
074: *  We give an example where N = 6.
075: *
076: *     AP is Upper             AP is Lower
077: *
078: *   00 01 02 03 04 05       00
079: *      11 12 13 14 15       10 11
080: *         22 23 24 25       20 21 22
081: *            33 34 35       30 31 32 33
082: *               44 45       40 41 42 43 44
083: *                  55       50 51 52 53 54 55
084: *
085: *
086: *  Let TRANSR = 'N'. RFP holds AP as follows:
087: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
088: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
089: *  conjugate-transpose of the first three columns of AP upper.
090: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
091: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
092: *  conjugate-transpose of the last three columns of AP lower.
093: *  To denote conjugate we place -- above the element. This covers the
094: *  case N even and TRANSR = 'N'.
095: *
096: *         RFP A                   RFP A
097: *
098: *                                -- -- --
099: *        03 04 05                33 43 53
100: *                                   -- --
101: *        13 14 15                00 44 54
102: *                                      --
103: *        23 24 25                10 11 55
104: *
105: *        33 34 35                20 21 22
106: *        --
107: *        00 44 45                30 31 32
108: *        -- --
109: *        01 11 55                40 41 42
110: *        -- -- --
111: *        02 12 22                50 51 52
112: *
113: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
114: *  transpose of RFP A above. One therefore gets:
115: *
116: *
117: *           RFP A                   RFP A
118: *
119: *     -- -- -- --                -- -- -- -- -- --
120: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
121: *     -- -- -- -- --                -- -- -- -- --
122: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
123: *     -- -- -- -- -- --                -- -- -- --
124: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
125: *
126: *
127: *  We next  consider Standard Packed Format when N is odd.
128: *  We give an example where N = 5.
129: *
130: *     AP is Upper                 AP is Lower
131: *
132: *   00 01 02 03 04              00
133: *      11 12 13 14              10 11
134: *         22 23 24              20 21 22
135: *            33 34              30 31 32 33
136: *               44              40 41 42 43 44
137: *
138: *
139: *  Let TRANSR = 'N'. RFP holds AP as follows:
140: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
141: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
142: *  conjugate-transpose of the first two   columns of AP upper.
143: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
144: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
145: *  conjugate-transpose of the last two   columns of AP lower.
146: *  To denote conjugate we place -- above the element. This covers the
147: *  case N odd  and TRANSR = 'N'.
148: *
149: *         RFP A                   RFP A
150: *
151: *                                   -- --
152: *        02 03 04                00 33 43
153: *                                      --
154: *        12 13 14                10 11 44
155: *
156: *        22 23 24                20 21 22
157: *        --
158: *        00 33 34                30 31 32
159: *        -- --
160: *        01 11 44                40 41 42
161: *
162: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
163: *  transpose of RFP A above. One therefore gets:
164: *
165: *
166: *           RFP A                   RFP A
167: *
168: *     -- -- --                   -- -- -- -- -- --
169: *     02 12 22 00 01             00 10 20 30 40 50
170: *     -- -- -- --                   -- -- -- -- --
171: *     03 13 23 33 11             33 11 21 31 41 51
172: *     -- -- -- -- --                   -- -- -- --
173: *     04 14 24 34 44             43 44 22 32 42 52
174: *
175: *  =====================================================================
176: *
177: *     .. Parameters ..
178:       REAL               ONE
179:       COMPLEX            CONE
180:       PARAMETER          ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) )
181: *     ..
182: *     .. Local Scalars ..
183:       LOGICAL            LOWER, NISODD, NORMALTRANSR
184:       INTEGER            N1, N2, K
185: *     ..
186: *     .. External Functions ..
187:       LOGICAL            LSAME
188:       EXTERNAL           LSAME
189: *     ..
190: *     .. External Subroutines ..
191:       EXTERNAL           XERBLA, CHERK, CPOTRF, CTRSM
192: *     ..
193: *     .. Intrinsic Functions ..
194:       INTRINSIC          MOD
195: *     ..
196: *     .. Executable Statements ..
197: *
198: *     Test the input parameters.
199: *
200:       INFO = 0
201:       NORMALTRANSR = LSAME( TRANSR, 'N' )
202:       LOWER = LSAME( UPLO, 'L' )
203:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
204:          INFO = -1
205:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
206:          INFO = -2
207:       ELSE IF( N.LT.0 ) THEN
208:          INFO = -3
209:       END IF
210:       IF( INFO.NE.0 ) THEN
211:          CALL XERBLA( 'CPFTRF', -INFO )
212:          RETURN
213:       END IF
214: *
215: *     Quick return if possible
216: *
217:       IF( N.EQ.0 )
218:      +   RETURN
219: *
220: *     If N is odd, set NISODD = .TRUE.
221: *     If N is even, set K = N/2 and NISODD = .FALSE.
222: *
223:       IF( MOD( N, 2 ).EQ.0 ) THEN
224:          K = N / 2
225:          NISODD = .FALSE.
226:       ELSE
227:          NISODD = .TRUE.
228:       END IF
229: *
230: *     Set N1 and N2 depending on LOWER
231: *
232:       IF( LOWER ) THEN
233:          N2 = N / 2
234:          N1 = N - N2
235:       ELSE
236:          N1 = N / 2
237:          N2 = N - N1
238:       END IF
239: *
240: *     start execution: there are eight cases
241: *
242:       IF( NISODD ) THEN
243: *
244: *        N is odd
245: *
246:          IF( NORMALTRANSR ) THEN
247: *
248: *           N is odd and TRANSR = 'N'
249: *
250:             IF( LOWER ) THEN
251: *
252: *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
253: *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
254: *             T1 -> a(0), T2 -> a(n), S -> a(n1)
255: *
256:                CALL CPOTRF( 'L', N1, A( 0 ), N, INFO )
257:                IF( INFO.GT.0 )
258:      +            RETURN
259:                CALL CTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, A( 0 ), N,
260:      +                     A( N1 ), N )
261:                CALL CHERK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE,
262:      +                     A( N ), N )
263:                CALL CPOTRF( 'U', N2, A( N ), N, INFO )
264:                IF( INFO.GT.0 )
265:      +            INFO = INFO + N1
266: *
267:             ELSE
268: *
269: *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
270: *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
271: *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
272: *
273:                CALL CPOTRF( 'L', N1, A( N2 ), N, INFO )
274:                IF( INFO.GT.0 )
275:      +            RETURN
276:                CALL CTRSM( 'L', 'L', 'N', 'N', N1, N2, CONE, A( N2 ), N,
277:      +                     A( 0 ), N )
278:                CALL CHERK( 'U', 'C', N2, N1, -ONE, A( 0 ), N, ONE,
279:      +                     A( N1 ), N )
280:                CALL CPOTRF( 'U', N2, A( N1 ), N, INFO )
281:                IF( INFO.GT.0 )
282:      +            INFO = INFO + N1
283: *
284:             END IF
285: *
286:          ELSE
287: *
288: *           N is odd and TRANSR = 'C'
289: *
290:             IF( LOWER ) THEN
291: *
292: *              SRPA for LOWER, TRANSPOSE and N is odd
293: *              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
294: *              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
295: *
296:                CALL CPOTRF( 'U', N1, A( 0 ), N1, INFO )
297:                IF( INFO.GT.0 )
298:      +            RETURN
299:                CALL CTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, A( 0 ), N1,
300:      +                     A( N1*N1 ), N1 )
301:                CALL CHERK( 'L', 'C', N2, N1, -ONE, A( N1*N1 ), N1, ONE,
302:      +                     A( 1 ), N1 )
303:                CALL CPOTRF( 'L', N2, A( 1 ), N1, INFO )
304:                IF( INFO.GT.0 )
305:      +            INFO = INFO + N1
306: *
307:             ELSE
308: *
309: *              SRPA for UPPER, TRANSPOSE and N is odd
310: *              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
311: *              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
312: *
313:                CALL CPOTRF( 'U', N1, A( N2*N2 ), N2, INFO )
314:                IF( INFO.GT.0 )
315:      +            RETURN
316:                CALL CTRSM( 'R', 'U', 'N', 'N', N2, N1, CONE, A( N2*N2 ),
317:      +                     N2, A( 0 ), N2 )
318:                CALL CHERK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE,
319:      +                     A( N1*N2 ), N2 )
320:                CALL CPOTRF( 'L', N2, A( N1*N2 ), N2, INFO )
321:                IF( INFO.GT.0 )
322:      +            INFO = INFO + N1
323: *
324:             END IF
325: *
326:          END IF
327: *
328:       ELSE
329: *
330: *        N is even
331: *
332:          IF( NORMALTRANSR ) THEN
333: *
334: *           N is even and TRANSR = 'N'
335: *
336:             IF( LOWER ) THEN
337: *
338: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
339: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
340: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
341: *
342:                CALL CPOTRF( 'L', K, A( 1 ), N+1, INFO )
343:                IF( INFO.GT.0 )
344:      +            RETURN
345:                CALL CTRSM( 'R', 'L', 'C', 'N', K, K, CONE, A( 1 ), N+1,
346:      +                     A( K+1 ), N+1 )
347:                CALL CHERK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE,
348:      +                     A( 0 ), N+1 )
349:                CALL CPOTRF( 'U', K, A( 0 ), N+1, INFO )
350:                IF( INFO.GT.0 )
351:      +            INFO = INFO + K
352: *
353:             ELSE
354: *
355: *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
356: *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
357: *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
358: *
359:                CALL CPOTRF( 'L', K, A( K+1 ), N+1, INFO )
360:                IF( INFO.GT.0 )
361:      +            RETURN
362:                CALL CTRSM( 'L', 'L', 'N', 'N', K, K, CONE, A( K+1 ),
363:      +                     N+1, A( 0 ), N+1 )
364:                CALL CHERK( 'U', 'C', K, K, -ONE, A( 0 ), N+1, ONE,
365:      +                     A( K ), N+1 )
366:                CALL CPOTRF( 'U', K, A( K ), N+1, INFO )
367:                IF( INFO.GT.0 )
368:      +            INFO = INFO + K
369: *
370:             END IF
371: *
372:          ELSE
373: *
374: *           N is even and TRANSR = 'C'
375: *
376:             IF( LOWER ) THEN
377: *
378: *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
379: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
380: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
381: *
382:                CALL CPOTRF( 'U', K, A( 0+K ), K, INFO )
383:                IF( INFO.GT.0 )
384:      +            RETURN
385:                CALL CTRSM( 'L', 'U', 'C', 'N', K, K, CONE, A( K ), N1,
386:      +                     A( K*( K+1 ) ), K )
387:                CALL CHERK( 'L', 'C', K, K, -ONE, A( K*( K+1 ) ), K, ONE,
388:      +                     A( 0 ), K )
389:                CALL CPOTRF( 'L', K, A( 0 ), K, INFO )
390:                IF( INFO.GT.0 )
391:      +            INFO = INFO + K
392: *
393:             ELSE
394: *
395: *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
396: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
397: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
398: *
399:                CALL CPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO )
400:                IF( INFO.GT.0 )
401:      +            RETURN
402:                CALL CTRSM( 'R', 'U', 'N', 'N', K, K, CONE,
403:      +                     A( K*( K+1 ) ), K, A( 0 ), K )
404:                CALL CHERK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE,
405:      +                     A( K*K ), K )
406:                CALL CPOTRF( 'L', K, A( K*K ), K, INFO )
407:                IF( INFO.GT.0 )
408:      +            INFO = INFO + K
409: *
410:             END IF
411: *
412:          END IF
413: *
414:       END IF
415: *
416:       RETURN
417: *
418: *     End of CPFTRF
419: *
420:       END
421: