![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPFTRF( TRANSR, UPLO, N, A, INFO ) 2: * 3: * -- LAPACK routine (version 3.3.0) -- 4: * 5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 6: * November 2010 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