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