![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, 2: + C ) 3: * 4: * -- LAPACK routine (version 3.3.0) -- 5: * 6: * -- Contributed by Julien Langou of the Univ. of Colorado Denver -- 7: * November 2010 8: * 9: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 11: * 12: * .. 13: * .. Scalar Arguments .. 14: DOUBLE PRECISION ALPHA, BETA 15: INTEGER K, LDA, N 16: CHARACTER TRANS, TRANSR, UPLO 17: * .. 18: * .. Array Arguments .. 19: DOUBLE PRECISION A( LDA, * ), C( * ) 20: * .. 21: * 22: * Purpose 23: * ======= 24: * 25: * Level 3 BLAS like routine for C in RFP Format. 26: * 27: * DSFRK performs one of the symmetric rank--k operations 28: * 29: * C := alpha*A*A' + beta*C, 30: * 31: * or 32: * 33: * C := alpha*A'*A + beta*C, 34: * 35: * where alpha and beta are real scalars, C is an n--by--n symmetric 36: * matrix and A is an n--by--k matrix in the first case and a k--by--n 37: * matrix in the second case. 38: * 39: * Arguments 40: * ========== 41: * 42: * TRANSR (input) CHARACTER*1 43: * = 'N': The Normal Form of RFP A is stored; 44: * = 'T': The Transpose Form of RFP A is stored. 45: * 46: * UPLO (input) CHARACTER*1 47: * On entry, UPLO specifies whether the upper or lower 48: * triangular part of the array C is to be referenced as 49: * follows: 50: * 51: * UPLO = 'U' or 'u' Only the upper triangular part of C 52: * is to be referenced. 53: * 54: * UPLO = 'L' or 'l' Only the lower triangular part of C 55: * is to be referenced. 56: * 57: * Unchanged on exit. 58: * 59: * TRANS (input) CHARACTER*1 60: * On entry, TRANS specifies the operation to be performed as 61: * follows: 62: * 63: * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. 64: * 65: * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. 66: * 67: * Unchanged on exit. 68: * 69: * N (input) INTEGER 70: * On entry, N specifies the order of the matrix C. N must be 71: * at least zero. 72: * Unchanged on exit. 73: * 74: * K (input) INTEGER 75: * On entry with TRANS = 'N' or 'n', K specifies the number 76: * of columns of the matrix A, and on entry with TRANS = 'T' 77: * or 't', K specifies the number of rows of the matrix A. K 78: * must be at least zero. 79: * Unchanged on exit. 80: * 81: * ALPHA (input) DOUBLE PRECISION 82: * On entry, ALPHA specifies the scalar alpha. 83: * Unchanged on exit. 84: * 85: * A (input) DOUBLE PRECISION array, dimension (LDA,ka) 86: * where KA 87: * is K when TRANS = 'N' or 'n', and is N otherwise. Before 88: * entry with TRANS = 'N' or 'n', the leading N--by--K part of 89: * the array A must contain the matrix A, otherwise the leading 90: * K--by--N part of the array A must contain the matrix A. 91: * Unchanged on exit. 92: * 93: * LDA (input) INTEGER 94: * On entry, LDA specifies the first dimension of A as declared 95: * in the calling (sub) program. When TRANS = 'N' or 'n' 96: * then LDA must be at least max( 1, n ), otherwise LDA must 97: * be at least max( 1, k ). 98: * Unchanged on exit. 99: * 100: * BETA (input) DOUBLE PRECISION 101: * On entry, BETA specifies the scalar beta. 102: * Unchanged on exit. 103: * 104: * 105: * C (input/output) DOUBLE PRECISION array, dimension (NT) 106: * NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP 107: * Format. RFP Format is described by TRANSR, UPLO and N. 108: * 109: * Arguments 110: * ========== 111: * 112: * .. 113: * .. Parameters .. 114: DOUBLE PRECISION ONE, ZERO 115: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 116: * .. 117: * .. Local Scalars .. 118: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS 119: INTEGER INFO, NROWA, J, NK, N1, N2 120: * .. 121: * .. External Functions .. 122: LOGICAL LSAME 123: EXTERNAL LSAME 124: * .. 125: * .. External Subroutines .. 126: EXTERNAL XERBLA, DGEMM, DSYRK 127: * .. 128: * .. Intrinsic Functions .. 129: INTRINSIC MAX 130: * .. 131: * .. Executable Statements .. 132: * 133: * Test the input parameters. 134: * 135: INFO = 0 136: NORMALTRANSR = LSAME( TRANSR, 'N' ) 137: LOWER = LSAME( UPLO, 'L' ) 138: NOTRANS = LSAME( TRANS, 'N' ) 139: * 140: IF( NOTRANS ) THEN 141: NROWA = N 142: ELSE 143: NROWA = K 144: END IF 145: * 146: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 147: INFO = -1 148: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 149: INFO = -2 150: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 151: INFO = -3 152: ELSE IF( N.LT.0 ) THEN 153: INFO = -4 154: ELSE IF( K.LT.0 ) THEN 155: INFO = -5 156: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN 157: INFO = -8 158: END IF 159: IF( INFO.NE.0 ) THEN 160: CALL XERBLA( 'DSFRK ', -INFO ) 161: RETURN 162: END IF 163: * 164: * Quick return if possible. 165: * 166: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not 167: * done (it is in DSYRK for example) and left in the general case. 168: * 169: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. 170: + ( BETA.EQ.ONE ) ) )RETURN 171: * 172: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN 173: DO J = 1, ( ( N*( N+1 ) ) / 2 ) 174: C( J ) = ZERO 175: END DO 176: RETURN 177: END IF 178: * 179: * C is N-by-N. 180: * If N is odd, set NISODD = .TRUE., and N1 and N2. 181: * If N is even, NISODD = .FALSE., and NK. 182: * 183: IF( MOD( N, 2 ).EQ.0 ) THEN 184: NISODD = .FALSE. 185: NK = N / 2 186: ELSE 187: NISODD = .TRUE. 188: IF( LOWER ) THEN 189: N2 = N / 2 190: N1 = N - N2 191: ELSE 192: N1 = N / 2 193: N2 = N - N1 194: END IF 195: END IF 196: * 197: IF( NISODD ) THEN 198: * 199: * N is odd 200: * 201: IF( NORMALTRANSR ) THEN 202: * 203: * N is odd and TRANSR = 'N' 204: * 205: IF( LOWER ) THEN 206: * 207: * N is odd, TRANSR = 'N', and UPLO = 'L' 208: * 209: IF( NOTRANS ) THEN 210: * 211: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 212: * 213: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 214: + BETA, C( 1 ), N ) 215: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 216: + BETA, C( N+1 ), N ) 217: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ), 218: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N ) 219: * 220: ELSE 221: * 222: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' 223: * 224: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 225: + BETA, C( 1 ), N ) 226: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 227: + BETA, C( N+1 ), N ) 228: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ), 229: + LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N ) 230: * 231: END IF 232: * 233: ELSE 234: * 235: * N is odd, TRANSR = 'N', and UPLO = 'U' 236: * 237: IF( NOTRANS ) THEN 238: * 239: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 240: * 241: CALL DSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 242: + BETA, C( N2+1 ), N ) 243: CALL DSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA, 244: + BETA, C( N1+1 ), N ) 245: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ), 246: + LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N ) 247: * 248: ELSE 249: * 250: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' 251: * 252: CALL DSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 253: + BETA, C( N2+1 ), N ) 254: CALL DSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA, 255: + BETA, C( N1+1 ), N ) 256: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ), 257: + LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N ) 258: * 259: END IF 260: * 261: END IF 262: * 263: ELSE 264: * 265: * N is odd, and TRANSR = 'T' 266: * 267: IF( LOWER ) THEN 268: * 269: * N is odd, TRANSR = 'T', and UPLO = 'L' 270: * 271: IF( NOTRANS ) THEN 272: * 273: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' 274: * 275: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 276: + BETA, C( 1 ), N1 ) 277: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 278: + BETA, C( 2 ), N1 ) 279: CALL DGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ), 280: + LDA, A( N1+1, 1 ), LDA, BETA, 281: + C( N1*N1+1 ), N1 ) 282: * 283: ELSE 284: * 285: * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' 286: * 287: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 288: + BETA, C( 1 ), N1 ) 289: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 290: + BETA, C( 2 ), N1 ) 291: CALL DGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ), 292: + LDA, A( 1, N1+1 ), LDA, BETA, 293: + C( N1*N1+1 ), N1 ) 294: * 295: END IF 296: * 297: ELSE 298: * 299: * N is odd, TRANSR = 'T', and UPLO = 'U' 300: * 301: IF( NOTRANS ) THEN 302: * 303: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' 304: * 305: CALL DSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 306: + BETA, C( N2*N2+1 ), N2 ) 307: CALL DSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 308: + BETA, C( N1*N2+1 ), N2 ) 309: CALL DGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ), 310: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 ) 311: * 312: ELSE 313: * 314: * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' 315: * 316: CALL DSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 317: + BETA, C( N2*N2+1 ), N2 ) 318: CALL DSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 319: + BETA, C( N1*N2+1 ), N2 ) 320: CALL DGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ), 321: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 ) 322: * 323: END IF 324: * 325: END IF 326: * 327: END IF 328: * 329: ELSE 330: * 331: * N is even 332: * 333: IF( NORMALTRANSR ) THEN 334: * 335: * N is even and TRANSR = 'N' 336: * 337: IF( LOWER ) THEN 338: * 339: * N is even, TRANSR = 'N', and UPLO = 'L' 340: * 341: IF( NOTRANS ) THEN 342: * 343: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 344: * 345: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 346: + BETA, C( 2 ), N+1 ) 347: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 348: + BETA, C( 1 ), N+1 ) 349: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ), 350: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ), 351: + N+1 ) 352: * 353: ELSE 354: * 355: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' 356: * 357: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 358: + BETA, C( 2 ), N+1 ) 359: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 360: + BETA, C( 1 ), N+1 ) 361: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ), 362: + LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ), 363: + N+1 ) 364: * 365: END IF 366: * 367: ELSE 368: * 369: * N is even, TRANSR = 'N', and UPLO = 'U' 370: * 371: IF( NOTRANS ) THEN 372: * 373: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 374: * 375: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 376: + BETA, C( NK+2 ), N+1 ) 377: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 378: + BETA, C( NK+1 ), N+1 ) 379: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ), 380: + LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ), 381: + N+1 ) 382: * 383: ELSE 384: * 385: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' 386: * 387: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 388: + BETA, C( NK+2 ), N+1 ) 389: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 390: + BETA, C( NK+1 ), N+1 ) 391: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ), 392: + LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ), 393: + N+1 ) 394: * 395: END IF 396: * 397: END IF 398: * 399: ELSE 400: * 401: * N is even, and TRANSR = 'T' 402: * 403: IF( LOWER ) THEN 404: * 405: * N is even, TRANSR = 'T', and UPLO = 'L' 406: * 407: IF( NOTRANS ) THEN 408: * 409: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' 410: * 411: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 412: + BETA, C( NK+1 ), NK ) 413: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 414: + BETA, C( 1 ), NK ) 415: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ), 416: + LDA, A( NK+1, 1 ), LDA, BETA, 417: + C( ( ( NK+1 )*NK )+1 ), NK ) 418: * 419: ELSE 420: * 421: * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' 422: * 423: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 424: + BETA, C( NK+1 ), NK ) 425: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 426: + BETA, C( 1 ), NK ) 427: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ), 428: + LDA, A( 1, NK+1 ), LDA, BETA, 429: + C( ( ( NK+1 )*NK )+1 ), NK ) 430: * 431: END IF 432: * 433: ELSE 434: * 435: * N is even, TRANSR = 'T', and UPLO = 'U' 436: * 437: IF( NOTRANS ) THEN 438: * 439: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' 440: * 441: CALL DSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 442: + BETA, C( NK*( NK+1 )+1 ), NK ) 443: CALL DSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 444: + BETA, C( NK*NK+1 ), NK ) 445: CALL DGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ), 446: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK ) 447: * 448: ELSE 449: * 450: * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' 451: * 452: CALL DSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 453: + BETA, C( NK*( NK+1 )+1 ), NK ) 454: CALL DSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 455: + BETA, C( NK*NK+1 ), NK ) 456: CALL DGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ), 457: + LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK ) 458: * 459: END IF 460: * 461: END IF 462: * 463: END IF 464: * 465: END IF 466: * 467: RETURN 468: * 469: * End of DSFRK 470: * 471: END