1 SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 *
11 * ..
12 * .. Scalar Arguments ..
13 CHARACTER TRANSR, UPLO
14 INTEGER N, INFO
15 * ..
16 * .. Array Arguments ..
17 DOUBLE PRECISION A( 0: * )
18 *
19 * Purpose
20 * =======
21 *
22 * DPFTRF computes the Cholesky factorization of a real symmetric
23 * positive definite matrix A.
24 *
25 * The factorization has the form
26 * A = U**T * U, if UPLO = 'U', or
27 * A = L * L**T, if UPLO = 'L',
28 * where U is an upper triangular matrix and L is lower triangular.
29 *
30 * This is the block version of the algorithm, calling Level 3 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * TRANSR (input) CHARACTER*1
36 * = 'N': The Normal TRANSR of RFP A is stored;
37 * = 'T': The Transpose TRANSR of RFP A is stored.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of RFP A is stored;
41 * = 'L': Lower triangle of RFP A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * A (input/output) DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
47 * On entry, the symmetric matrix A in RFP format. RFP format is
48 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
49 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
50 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
51 * the transpose of RFP A as defined when
52 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
53 * follows: If UPLO = 'U' the RFP A contains the NT elements of
54 * upper packed A. If UPLO = 'L' the RFP A contains the elements
55 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
56 * 'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
57 * is odd. See the Note below for more details.
58 *
59 * On exit, if INFO = 0, the factor U or L from the Cholesky
60 * factorization RFP A = U**T*U or RFP A = L*L**T.
61 *
62 * INFO (output) INTEGER
63 * = 0: successful exit
64 * < 0: if INFO = -i, the i-th argument had an illegal value
65 * > 0: if INFO = i, the leading minor of order i is not
66 * positive definite, and the factorization could not be
67 * completed.
68 *
69 * Further Details
70 * ===============
71 *
72 * We first consider Rectangular Full Packed (RFP) Format when N is
73 * even. We give an example where N = 6.
74 *
75 * AP is Upper AP is Lower
76 *
77 * 00 01 02 03 04 05 00
78 * 11 12 13 14 15 10 11
79 * 22 23 24 25 20 21 22
80 * 33 34 35 30 31 32 33
81 * 44 45 40 41 42 43 44
82 * 55 50 51 52 53 54 55
83 *
84 *
85 * Let TRANSR = 'N'. RFP holds AP as follows:
86 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
87 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
88 * the transpose of the first three columns of AP upper.
89 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
90 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
91 * the transpose of the last three columns of AP lower.
92 * This covers the case N even and TRANSR = 'N'.
93 *
94 * RFP A RFP A
95 *
96 * 03 04 05 33 43 53
97 * 13 14 15 00 44 54
98 * 23 24 25 10 11 55
99 * 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 then 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 DOUBLE PRECISION ONE
157 PARAMETER ( ONE = 1.0D+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, DSYRK, DPOTRF, DTRSM
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( 'DPFTRF', -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 DPOTRF( 'L', N1, A( 0 ), N, INFO )
234 IF( INFO.GT.0 )
235 $ RETURN
236 CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, A( 0 ), N,
237 $ A( N1 ), N )
238 CALL DSYRK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE,
239 $ A( N ), N )
240 CALL DPOTRF( '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 DPOTRF( 'L', N1, A( N2 ), N, INFO )
251 IF( INFO.GT.0 )
252 $ RETURN
253 CALL DTRSM( 'L', 'L', 'N', 'N', N1, N2, ONE, A( N2 ), N,
254 $ A( 0 ), N )
255 CALL DSYRK( 'U', 'T', N2, N1, -ONE, A( 0 ), N, ONE,
256 $ A( N1 ), N )
257 CALL DPOTRF( '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 DPOTRF( 'U', N1, A( 0 ), N1, INFO )
274 IF( INFO.GT.0 )
275 $ RETURN
276 CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, A( 0 ), N1,
277 $ A( N1*N1 ), N1 )
278 CALL DSYRK( 'L', 'T', N2, N1, -ONE, A( N1*N1 ), N1, ONE,
279 $ A( 1 ), N1 )
280 CALL DPOTRF( '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 DPOTRF( 'U', N1, A( N2*N2 ), N2, INFO )
291 IF( INFO.GT.0 )
292 $ RETURN
293 CALL DTRSM( 'R', 'U', 'N', 'N', N2, N1, ONE, A( N2*N2 ),
294 $ N2, A( 0 ), N2 )
295 CALL DSYRK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE,
296 $ A( N1*N2 ), N2 )
297 CALL DPOTRF( '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 DPOTRF( 'L', K, A( 1 ), N+1, INFO )
320 IF( INFO.GT.0 )
321 $ RETURN
322 CALL DTRSM( 'R', 'L', 'T', 'N', K, K, ONE, A( 1 ), N+1,
323 $ A( K+1 ), N+1 )
324 CALL DSYRK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE,
325 $ A( 0 ), N+1 )
326 CALL DPOTRF( '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 DPOTRF( 'L', K, A( K+1 ), N+1, INFO )
337 IF( INFO.GT.0 )
338 $ RETURN
339 CALL DTRSM( 'L', 'L', 'N', 'N', K, K, ONE, A( K+1 ),
340 $ N+1, A( 0 ), N+1 )
341 CALL DSYRK( 'U', 'T', K, K, -ONE, A( 0 ), N+1, ONE,
342 $ A( K ), N+1 )
343 CALL DPOTRF( '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 DPOTRF( 'U', K, A( 0+K ), K, INFO )
360 IF( INFO.GT.0 )
361 $ RETURN
362 CALL DTRSM( 'L', 'U', 'T', 'N', K, K, ONE, A( K ), N1,
363 $ A( K*( K+1 ) ), K )
364 CALL DSYRK( 'L', 'T', K, K, -ONE, A( K*( K+1 ) ), K, ONE,
365 $ A( 0 ), K )
366 CALL DPOTRF( '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 DPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO )
377 IF( INFO.GT.0 )
378 $ RETURN
379 CALL DTRSM( 'R', 'U', 'N', 'N', K, K, ONE,
380 $ A( K*( K+1 ) ), K, A( 0 ), K )
381 CALL DSYRK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE,
382 $ A( K*K ), K )
383 CALL DPOTRF( '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 DPFTRF
396 *
397 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 *
5 * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6 * -- April 2011 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10 *
11 * ..
12 * .. Scalar Arguments ..
13 CHARACTER TRANSR, UPLO
14 INTEGER N, INFO
15 * ..
16 * .. Array Arguments ..
17 DOUBLE PRECISION A( 0: * )
18 *
19 * Purpose
20 * =======
21 *
22 * DPFTRF computes the Cholesky factorization of a real symmetric
23 * positive definite matrix A.
24 *
25 * The factorization has the form
26 * A = U**T * U, if UPLO = 'U', or
27 * A = L * L**T, if UPLO = 'L',
28 * where U is an upper triangular matrix and L is lower triangular.
29 *
30 * This is the block version of the algorithm, calling Level 3 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * TRANSR (input) CHARACTER*1
36 * = 'N': The Normal TRANSR of RFP A is stored;
37 * = 'T': The Transpose TRANSR of RFP A is stored.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of RFP A is stored;
41 * = 'L': Lower triangle of RFP A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * A (input/output) DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
47 * On entry, the symmetric matrix A in RFP format. RFP format is
48 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
49 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
50 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
51 * the transpose of RFP A as defined when
52 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
53 * follows: If UPLO = 'U' the RFP A contains the NT elements of
54 * upper packed A. If UPLO = 'L' the RFP A contains the elements
55 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
56 * 'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
57 * is odd. See the Note below for more details.
58 *
59 * On exit, if INFO = 0, the factor U or L from the Cholesky
60 * factorization RFP A = U**T*U or RFP A = L*L**T.
61 *
62 * INFO (output) INTEGER
63 * = 0: successful exit
64 * < 0: if INFO = -i, the i-th argument had an illegal value
65 * > 0: if INFO = i, the leading minor of order i is not
66 * positive definite, and the factorization could not be
67 * completed.
68 *
69 * Further Details
70 * ===============
71 *
72 * We first consider Rectangular Full Packed (RFP) Format when N is
73 * even. We give an example where N = 6.
74 *
75 * AP is Upper AP is Lower
76 *
77 * 00 01 02 03 04 05 00
78 * 11 12 13 14 15 10 11
79 * 22 23 24 25 20 21 22
80 * 33 34 35 30 31 32 33
81 * 44 45 40 41 42 43 44
82 * 55 50 51 52 53 54 55
83 *
84 *
85 * Let TRANSR = 'N'. RFP holds AP as follows:
86 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
87 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
88 * the transpose of the first three columns of AP upper.
89 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
90 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
91 * the transpose of the last three columns of AP lower.
92 * This covers the case N even and TRANSR = 'N'.
93 *
94 * RFP A RFP A
95 *
96 * 03 04 05 33 43 53
97 * 13 14 15 00 44 54
98 * 23 24 25 10 11 55
99 * 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 then 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 DOUBLE PRECISION ONE
157 PARAMETER ( ONE = 1.0D+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, DSYRK, DPOTRF, DTRSM
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( 'DPFTRF', -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 DPOTRF( 'L', N1, A( 0 ), N, INFO )
234 IF( INFO.GT.0 )
235 $ RETURN
236 CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, A( 0 ), N,
237 $ A( N1 ), N )
238 CALL DSYRK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE,
239 $ A( N ), N )
240 CALL DPOTRF( '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 DPOTRF( 'L', N1, A( N2 ), N, INFO )
251 IF( INFO.GT.0 )
252 $ RETURN
253 CALL DTRSM( 'L', 'L', 'N', 'N', N1, N2, ONE, A( N2 ), N,
254 $ A( 0 ), N )
255 CALL DSYRK( 'U', 'T', N2, N1, -ONE, A( 0 ), N, ONE,
256 $ A( N1 ), N )
257 CALL DPOTRF( '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 DPOTRF( 'U', N1, A( 0 ), N1, INFO )
274 IF( INFO.GT.0 )
275 $ RETURN
276 CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, A( 0 ), N1,
277 $ A( N1*N1 ), N1 )
278 CALL DSYRK( 'L', 'T', N2, N1, -ONE, A( N1*N1 ), N1, ONE,
279 $ A( 1 ), N1 )
280 CALL DPOTRF( '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 DPOTRF( 'U', N1, A( N2*N2 ), N2, INFO )
291 IF( INFO.GT.0 )
292 $ RETURN
293 CALL DTRSM( 'R', 'U', 'N', 'N', N2, N1, ONE, A( N2*N2 ),
294 $ N2, A( 0 ), N2 )
295 CALL DSYRK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE,
296 $ A( N1*N2 ), N2 )
297 CALL DPOTRF( '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 DPOTRF( 'L', K, A( 1 ), N+1, INFO )
320 IF( INFO.GT.0 )
321 $ RETURN
322 CALL DTRSM( 'R', 'L', 'T', 'N', K, K, ONE, A( 1 ), N+1,
323 $ A( K+1 ), N+1 )
324 CALL DSYRK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE,
325 $ A( 0 ), N+1 )
326 CALL DPOTRF( '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 DPOTRF( 'L', K, A( K+1 ), N+1, INFO )
337 IF( INFO.GT.0 )
338 $ RETURN
339 CALL DTRSM( 'L', 'L', 'N', 'N', K, K, ONE, A( K+1 ),
340 $ N+1, A( 0 ), N+1 )
341 CALL DSYRK( 'U', 'T', K, K, -ONE, A( 0 ), N+1, ONE,
342 $ A( K ), N+1 )
343 CALL DPOTRF( '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 DPOTRF( 'U', K, A( 0+K ), K, INFO )
360 IF( INFO.GT.0 )
361 $ RETURN
362 CALL DTRSM( 'L', 'U', 'T', 'N', K, K, ONE, A( K ), N1,
363 $ A( K*( K+1 ) ), K )
364 CALL DSYRK( 'L', 'T', K, K, -ONE, A( K*( K+1 ) ), K, ONE,
365 $ A( 0 ), K )
366 CALL DPOTRF( '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 DPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO )
377 IF( INFO.GT.0 )
378 $ RETURN
379 CALL DTRSM( 'R', 'U', 'N', 'N', K, K, ONE,
380 $ A( K*( K+1 ) ), K, A( 0 ), K )
381 CALL DSYRK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE,
382 $ A( K*K ), K )
383 CALL DPOTRF( '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 DPFTRF
396 *
397 END