![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHFRK( 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: COMPLEX*16 A( LDA, * ), C( * ) 20: * .. 21: * 22: * Purpose 23: * ======= 24: * 25: * Level 3 BLAS like routine for C in RFP Format. 26: * 27: * ZHFRK performs one of the Hermitian rank--k operations 28: * 29: * C := alpha*A*conjg( A' ) + beta*C, 30: * 31: * or 32: * 33: * C := alpha*conjg( A' )*A + beta*C, 34: * 35: * where alpha and beta are real scalars, C is an n--by--n Hermitian 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: * = 'C': The Conjugate-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*conjg( A' ) + beta*C. 64: * 65: * TRANS = 'C' or 'c' C := alpha*conjg( 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 77: * TRANS = 'C' or 'c', K specifies the number of rows of the 78: * matrix A. K 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) COMPLEX*16 array of 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: * C (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 105: * On entry, the matrix A in RFP Format. RFP Format is 106: * described by TRANSR, UPLO and N. Note that the imaginary 107: * parts of the diagonal elements need not be set, they are 108: * assumed to be zero, and on exit they are set to zero. 109: * 110: * Arguments 111: * ========== 112: * 113: * .. 114: * .. Parameters .. 115: DOUBLE PRECISION ONE, ZERO 116: COMPLEX*16 CZERO 117: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 118: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 119: * .. 120: * .. Local Scalars .. 121: LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS 122: INTEGER INFO, NROWA, J, NK, N1, N2 123: COMPLEX*16 CALPHA, CBETA 124: * .. 125: * .. External Functions .. 126: LOGICAL LSAME 127: EXTERNAL LSAME 128: * .. 129: * .. External Subroutines .. 130: EXTERNAL XERBLA, ZGEMM, ZHERK 131: * .. 132: * .. Intrinsic Functions .. 133: INTRINSIC MAX, DCMPLX 134: * .. 135: * .. Executable Statements .. 136: * 137: * 138: * Test the input parameters. 139: * 140: INFO = 0 141: NORMALTRANSR = LSAME( TRANSR, 'N' ) 142: LOWER = LSAME( UPLO, 'L' ) 143: NOTRANS = LSAME( TRANS, 'N' ) 144: * 145: IF( NOTRANS ) THEN 146: NROWA = N 147: ELSE 148: NROWA = K 149: END IF 150: * 151: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 152: INFO = -1 153: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 154: INFO = -2 155: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 156: INFO = -3 157: ELSE IF( N.LT.0 ) THEN 158: INFO = -4 159: ELSE IF( K.LT.0 ) THEN 160: INFO = -5 161: ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN 162: INFO = -8 163: END IF 164: IF( INFO.NE.0 ) THEN 165: CALL XERBLA( 'ZHFRK ', -INFO ) 166: RETURN 167: END IF 168: * 169: * Quick return if possible. 170: * 171: * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not 172: * done (it is in ZHERK for example) and left in the general case. 173: * 174: IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. 175: + ( BETA.EQ.ONE ) ) )RETURN 176: * 177: IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN 178: DO J = 1, ( ( N*( N+1 ) ) / 2 ) 179: C( J ) = CZERO 180: END DO 181: RETURN 182: END IF 183: * 184: CALPHA = DCMPLX( ALPHA, ZERO ) 185: CBETA = DCMPLX( BETA, ZERO ) 186: * 187: * C is N-by-N. 188: * If N is odd, set NISODD = .TRUE., and N1 and N2. 189: * If N is even, NISODD = .FALSE., and NK. 190: * 191: IF( MOD( N, 2 ).EQ.0 ) THEN 192: NISODD = .FALSE. 193: NK = N / 2 194: ELSE 195: NISODD = .TRUE. 196: IF( LOWER ) THEN 197: N2 = N / 2 198: N1 = N - N2 199: ELSE 200: N1 = N / 2 201: N2 = N - N1 202: END IF 203: END IF 204: * 205: IF( NISODD ) THEN 206: * 207: * N is odd 208: * 209: IF( NORMALTRANSR ) THEN 210: * 211: * N is odd and TRANSR = 'N' 212: * 213: IF( LOWER ) THEN 214: * 215: * N is odd, TRANSR = 'N', and UPLO = 'L' 216: * 217: IF( NOTRANS ) THEN 218: * 219: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 220: * 221: CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 222: + BETA, C( 1 ), N ) 223: CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 224: + BETA, C( N+1 ), N ) 225: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 226: + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 227: * 228: ELSE 229: * 230: * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 231: * 232: CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 233: + BETA, C( 1 ), N ) 234: CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 235: + BETA, C( N+1 ), N ) 236: CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 237: + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 238: * 239: END IF 240: * 241: ELSE 242: * 243: * N is odd, TRANSR = 'N', and UPLO = 'U' 244: * 245: IF( NOTRANS ) THEN 246: * 247: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 248: * 249: CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 250: + BETA, C( N2+1 ), N ) 251: CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA, 252: + BETA, C( N1+1 ), N ) 253: CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 254: + LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N ) 255: * 256: ELSE 257: * 258: * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 259: * 260: CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 261: + BETA, C( N2+1 ), N ) 262: CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N2 ), LDA, 263: + BETA, C( N1+1 ), N ) 264: CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 265: + LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N ) 266: * 267: END IF 268: * 269: END IF 270: * 271: ELSE 272: * 273: * N is odd, and TRANSR = 'C' 274: * 275: IF( LOWER ) THEN 276: * 277: * N is odd, TRANSR = 'C', and UPLO = 'L' 278: * 279: IF( NOTRANS ) THEN 280: * 281: * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 282: * 283: CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 284: + BETA, C( 1 ), N1 ) 285: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 286: + BETA, C( 2 ), N1 ) 287: CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 288: + LDA, A( N1+1, 1 ), LDA, CBETA, 289: + C( N1*N1+1 ), N1 ) 290: * 291: ELSE 292: * 293: * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 294: * 295: CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 296: + BETA, C( 1 ), N1 ) 297: CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 298: + BETA, C( 2 ), N1 ) 299: CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 300: + LDA, A( 1, N1+1 ), LDA, CBETA, 301: + C( N1*N1+1 ), N1 ) 302: * 303: END IF 304: * 305: ELSE 306: * 307: * N is odd, TRANSR = 'C', and UPLO = 'U' 308: * 309: IF( NOTRANS ) THEN 310: * 311: * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 312: * 313: CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 314: + BETA, C( N2*N2+1 ), N2 ) 315: CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 316: + BETA, C( N1*N2+1 ), N2 ) 317: CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 318: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 319: * 320: ELSE 321: * 322: * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 323: * 324: CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 325: + BETA, C( N2*N2+1 ), N2 ) 326: CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 327: + BETA, C( N1*N2+1 ), N2 ) 328: CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 329: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 330: * 331: END IF 332: * 333: END IF 334: * 335: END IF 336: * 337: ELSE 338: * 339: * N is even 340: * 341: IF( NORMALTRANSR ) THEN 342: * 343: * N is even and TRANSR = 'N' 344: * 345: IF( LOWER ) THEN 346: * 347: * N is even, TRANSR = 'N', and UPLO = 'L' 348: * 349: IF( NOTRANS ) THEN 350: * 351: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 352: * 353: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 354: + BETA, C( 2 ), N+1 ) 355: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 356: + BETA, C( 1 ), N+1 ) 357: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 358: + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 359: + N+1 ) 360: * 361: ELSE 362: * 363: * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 364: * 365: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 366: + BETA, C( 2 ), N+1 ) 367: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 368: + BETA, C( 1 ), N+1 ) 369: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 370: + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 371: + N+1 ) 372: * 373: END IF 374: * 375: ELSE 376: * 377: * N is even, TRANSR = 'N', and UPLO = 'U' 378: * 379: IF( NOTRANS ) THEN 380: * 381: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 382: * 383: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 384: + BETA, C( NK+2 ), N+1 ) 385: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 386: + BETA, C( NK+1 ), N+1 ) 387: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 388: + LDA, A( NK+1, 1 ), LDA, CBETA, C( 1 ), 389: + N+1 ) 390: * 391: ELSE 392: * 393: * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 394: * 395: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 396: + BETA, C( NK+2 ), N+1 ) 397: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 398: + BETA, C( NK+1 ), N+1 ) 399: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 400: + LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ), 401: + N+1 ) 402: * 403: END IF 404: * 405: END IF 406: * 407: ELSE 408: * 409: * N is even, and TRANSR = 'C' 410: * 411: IF( LOWER ) THEN 412: * 413: * N is even, TRANSR = 'C', and UPLO = 'L' 414: * 415: IF( NOTRANS ) THEN 416: * 417: * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 418: * 419: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 420: + BETA, C( NK+1 ), NK ) 421: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 422: + BETA, C( 1 ), NK ) 423: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 424: + LDA, A( NK+1, 1 ), LDA, CBETA, 425: + C( ( ( NK+1 )*NK )+1 ), NK ) 426: * 427: ELSE 428: * 429: * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 430: * 431: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 432: + BETA, C( NK+1 ), NK ) 433: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 434: + BETA, C( 1 ), NK ) 435: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 436: + LDA, A( 1, NK+1 ), LDA, CBETA, 437: + C( ( ( NK+1 )*NK )+1 ), NK ) 438: * 439: END IF 440: * 441: ELSE 442: * 443: * N is even, TRANSR = 'C', and UPLO = 'U' 444: * 445: IF( NOTRANS ) THEN 446: * 447: * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 448: * 449: CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 450: + BETA, C( NK*( NK+1 )+1 ), NK ) 451: CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 452: + BETA, C( NK*NK+1 ), NK ) 453: CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 454: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 455: * 456: ELSE 457: * 458: * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 459: * 460: CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 461: + BETA, C( NK*( NK+1 )+1 ), NK ) 462: CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 463: + BETA, C( NK*NK+1 ), NK ) 464: CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 465: + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 466: * 467: END IF 468: * 469: END IF 470: * 471: END IF 472: * 473: END IF 474: * 475: RETURN 476: * 477: * End of ZHFRK 478: * 479: END