 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * 00005 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00006 * -- April 2011 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * 00011 * .. 00012 * .. Scalar Arguments .. 00013 CHARACTER TRANSR, UPLO 00014 INTEGER N, INFO 00015 * .. 00016 * .. Array Arguments .. 00017 DOUBLE PRECISION A( 0: * ) 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DPFTRF computes the Cholesky factorization of a real symmetric 00023 * positive definite matrix A. 00024 * 00025 * The factorization has the form 00026 * A = U**T * U, if UPLO = 'U', or 00027 * A = L * L**T, if UPLO = 'L', 00028 * where U is an upper triangular matrix and L is lower triangular. 00029 * 00030 * This is the block version of the algorithm, calling Level 3 BLAS. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * TRANSR (input) CHARACTER*1 00036 * = 'N': The Normal TRANSR of RFP A is stored; 00037 * = 'T': The Transpose TRANSR of RFP A is stored. 00038 * 00039 * UPLO (input) CHARACTER*1 00040 * = 'U': Upper triangle of RFP A is stored; 00041 * = 'L': Lower triangle of RFP A is stored. 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrix A. N >= 0. 00045 * 00046 * A (input/output) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ); 00047 * On entry, the symmetric matrix A in RFP format. RFP format is 00048 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' 00049 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is 00050 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is 00051 * the transpose of RFP A as defined when 00052 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as 00053 * follows: If UPLO = 'U' the RFP A contains the NT elements of 00054 * upper packed A. If UPLO = 'L' the RFP A contains the elements 00055 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 00056 * 'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N 00057 * is odd. See the Note below for more details. 00058 * 00059 * On exit, if INFO = 0, the factor U or L from the Cholesky 00060 * factorization RFP A = U**T*U or RFP A = L*L**T. 00061 * 00062 * INFO (output) INTEGER 00063 * = 0: successful exit 00064 * < 0: if INFO = -i, the i-th argument had an illegal value 00065 * > 0: if INFO = i, the leading minor of order i is not 00066 * positive definite, and the factorization could not be 00067 * completed. 00068 * 00069 * Further Details 00070 * =============== 00071 * 00072 * We first consider Rectangular Full Packed (RFP) Format when N is 00073 * even. We give an example where N = 6. 00074 * 00075 * AP is Upper AP is Lower 00076 * 00077 * 00 01 02 03 04 05 00 00078 * 11 12 13 14 15 10 11 00079 * 22 23 24 25 20 21 22 00080 * 33 34 35 30 31 32 33 00081 * 44 45 40 41 42 43 44 00082 * 55 50 51 52 53 54 55 00083 * 00084 * 00085 * Let TRANSR = 'N'. RFP holds AP as follows: 00086 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00087 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00088 * the transpose of the first three columns of AP upper. 00089 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00090 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00091 * the transpose of the last three columns of AP lower. 00092 * This covers the case N even and TRANSR = 'N'. 00093 * 00094 * RFP A RFP A 00095 * 00096 * 03 04 05 33 43 53 00097 * 13 14 15 00 44 54 00098 * 23 24 25 10 11 55 00099 * 33 34 35 20 21 22 00100 * 00 44 45 30 31 32 00101 * 01 11 55 40 41 42 00102 * 02 12 22 50 51 52 00103 * 00104 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00105 * transpose of RFP A above. One therefore gets: 00106 * 00107 * 00108 * RFP A RFP A 00109 * 00110 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00111 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00112 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00113 * 00114 * 00115 * We then consider Rectangular Full Packed (RFP) Format when N is 00116 * odd. We give an example where N = 5. 00117 * 00118 * AP is Upper AP is Lower 00119 * 00120 * 00 01 02 03 04 00 00121 * 11 12 13 14 10 11 00122 * 22 23 24 20 21 22 00123 * 33 34 30 31 32 33 00124 * 44 40 41 42 43 44 00125 * 00126 * 00127 * Let TRANSR = 'N'. RFP holds AP as follows: 00128 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00129 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00130 * the transpose of the first two columns of AP upper. 00131 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00132 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00133 * the transpose of the last two columns of AP lower. 00134 * This covers the case N odd and TRANSR = 'N'. 00135 * 00136 * RFP A RFP A 00137 * 00138 * 02 03 04 00 33 43 00139 * 12 13 14 10 11 44 00140 * 22 23 24 20 21 22 00141 * 00 33 34 30 31 32 00142 * 01 11 44 40 41 42 00143 * 00144 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00145 * transpose of RFP A above. One therefore gets: 00146 * 00147 * RFP A RFP A 00148 * 00149 * 02 12 22 00 01 00 10 20 30 40 50 00150 * 03 13 23 33 11 33 11 21 31 41 51 00151 * 04 14 24 34 44 43 44 22 32 42 52 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 DOUBLE PRECISION ONE 00157 PARAMETER ( ONE = 1.0D+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 LOGICAL LOWER, NISODD, NORMALTRANSR 00161 INTEGER N1, N2, K 00162 * .. 00163 * .. External Functions .. 00164 LOGICAL LSAME 00165 EXTERNAL LSAME 00166 * .. 00167 * .. External Subroutines .. 00168 EXTERNAL XERBLA, DSYRK, DPOTRF, DTRSM 00169 * .. 00170 * .. Intrinsic Functions .. 00171 INTRINSIC MOD 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Test the input parameters. 00176 * 00177 INFO = 0 00178 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00179 LOWER = LSAME( UPLO, 'L' ) 00180 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 00181 INFO = -1 00182 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00183 INFO = -2 00184 ELSE IF( N.LT.0 ) THEN 00185 INFO = -3 00186 END IF 00187 IF( INFO.NE.0 ) THEN 00188 CALL XERBLA( 'DPFTRF', -INFO ) 00189 RETURN 00190 END IF 00191 * 00192 * Quick return if possible 00193 * 00194 IF( N.EQ.0 ) 00195 $ RETURN 00196 * 00197 * If N is odd, set NISODD = .TRUE. 00198 * If N is even, set K = N/2 and NISODD = .FALSE. 00199 * 00200 IF( MOD( N, 2 ).EQ.0 ) THEN 00201 K = N / 2 00202 NISODD = .FALSE. 00203 ELSE 00204 NISODD = .TRUE. 00205 END IF 00206 * 00207 * Set N1 and N2 depending on LOWER 00208 * 00209 IF( LOWER ) THEN 00210 N2 = N / 2 00211 N1 = N - N2 00212 ELSE 00213 N1 = N / 2 00214 N2 = N - N1 00215 END IF 00216 * 00217 * start execution: there are eight cases 00218 * 00219 IF( NISODD ) THEN 00220 * 00221 * N is odd 00222 * 00223 IF( NORMALTRANSR ) THEN 00224 * 00225 * N is odd and TRANSR = 'N' 00226 * 00227 IF( LOWER ) THEN 00228 * 00229 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 00230 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 00231 * T1 -> a(0), T2 -> a(n), S -> a(n1) 00232 * 00233 CALL DPOTRF( 'L', N1, A( 0 ), N, INFO ) 00234 IF( INFO.GT.0 ) 00235 $ RETURN 00236 CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, A( 0 ), N, 00237 $ A( N1 ), N ) 00238 CALL DSYRK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE, 00239 $ A( N ), N ) 00240 CALL DPOTRF( 'U', N2, A( N ), N, INFO ) 00241 IF( INFO.GT.0 ) 00242 $ INFO = INFO + N1 00243 * 00244 ELSE 00245 * 00246 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 00247 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 00248 * T1 -> a(n2), T2 -> a(n1), S -> a(0) 00249 * 00250 CALL DPOTRF( 'L', N1, A( N2 ), N, INFO ) 00251 IF( INFO.GT.0 ) 00252 $ RETURN 00253 CALL DTRSM( 'L', 'L', 'N', 'N', N1, N2, ONE, A( N2 ), N, 00254 $ A( 0 ), N ) 00255 CALL DSYRK( 'U', 'T', N2, N1, -ONE, A( 0 ), N, ONE, 00256 $ A( N1 ), N ) 00257 CALL DPOTRF( 'U', N2, A( N1 ), N, INFO ) 00258 IF( INFO.GT.0 ) 00259 $ INFO = INFO + N1 00260 * 00261 END IF 00262 * 00263 ELSE 00264 * 00265 * N is odd and TRANSR = 'T' 00266 * 00267 IF( LOWER ) THEN 00268 * 00269 * SRPA for LOWER, TRANSPOSE and N is odd 00270 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 00271 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 00272 * 00273 CALL DPOTRF( 'U', N1, A( 0 ), N1, INFO ) 00274 IF( INFO.GT.0 ) 00275 $ RETURN 00276 CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, A( 0 ), N1, 00277 $ A( N1*N1 ), N1 ) 00278 CALL DSYRK( 'L', 'T', N2, N1, -ONE, A( N1*N1 ), N1, ONE, 00279 $ A( 1 ), N1 ) 00280 CALL DPOTRF( 'L', N2, A( 1 ), N1, INFO ) 00281 IF( INFO.GT.0 ) 00282 $ INFO = INFO + N1 00283 * 00284 ELSE 00285 * 00286 * SRPA for UPPER, TRANSPOSE and N is odd 00287 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 00288 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 00289 * 00290 CALL DPOTRF( 'U', N1, A( N2*N2 ), N2, INFO ) 00291 IF( INFO.GT.0 ) 00292 $ RETURN 00293 CALL DTRSM( 'R', 'U', 'N', 'N', N2, N1, ONE, A( N2*N2 ), 00294 $ N2, A( 0 ), N2 ) 00295 CALL DSYRK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE, 00296 $ A( N1*N2 ), N2 ) 00297 CALL DPOTRF( 'L', N2, A( N1*N2 ), N2, INFO ) 00298 IF( INFO.GT.0 ) 00299 $ INFO = INFO + N1 00300 * 00301 END IF 00302 * 00303 END IF 00304 * 00305 ELSE 00306 * 00307 * N is even 00308 * 00309 IF( NORMALTRANSR ) THEN 00310 * 00311 * N is even and TRANSR = 'N' 00312 * 00313 IF( LOWER ) THEN 00314 * 00315 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00316 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 00317 * T1 -> a(1), T2 -> a(0), S -> a(k+1) 00318 * 00319 CALL DPOTRF( 'L', K, A( 1 ), N+1, INFO ) 00320 IF( INFO.GT.0 ) 00321 $ RETURN 00322 CALL DTRSM( 'R', 'L', 'T', 'N', K, K, ONE, A( 1 ), N+1, 00323 $ A( K+1 ), N+1 ) 00324 CALL DSYRK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE, 00325 $ A( 0 ), N+1 ) 00326 CALL DPOTRF( 'U', K, A( 0 ), N+1, INFO ) 00327 IF( INFO.GT.0 ) 00328 $ INFO = INFO + K 00329 * 00330 ELSE 00331 * 00332 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00333 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 00334 * T1 -> a(k+1), T2 -> a(k), S -> a(0) 00335 * 00336 CALL DPOTRF( 'L', K, A( K+1 ), N+1, INFO ) 00337 IF( INFO.GT.0 ) 00338 $ RETURN 00339 CALL DTRSM( 'L', 'L', 'N', 'N', K, K, ONE, A( K+1 ), 00340 $ N+1, A( 0 ), N+1 ) 00341 CALL DSYRK( 'U', 'T', K, K, -ONE, A( 0 ), N+1, ONE, 00342 $ A( K ), N+1 ) 00343 CALL DPOTRF( 'U', K, A( K ), N+1, INFO ) 00344 IF( INFO.GT.0 ) 00345 $ INFO = INFO + K 00346 * 00347 END IF 00348 * 00349 ELSE 00350 * 00351 * N is even and TRANSR = 'T' 00352 * 00353 IF( LOWER ) THEN 00354 * 00355 * SRPA for LOWER, TRANSPOSE and N is even (see paper) 00356 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) 00357 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 00358 * 00359 CALL DPOTRF( 'U', K, A( 0+K ), K, INFO ) 00360 IF( INFO.GT.0 ) 00361 $ RETURN 00362 CALL DTRSM( 'L', 'U', 'T', 'N', K, K, ONE, A( K ), N1, 00363 $ A( K*( K+1 ) ), K ) 00364 CALL DSYRK( 'L', 'T', K, K, -ONE, A( K*( K+1 ) ), K, ONE, 00365 $ A( 0 ), K ) 00366 CALL DPOTRF( 'L', K, A( 0 ), K, INFO ) 00367 IF( INFO.GT.0 ) 00368 $ INFO = INFO + K 00369 * 00370 ELSE 00371 * 00372 * SRPA for UPPER, TRANSPOSE and N is even (see paper) 00373 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0) 00374 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 00375 * 00376 CALL DPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO ) 00377 IF( INFO.GT.0 ) 00378 $ RETURN 00379 CALL DTRSM( 'R', 'U', 'N', 'N', K, K, ONE, 00380 $ A( K*( K+1 ) ), K, A( 0 ), K ) 00381 CALL DSYRK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE, 00382 $ A( K*K ), K ) 00383 CALL DPOTRF( 'L', K, A( K*K ), K, INFO ) 00384 IF( INFO.GT.0 ) 00385 $ INFO = INFO + K 00386 * 00387 END IF 00388 * 00389 END IF 00390 * 00391 END IF 00392 * 00393 RETURN 00394 * 00395 * End of DPFTRF 00396 * 00397 END
1.7.3 
