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