![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, 2: $ T, LDT, C, LDC, WORK, LDWORK ) 3: IMPLICIT NONE 4: * 5: * -- LAPACK auxiliary routine (version 3.2) -- 6: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8: * November 2006 9: * 10: * .. Scalar Arguments .. 11: CHARACTER DIRECT, SIDE, STOREV, TRANS 12: INTEGER K, LDC, LDT, LDV, LDWORK, M, N 13: * .. 14: * .. Array Arguments .. 15: DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), 16: $ WORK( LDWORK, * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * DLARFB applies a real block reflector H or its transpose H' to a 23: * real m by n matrix C, from either the left or the right. 24: * 25: * Arguments 26: * ========= 27: * 28: * SIDE (input) CHARACTER*1 29: * = 'L': apply H or H' from the Left 30: * = 'R': apply H or H' from the Right 31: * 32: * TRANS (input) CHARACTER*1 33: * = 'N': apply H (No transpose) 34: * = 'T': apply H' (Transpose) 35: * 36: * DIRECT (input) CHARACTER*1 37: * Indicates how H is formed from a product of elementary 38: * reflectors 39: * = 'F': H = H(1) H(2) . . . H(k) (Forward) 40: * = 'B': H = H(k) . . . H(2) H(1) (Backward) 41: * 42: * STOREV (input) CHARACTER*1 43: * Indicates how the vectors which define the elementary 44: * reflectors are stored: 45: * = 'C': Columnwise 46: * = 'R': Rowwise 47: * 48: * M (input) INTEGER 49: * The number of rows of the matrix C. 50: * 51: * N (input) INTEGER 52: * The number of columns of the matrix C. 53: * 54: * K (input) INTEGER 55: * The order of the matrix T (= the number of elementary 56: * reflectors whose product defines the block reflector). 57: * 58: * V (input) DOUBLE PRECISION array, dimension 59: * (LDV,K) if STOREV = 'C' 60: * (LDV,M) if STOREV = 'R' and SIDE = 'L' 61: * (LDV,N) if STOREV = 'R' and SIDE = 'R' 62: * The matrix V. See further details. 63: * 64: * LDV (input) INTEGER 65: * The leading dimension of the array V. 66: * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); 67: * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); 68: * if STOREV = 'R', LDV >= K. 69: * 70: * T (input) DOUBLE PRECISION array, dimension (LDT,K) 71: * The triangular k by k matrix T in the representation of the 72: * block reflector. 73: * 74: * LDT (input) INTEGER 75: * The leading dimension of the array T. LDT >= K. 76: * 77: * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 78: * On entry, the m by n matrix C. 79: * On exit, C is overwritten by H*C or H'*C or C*H or C*H'. 80: * 81: * LDC (input) INTEGER 82: * The leading dimension of the array C. LDA >= max(1,M). 83: * 84: * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) 85: * 86: * LDWORK (input) INTEGER 87: * The leading dimension of the array WORK. 88: * If SIDE = 'L', LDWORK >= max(1,N); 89: * if SIDE = 'R', LDWORK >= max(1,M). 90: * 91: * ===================================================================== 92: * 93: * .. Parameters .. 94: DOUBLE PRECISION ONE 95: PARAMETER ( ONE = 1.0D+0 ) 96: * .. 97: * .. Local Scalars .. 98: CHARACTER TRANST 99: INTEGER I, J, LASTV, LASTC 100: * .. 101: * .. External Functions .. 102: LOGICAL LSAME 103: INTEGER ILADLR, ILADLC 104: EXTERNAL LSAME, ILADLR, ILADLC 105: * .. 106: * .. External Subroutines .. 107: EXTERNAL DCOPY, DGEMM, DTRMM 108: * .. 109: * .. Executable Statements .. 110: * 111: * Quick return if possible 112: * 113: IF( M.LE.0 .OR. N.LE.0 ) 114: $ RETURN 115: * 116: IF( LSAME( TRANS, 'N' ) ) THEN 117: TRANST = 'T' 118: ELSE 119: TRANST = 'N' 120: END IF 121: * 122: IF( LSAME( STOREV, 'C' ) ) THEN 123: * 124: IF( LSAME( DIRECT, 'F' ) ) THEN 125: * 126: * Let V = ( V1 ) (first K rows) 127: * ( V2 ) 128: * where V1 is unit lower triangular. 129: * 130: IF( LSAME( SIDE, 'L' ) ) THEN 131: * 132: * Form H * C or H' * C where C = ( C1 ) 133: * ( C2 ) 134: * 135: LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 136: LASTC = ILADLC( LASTV, N, C, LDC ) 137: * 138: * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 139: * 140: * W := C1' 141: * 142: DO 10 J = 1, K 143: CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 144: 10 CONTINUE 145: * 146: * W := W * V1 147: * 148: CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 149: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 150: IF( LASTV.GT.K ) THEN 151: * 152: * W := W + C2'*V2 153: * 154: CALL DGEMM( 'Transpose', 'No transpose', 155: $ LASTC, K, LASTV-K, 156: $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, 157: $ ONE, WORK, LDWORK ) 158: END IF 159: * 160: * W := W * T' or W * T 161: * 162: CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 163: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 164: * 165: * C := C - V * W' 166: * 167: IF( LASTV.GT.K ) THEN 168: * 169: * C2 := C2 - V2 * W' 170: * 171: CALL DGEMM( 'No transpose', 'Transpose', 172: $ LASTV-K, LASTC, K, 173: $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, 174: $ C( K+1, 1 ), LDC ) 175: END IF 176: * 177: * W := W * V1' 178: * 179: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 180: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 181: * 182: * C1 := C1 - W' 183: * 184: DO 30 J = 1, K 185: DO 20 I = 1, LASTC 186: C( J, I ) = C( J, I ) - WORK( I, J ) 187: 20 CONTINUE 188: 30 CONTINUE 189: * 190: ELSE IF( LSAME( SIDE, 'R' ) ) THEN 191: * 192: * Form C * H or C * H' where C = ( C1 C2 ) 193: * 194: LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 195: LASTC = ILADLR( M, LASTV, C, LDC ) 196: * 197: * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 198: * 199: * W := C1 200: * 201: DO 40 J = 1, K 202: CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 203: 40 CONTINUE 204: * 205: * W := W * V1 206: * 207: CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 208: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 209: IF( LASTV.GT.K ) THEN 210: * 211: * W := W + C2 * V2 212: * 213: CALL DGEMM( 'No transpose', 'No transpose', 214: $ LASTC, K, LASTV-K, 215: $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 216: $ ONE, WORK, LDWORK ) 217: END IF 218: * 219: * W := W * T or W * T' 220: * 221: CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 222: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 223: * 224: * C := C - W * V' 225: * 226: IF( LASTV.GT.K ) THEN 227: * 228: * C2 := C2 - W * V2' 229: * 230: CALL DGEMM( 'No transpose', 'Transpose', 231: $ LASTC, LASTV-K, K, 232: $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, 233: $ C( 1, K+1 ), LDC ) 234: END IF 235: * 236: * W := W * V1' 237: * 238: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 239: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 240: * 241: * C1 := C1 - W 242: * 243: DO 60 J = 1, K 244: DO 50 I = 1, LASTC 245: C( I, J ) = C( I, J ) - WORK( I, J ) 246: 50 CONTINUE 247: 60 CONTINUE 248: END IF 249: * 250: ELSE 251: * 252: * Let V = ( V1 ) 253: * ( V2 ) (last K rows) 254: * where V2 is unit upper triangular. 255: * 256: IF( LSAME( SIDE, 'L' ) ) THEN 257: * 258: * Form H * C or H' * C where C = ( C1 ) 259: * ( C2 ) 260: * 261: LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 262: LASTC = ILADLC( LASTV, N, C, LDC ) 263: * 264: * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 265: * 266: * W := C2' 267: * 268: DO 70 J = 1, K 269: CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 270: $ WORK( 1, J ), 1 ) 271: 70 CONTINUE 272: * 273: * W := W * V2 274: * 275: CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 276: $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 277: $ WORK, LDWORK ) 278: IF( LASTV.GT.K ) THEN 279: * 280: * W := W + C1'*V1 281: * 282: CALL DGEMM( 'Transpose', 'No transpose', 283: $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 284: $ ONE, WORK, LDWORK ) 285: END IF 286: * 287: * W := W * T' or W * T 288: * 289: CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 290: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 291: * 292: * C := C - V * W' 293: * 294: IF( LASTV.GT.K ) THEN 295: * 296: * C1 := C1 - V1 * W' 297: * 298: CALL DGEMM( 'No transpose', 'Transpose', 299: $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 300: $ ONE, C, LDC ) 301: END IF 302: * 303: * W := W * V2' 304: * 305: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 306: $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 307: $ WORK, LDWORK ) 308: * 309: * C2 := C2 - W' 310: * 311: DO 90 J = 1, K 312: DO 80 I = 1, LASTC 313: C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 314: 80 CONTINUE 315: 90 CONTINUE 316: * 317: ELSE IF( LSAME( SIDE, 'R' ) ) THEN 318: * 319: * Form C * H or C * H' where C = ( C1 C2 ) 320: * 321: LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 322: LASTC = ILADLR( M, LASTV, C, LDC ) 323: * 324: * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 325: * 326: * W := C2 327: * 328: DO 100 J = 1, K 329: CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) 330: 100 CONTINUE 331: * 332: * W := W * V2 333: * 334: CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 335: $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 336: $ WORK, LDWORK ) 337: IF( LASTV.GT.K ) THEN 338: * 339: * W := W + C1 * V1 340: * 341: CALL DGEMM( 'No transpose', 'No transpose', 342: $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 343: $ ONE, WORK, LDWORK ) 344: END IF 345: * 346: * W := W * T or W * T' 347: * 348: CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 349: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 350: * 351: * C := C - W * V' 352: * 353: IF( LASTV.GT.K ) THEN 354: * 355: * C1 := C1 - W * V1' 356: * 357: CALL DGEMM( 'No transpose', 'Transpose', 358: $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 359: $ ONE, C, LDC ) 360: END IF 361: * 362: * W := W * V2' 363: * 364: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 365: $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 366: $ WORK, LDWORK ) 367: * 368: * C2 := C2 - W 369: * 370: DO 120 J = 1, K 371: DO 110 I = 1, LASTC 372: C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 373: 110 CONTINUE 374: 120 CONTINUE 375: END IF 376: END IF 377: * 378: ELSE IF( LSAME( STOREV, 'R' ) ) THEN 379: * 380: IF( LSAME( DIRECT, 'F' ) ) THEN 381: * 382: * Let V = ( V1 V2 ) (V1: first K columns) 383: * where V1 is unit upper triangular. 384: * 385: IF( LSAME( SIDE, 'L' ) ) THEN 386: * 387: * Form H * C or H' * C where C = ( C1 ) 388: * ( C2 ) 389: * 390: LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 391: LASTC = ILADLC( LASTV, N, C, LDC ) 392: * 393: * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 394: * 395: * W := C1' 396: * 397: DO 130 J = 1, K 398: CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 399: 130 CONTINUE 400: * 401: * W := W * V1' 402: * 403: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 404: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 405: IF( LASTV.GT.K ) THEN 406: * 407: * W := W + C2'*V2' 408: * 409: CALL DGEMM( 'Transpose', 'Transpose', 410: $ LASTC, K, LASTV-K, 411: $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, 412: $ ONE, WORK, LDWORK ) 413: END IF 414: * 415: * W := W * T' or W * T 416: * 417: CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 418: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 419: * 420: * C := C - V' * W' 421: * 422: IF( LASTV.GT.K ) THEN 423: * 424: * C2 := C2 - V2' * W' 425: * 426: CALL DGEMM( 'Transpose', 'Transpose', 427: $ LASTV-K, LASTC, K, 428: $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, 429: $ ONE, C( K+1, 1 ), LDC ) 430: END IF 431: * 432: * W := W * V1 433: * 434: CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 435: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 436: * 437: * C1 := C1 - W' 438: * 439: DO 150 J = 1, K 440: DO 140 I = 1, LASTC 441: C( J, I ) = C( J, I ) - WORK( I, J ) 442: 140 CONTINUE 443: 150 CONTINUE 444: * 445: ELSE IF( LSAME( SIDE, 'R' ) ) THEN 446: * 447: * Form C * H or C * H' where C = ( C1 C2 ) 448: * 449: LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 450: LASTC = ILADLR( M, LASTV, C, LDC ) 451: * 452: * W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 453: * 454: * W := C1 455: * 456: DO 160 J = 1, K 457: CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 458: 160 CONTINUE 459: * 460: * W := W * V1' 461: * 462: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 463: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 464: IF( LASTV.GT.K ) THEN 465: * 466: * W := W + C2 * V2' 467: * 468: CALL DGEMM( 'No transpose', 'Transpose', 469: $ LASTC, K, LASTV-K, 470: $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, 471: $ ONE, WORK, LDWORK ) 472: END IF 473: * 474: * W := W * T or W * T' 475: * 476: CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 477: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 478: * 479: * C := C - W * V 480: * 481: IF( LASTV.GT.K ) THEN 482: * 483: * C2 := C2 - W * V2 484: * 485: CALL DGEMM( 'No transpose', 'No transpose', 486: $ LASTC, LASTV-K, K, 487: $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, 488: $ ONE, C( 1, K+1 ), LDC ) 489: END IF 490: * 491: * W := W * V1 492: * 493: CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 494: $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 495: * 496: * C1 := C1 - W 497: * 498: DO 180 J = 1, K 499: DO 170 I = 1, LASTC 500: C( I, J ) = C( I, J ) - WORK( I, J ) 501: 170 CONTINUE 502: 180 CONTINUE 503: * 504: END IF 505: * 506: ELSE 507: * 508: * Let V = ( V1 V2 ) (V2: last K columns) 509: * where V2 is unit lower triangular. 510: * 511: IF( LSAME( SIDE, 'L' ) ) THEN 512: * 513: * Form H * C or H' * C where C = ( C1 ) 514: * ( C2 ) 515: * 516: LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 517: LASTC = ILADLC( LASTV, N, C, LDC ) 518: * 519: * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 520: * 521: * W := C2' 522: * 523: DO 190 J = 1, K 524: CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 525: $ WORK( 1, J ), 1 ) 526: 190 CONTINUE 527: * 528: * W := W * V2' 529: * 530: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 531: $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 532: $ WORK, LDWORK ) 533: IF( LASTV.GT.K ) THEN 534: * 535: * W := W + C1'*V1' 536: * 537: CALL DGEMM( 'Transpose', 'Transpose', 538: $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 539: $ ONE, WORK, LDWORK ) 540: END IF 541: * 542: * W := W * T' or W * T 543: * 544: CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 545: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 546: * 547: * C := C - V' * W' 548: * 549: IF( LASTV.GT.K ) THEN 550: * 551: * C1 := C1 - V1' * W' 552: * 553: CALL DGEMM( 'Transpose', 'Transpose', 554: $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 555: $ ONE, C, LDC ) 556: END IF 557: * 558: * W := W * V2 559: * 560: CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 561: $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 562: $ WORK, LDWORK ) 563: * 564: * C2 := C2 - W' 565: * 566: DO 210 J = 1, K 567: DO 200 I = 1, LASTC 568: C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 569: 200 CONTINUE 570: 210 CONTINUE 571: * 572: ELSE IF( LSAME( SIDE, 'R' ) ) THEN 573: * 574: * Form C * H or C * H' where C = ( C1 C2 ) 575: * 576: LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 577: LASTC = ILADLR( M, LASTV, C, LDC ) 578: * 579: * W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 580: * 581: * W := C2 582: * 583: DO 220 J = 1, K 584: CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1, 585: $ WORK( 1, J ), 1 ) 586: 220 CONTINUE 587: * 588: * W := W * V2' 589: * 590: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 591: $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 592: $ WORK, LDWORK ) 593: IF( LASTV.GT.K ) THEN 594: * 595: * W := W + C1 * V1' 596: * 597: CALL DGEMM( 'No transpose', 'Transpose', 598: $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 599: $ ONE, WORK, LDWORK ) 600: END IF 601: * 602: * W := W * T or W * T' 603: * 604: CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 605: $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 606: * 607: * C := C - W * V 608: * 609: IF( LASTV.GT.K ) THEN 610: * 611: * C1 := C1 - W * V1 612: * 613: CALL DGEMM( 'No transpose', 'No transpose', 614: $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 615: $ ONE, C, LDC ) 616: END IF 617: * 618: * W := W * V2 619: * 620: CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 621: $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 622: $ WORK, LDWORK ) 623: * 624: * C1 := C1 - W 625: * 626: DO 240 J = 1, K 627: DO 230 I = 1, LASTC 628: C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 629: 230 CONTINUE 630: 240 CONTINUE 631: * 632: END IF 633: * 634: END IF 635: END IF 636: * 637: RETURN 638: * 639: * End of DLARFB 640: * 641: END