Annotation of rpl/lapack/lapack/dtprfb.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DTPRFB
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTPRFB + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfb.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfb.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfb.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
! 22: * V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER DIRECT, SIDE, STOREV, TRANS
! 26: * INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 30: * $ V( LDV, * ), WORK( LDWORK, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DTPRFB applies a real "triangular-pentagonal" block reflector H or its
! 40: *> transpose H**T to a real matrix C, which is composed of two
! 41: *> blocks A and B, either from the left or right.
! 42: *>
! 43: *> \endverbatim
! 44: *
! 45: * Arguments:
! 46: * ==========
! 47: *
! 48: *> \param[in] SIDE
! 49: *> \verbatim
! 50: *> SIDE is CHARACTER*1
! 51: *> = 'L': apply H or H**T from the Left
! 52: *> = 'R': apply H or H**T from the Right
! 53: *> \endverbatim
! 54: *>
! 55: *> \param[in] TRANS
! 56: *> \verbatim
! 57: *> TRANS is CHARACTER*1
! 58: *> = 'N': apply H (No transpose)
! 59: *> = 'T': apply H**T (Transpose)
! 60: *> \endverbatim
! 61: *>
! 62: *> \param[in] DIRECT
! 63: *> \verbatim
! 64: *> DIRECT is CHARACTER*1
! 65: *> Indicates how H is formed from a product of elementary
! 66: *> reflectors
! 67: *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
! 68: *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] STOREV
! 72: *> \verbatim
! 73: *> STOREV is CHARACTER*1
! 74: *> Indicates how the vectors which define the elementary
! 75: *> reflectors are stored:
! 76: *> = 'C': Columns
! 77: *> = 'R': Rows
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] M
! 81: *> \verbatim
! 82: *> M is INTEGER
! 83: *> The number of rows of the matrix B.
! 84: *> M >= 0.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in] N
! 88: *> \verbatim
! 89: *> N is INTEGER
! 90: *> The number of columns of the matrix B.
! 91: *> N >= 0.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] K
! 95: *> \verbatim
! 96: *> K is INTEGER
! 97: *> The order of the matrix T, i.e. the number of elementary
! 98: *> reflectors whose product defines the block reflector.
! 99: *> K >= 0.
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in] L
! 103: *> \verbatim
! 104: *> L is INTEGER
! 105: *> The order of the trapezoidal part of V.
! 106: *> K >= L >= 0. See Further Details.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in] V
! 110: *> \verbatim
! 111: *> V is DOUBLE PRECISION array, dimension
! 112: *> (LDV,K) if STOREV = 'C'
! 113: *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
! 114: *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
! 115: *> The pentagonal matrix V, which contains the elementary reflectors
! 116: *> H(1), H(2), ..., H(K). See Further Details.
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] LDV
! 120: *> \verbatim
! 121: *> LDV is INTEGER
! 122: *> The leading dimension of the array V.
! 123: *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
! 124: *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
! 125: *> if STOREV = 'R', LDV >= K.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] T
! 129: *> \verbatim
! 130: *> T is DOUBLE PRECISION array, dimension (LDT,K)
! 131: *> The triangular K-by-K matrix T in the representation of the
! 132: *> block reflector.
! 133: *> \endverbatim
! 134: *>
! 135: *> \param[in] LDT
! 136: *> \verbatim
! 137: *> LDT is INTEGER
! 138: *> The leading dimension of the array T.
! 139: *> LDT >= K.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[in,out] A
! 143: *> \verbatim
! 144: *> A is DOUBLE PRECISION array, dimension
! 145: *> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
! 146: *> On entry, the K-by-N or M-by-K matrix A.
! 147: *> On exit, A is overwritten by the corresponding block of
! 148: *> H*C or H**T*C or C*H or C*H**T. See Futher Details.
! 149: *> \endverbatim
! 150: *>
! 151: *> \param[in] LDA
! 152: *> \verbatim
! 153: *> LDA is INTEGER
! 154: *> The leading dimension of the array A.
! 155: *> If SIDE = 'L', LDC >= max(1,K);
! 156: *> If SIDE = 'R', LDC >= max(1,M).
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[in,out] B
! 160: *> \verbatim
! 161: *> B is DOUBLE PRECISION array, dimension (LDB,N)
! 162: *> On entry, the M-by-N matrix B.
! 163: *> On exit, B is overwritten by the corresponding block of
! 164: *> H*C or H**T*C or C*H or C*H**T. See Further Details.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[in] LDB
! 168: *> \verbatim
! 169: *> LDB is INTEGER
! 170: *> The leading dimension of the array B.
! 171: *> LDB >= max(1,M).
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[out] WORK
! 175: *> \verbatim
! 176: *> WORK is DOUBLE PRECISION array, dimension
! 177: *> (LDWORK,N) if SIDE = 'L',
! 178: *> (LDWORK,K) if SIDE = 'R'.
! 179: *> \endverbatim
! 180: *>
! 181: *> \param[in] LDWORK
! 182: *> \verbatim
! 183: *> LDWORK is INTEGER
! 184: *> The leading dimension of the array WORK.
! 185: *> If SIDE = 'L', LDWORK >= K;
! 186: *> if SIDE = 'R', LDWORK >= M.
! 187: *> \endverbatim
! 188: *
! 189: * Authors:
! 190: * ========
! 191: *
! 192: *> \author Univ. of Tennessee
! 193: *> \author Univ. of California Berkeley
! 194: *> \author Univ. of Colorado Denver
! 195: *> \author NAG Ltd.
! 196: *
! 197: *> \date November 2011
! 198: *
! 199: *> \ingroup doubleOTHERauxiliary
! 200: *
! 201: *> \par Further Details:
! 202: * =====================
! 203: *>
! 204: *> \verbatim
! 205: *>
! 206: *> The matrix C is a composite matrix formed from blocks A and B.
! 207: *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
! 208: *> and if SIDE = 'L', A is of size K-by-N.
! 209: *>
! 210: *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
! 211: *>
! 212: *> If SIDE = 'L' and DIRECT = 'F', C = [A]
! 213: *> [B].
! 214: *>
! 215: *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
! 216: *>
! 217: *> If SIDE = 'L' and DIRECT = 'B', C = [B]
! 218: *> [A].
! 219: *>
! 220: *> The pentagonal matrix V is composed of a rectangular block V1 and a
! 221: *> trapezoidal block V2. The size of the trapezoidal block is determined by
! 222: *> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
! 223: *> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
! 224: *>
! 225: *> If DIRECT = 'F' and STOREV = 'C': V = [V1]
! 226: *> [V2]
! 227: *> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
! 228: *>
! 229: *> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
! 230: *>
! 231: *> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
! 232: *>
! 233: *> If DIRECT = 'B' and STOREV = 'C': V = [V2]
! 234: *> [V1]
! 235: *> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
! 236: *>
! 237: *> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
! 238: *>
! 239: *> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
! 240: *>
! 241: *> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
! 242: *>
! 243: *> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
! 244: *>
! 245: *> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
! 246: *>
! 247: *> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
! 248: *> \endverbatim
! 249: *>
! 250: * =====================================================================
! 251: SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
! 252: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
! 253: *
! 254: * -- LAPACK auxiliary routine (version 3.4.0) --
! 255: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 256: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 257: * November 2011
! 258: *
! 259: * .. Scalar Arguments ..
! 260: CHARACTER DIRECT, SIDE, STOREV, TRANS
! 261: INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
! 262: * ..
! 263: * .. Array Arguments ..
! 264: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 265: $ V( LDV, * ), WORK( LDWORK, * )
! 266: * ..
! 267: *
! 268: * ==========================================================================
! 269: *
! 270: * .. Parameters ..
! 271: DOUBLE PRECISION ONE, ZERO
! 272: PARAMETER ( ONE = 1.0, ZERO = 0.0 )
! 273: * ..
! 274: * .. Local Scalars ..
! 275: INTEGER I, J, MP, NP, KP
! 276: LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
! 277: * ..
! 278: * .. External Functions ..
! 279: LOGICAL LSAME
! 280: EXTERNAL LSAME
! 281: * ..
! 282: * .. External Subroutines ..
! 283: EXTERNAL DGEMM, DTRMM
! 284: * ..
! 285: * .. Executable Statements ..
! 286: *
! 287: * Quick return if possible
! 288: *
! 289: IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
! 290: *
! 291: IF( LSAME( STOREV, 'C' ) ) THEN
! 292: COLUMN = .TRUE.
! 293: ROW = .FALSE.
! 294: ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
! 295: COLUMN = .FALSE.
! 296: ROW = .TRUE.
! 297: ELSE
! 298: COLUMN = .FALSE.
! 299: ROW = .FALSE.
! 300: END IF
! 301: *
! 302: IF( LSAME( SIDE, 'L' ) ) THEN
! 303: LEFT = .TRUE.
! 304: RIGHT = .FALSE.
! 305: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
! 306: LEFT = .FALSE.
! 307: RIGHT = .TRUE.
! 308: ELSE
! 309: LEFT = .FALSE.
! 310: RIGHT = .FALSE.
! 311: END IF
! 312: *
! 313: IF( LSAME( DIRECT, 'F' ) ) THEN
! 314: FORWARD = .TRUE.
! 315: BACKWARD = .FALSE.
! 316: ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
! 317: FORWARD = .FALSE.
! 318: BACKWARD = .TRUE.
! 319: ELSE
! 320: FORWARD = .FALSE.
! 321: BACKWARD = .FALSE.
! 322: END IF
! 323: *
! 324: * ---------------------------------------------------------------------------
! 325: *
! 326: IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
! 327: *
! 328: * ---------------------------------------------------------------------------
! 329: *
! 330: * Let W = [ I ] (K-by-K)
! 331: * [ V ] (M-by-K)
! 332: *
! 333: * Form H C or H**T C where C = [ A ] (K-by-N)
! 334: * [ B ] (M-by-N)
! 335: *
! 336: * H = I - W T W**T or H**T = I - W T**T W**T
! 337: *
! 338: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
! 339: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
! 340: *
! 341: * ---------------------------------------------------------------------------
! 342: *
! 343: MP = MIN( M-L+1, M )
! 344: KP = MIN( L+1, K )
! 345: *
! 346: DO J = 1, N
! 347: DO I = 1, L
! 348: WORK( I, J ) = B( M-L+I, J )
! 349: END DO
! 350: END DO
! 351: CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
! 352: $ WORK, LDWORK )
! 353: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
! 354: $ ONE, WORK, LDWORK )
! 355: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
! 356: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
! 357: *
! 358: DO J = 1, N
! 359: DO I = 1, K
! 360: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 361: END DO
! 362: END DO
! 363: *
! 364: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
! 365: $ WORK, LDWORK )
! 366: *
! 367: DO J = 1, N
! 368: DO I = 1, K
! 369: A( I, J ) = A( I, J ) - WORK( I, J )
! 370: END DO
! 371: END DO
! 372: *
! 373: CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
! 374: $ ONE, B, LDB )
! 375: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
! 376: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
! 377: CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
! 378: $ WORK, LDWORK )
! 379: DO J = 1, N
! 380: DO I = 1, L
! 381: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
! 382: END DO
! 383: END DO
! 384: *
! 385: * ---------------------------------------------------------------------------
! 386: *
! 387: ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
! 388: *
! 389: * ---------------------------------------------------------------------------
! 390: *
! 391: * Let W = [ I ] (K-by-K)
! 392: * [ V ] (N-by-K)
! 393: *
! 394: * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
! 395: *
! 396: * H = I - W T W**T or H**T = I - W T**T W**T
! 397: *
! 398: * A = A - (A + B V) T or A = A - (A + B V) T**T
! 399: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
! 400: *
! 401: * ---------------------------------------------------------------------------
! 402: *
! 403: NP = MIN( N-L+1, N )
! 404: KP = MIN( L+1, K )
! 405: *
! 406: DO J = 1, L
! 407: DO I = 1, M
! 408: WORK( I, J ) = B( I, N-L+J )
! 409: END DO
! 410: END DO
! 411: CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
! 412: $ WORK, LDWORK )
! 413: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
! 414: $ V, LDV, ONE, WORK, LDWORK )
! 415: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
! 416: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
! 417: *
! 418: DO J = 1, K
! 419: DO I = 1, M
! 420: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 421: END DO
! 422: END DO
! 423: *
! 424: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
! 425: $ WORK, LDWORK )
! 426: *
! 427: DO J = 1, K
! 428: DO I = 1, M
! 429: A( I, J ) = A( I, J ) - WORK( I, J )
! 430: END DO
! 431: END DO
! 432: *
! 433: CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
! 434: $ V, LDV, ONE, B, LDB )
! 435: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
! 436: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
! 437: CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV,
! 438: $ WORK, LDWORK )
! 439: DO J = 1, L
! 440: DO I = 1, M
! 441: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
! 442: END DO
! 443: END DO
! 444: *
! 445: * ---------------------------------------------------------------------------
! 446: *
! 447: ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
! 448: *
! 449: * ---------------------------------------------------------------------------
! 450: *
! 451: * Let W = [ V ] (M-by-K)
! 452: * [ I ] (K-by-K)
! 453: *
! 454: * Form H C or H**T C where C = [ B ] (M-by-N)
! 455: * [ A ] (K-by-N)
! 456: *
! 457: * H = I - W T W**T or H**T = I - W T**T W**T
! 458: *
! 459: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
! 460: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
! 461: *
! 462: * ---------------------------------------------------------------------------
! 463: *
! 464: MP = MIN( L+1, M )
! 465: KP = MIN( K-L+1, K )
! 466: *
! 467: DO J = 1, N
! 468: DO I = 1, L
! 469: WORK( K-L+I, J ) = B( I, J )
! 470: END DO
! 471: END DO
! 472: *
! 473: CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
! 474: $ WORK( KP, 1 ), LDWORK )
! 475: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
! 476: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
! 477: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV,
! 478: $ B, LDB, ZERO, WORK, LDWORK )
! 479: *
! 480: DO J = 1, N
! 481: DO I = 1, K
! 482: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 483: END DO
! 484: END DO
! 485: *
! 486: CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
! 487: $ WORK, LDWORK )
! 488: *
! 489: DO J = 1, N
! 490: DO I = 1, K
! 491: A( I, J ) = A( I, J ) - WORK( I, J )
! 492: END DO
! 493: END DO
! 494: *
! 495: CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
! 496: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
! 497: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
! 498: $ WORK, LDWORK, ONE, B, LDB )
! 499: CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
! 500: $ WORK( KP, 1 ), LDWORK )
! 501: DO J = 1, N
! 502: DO I = 1, L
! 503: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
! 504: END DO
! 505: END DO
! 506: *
! 507: * ---------------------------------------------------------------------------
! 508: *
! 509: ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
! 510: *
! 511: * ---------------------------------------------------------------------------
! 512: *
! 513: * Let W = [ V ] (N-by-K)
! 514: * [ I ] (K-by-K)
! 515: *
! 516: * Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K)
! 517: *
! 518: * H = I - W T W**T or H**T = I - W T**T W**T
! 519: *
! 520: * A = A - (A + B V) T or A = A - (A + B V) T**T
! 521: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
! 522: *
! 523: * ---------------------------------------------------------------------------
! 524: *
! 525: NP = MIN( L+1, N )
! 526: KP = MIN( K-L+1, K )
! 527: *
! 528: DO J = 1, L
! 529: DO I = 1, M
! 530: WORK( I, K-L+J ) = B( I, J )
! 531: END DO
! 532: END DO
! 533: CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
! 534: $ WORK( 1, KP ), LDWORK )
! 535: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
! 536: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
! 537: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
! 538: $ V, LDV, ZERO, WORK, LDWORK )
! 539: *
! 540: DO J = 1, K
! 541: DO I = 1, M
! 542: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 543: END DO
! 544: END DO
! 545: *
! 546: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
! 547: $ WORK, LDWORK )
! 548: *
! 549: DO J = 1, K
! 550: DO I = 1, M
! 551: A( I, J ) = A( I, J ) - WORK( I, J )
! 552: END DO
! 553: END DO
! 554: *
! 555: CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
! 556: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
! 557: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
! 558: $ V, LDV, ONE, B, LDB )
! 559: CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV,
! 560: $ WORK( 1, KP ), LDWORK )
! 561: DO J = 1, L
! 562: DO I = 1, M
! 563: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
! 564: END DO
! 565: END DO
! 566: *
! 567: * ---------------------------------------------------------------------------
! 568: *
! 569: ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
! 570: *
! 571: * ---------------------------------------------------------------------------
! 572: *
! 573: * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
! 574: *
! 575: * Form H C or H**T C where C = [ A ] (K-by-N)
! 576: * [ B ] (M-by-N)
! 577: *
! 578: * H = I - W**T T W or H**T = I - W**T T**T W
! 579: *
! 580: * A = A - T (A + V B) or A = A - T**T (A + V B)
! 581: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
! 582: *
! 583: * ---------------------------------------------------------------------------
! 584: *
! 585: MP = MIN( M-L+1, M )
! 586: KP = MIN( L+1, K )
! 587: *
! 588: DO J = 1, N
! 589: DO I = 1, L
! 590: WORK( I, J ) = B( M-L+I, J )
! 591: END DO
! 592: END DO
! 593: CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
! 594: $ WORK, LDB )
! 595: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
! 596: $ ONE, WORK, LDWORK )
! 597: CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
! 598: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
! 599: *
! 600: DO J = 1, N
! 601: DO I = 1, K
! 602: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 603: END DO
! 604: END DO
! 605: *
! 606: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
! 607: $ WORK, LDWORK )
! 608: *
! 609: DO J = 1, N
! 610: DO I = 1, K
! 611: A( I, J ) = A( I, J ) - WORK( I, J )
! 612: END DO
! 613: END DO
! 614: *
! 615: CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
! 616: $ ONE, B, LDB )
! 617: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
! 618: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
! 619: CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV,
! 620: $ WORK, LDWORK )
! 621: DO J = 1, N
! 622: DO I = 1, L
! 623: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
! 624: END DO
! 625: END DO
! 626: *
! 627: * ---------------------------------------------------------------------------
! 628: *
! 629: ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
! 630: *
! 631: * ---------------------------------------------------------------------------
! 632: *
! 633: * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
! 634: *
! 635: * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
! 636: *
! 637: * H = I - W**T T W or H**T = I - W**T T**T W
! 638: *
! 639: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
! 640: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
! 641: *
! 642: * ---------------------------------------------------------------------------
! 643: *
! 644: NP = MIN( N-L+1, N )
! 645: KP = MIN( L+1, K )
! 646: *
! 647: DO J = 1, L
! 648: DO I = 1, M
! 649: WORK( I, J ) = B( I, N-L+J )
! 650: END DO
! 651: END DO
! 652: CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
! 653: $ WORK, LDWORK )
! 654: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
! 655: $ ONE, WORK, LDWORK )
! 656: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB,
! 657: $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
! 658: *
! 659: DO J = 1, K
! 660: DO I = 1, M
! 661: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 662: END DO
! 663: END DO
! 664: *
! 665: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
! 666: $ WORK, LDWORK )
! 667: *
! 668: DO J = 1, K
! 669: DO I = 1, M
! 670: A( I, J ) = A( I, J ) - WORK( I, J )
! 671: END DO
! 672: END DO
! 673: *
! 674: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
! 675: $ V, LDV, ONE, B, LDB )
! 676: CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
! 677: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
! 678: CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
! 679: $ WORK, LDWORK )
! 680: DO J = 1, L
! 681: DO I = 1, M
! 682: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
! 683: END DO
! 684: END DO
! 685: *
! 686: * ---------------------------------------------------------------------------
! 687: *
! 688: ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
! 689: *
! 690: * ---------------------------------------------------------------------------
! 691: *
! 692: * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
! 693: *
! 694: * Form H C or H**T C where C = [ B ] (M-by-N)
! 695: * [ A ] (K-by-N)
! 696: *
! 697: * H = I - W**T T W or H**T = I - W**T T**T W
! 698: *
! 699: * A = A - T (A + V B) or A = A - T**T (A + V B)
! 700: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
! 701: *
! 702: * ---------------------------------------------------------------------------
! 703: *
! 704: MP = MIN( L+1, M )
! 705: KP = MIN( K-L+1, K )
! 706: *
! 707: DO J = 1, N
! 708: DO I = 1, L
! 709: WORK( K-L+I, J ) = B( I, J )
! 710: END DO
! 711: END DO
! 712: CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
! 713: $ WORK( KP, 1 ), LDWORK )
! 714: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
! 715: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
! 716: CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
! 717: $ ZERO, WORK, LDWORK )
! 718: *
! 719: DO J = 1, N
! 720: DO I = 1, K
! 721: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 722: END DO
! 723: END DO
! 724: *
! 725: CALL DTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
! 726: $ WORK, LDWORK )
! 727: *
! 728: DO J = 1, N
! 729: DO I = 1, K
! 730: A( I, J ) = A( I, J ) - WORK( I, J )
! 731: END DO
! 732: END DO
! 733: *
! 734: CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
! 735: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
! 736: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV,
! 737: $ WORK, LDWORK, ONE, B, LDB )
! 738: CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV,
! 739: $ WORK( KP, 1 ), LDWORK )
! 740: DO J = 1, N
! 741: DO I = 1, L
! 742: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
! 743: END DO
! 744: END DO
! 745: *
! 746: * ---------------------------------------------------------------------------
! 747: *
! 748: ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
! 749: *
! 750: * ---------------------------------------------------------------------------
! 751: *
! 752: * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
! 753: *
! 754: * Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N)
! 755: *
! 756: * H = I - W**T T W or H**T = I - W**T T**T W
! 757: *
! 758: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
! 759: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
! 760: *
! 761: * ---------------------------------------------------------------------------
! 762: *
! 763: NP = MIN( L+1, N )
! 764: KP = MIN( K-L+1, K )
! 765: *
! 766: DO J = 1, L
! 767: DO I = 1, M
! 768: WORK( I, K-L+J ) = B( I, J )
! 769: END DO
! 770: END DO
! 771: CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
! 772: $ WORK( 1, KP ), LDWORK )
! 773: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
! 774: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
! 775: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV,
! 776: $ ZERO, WORK, LDWORK )
! 777: *
! 778: DO J = 1, K
! 779: DO I = 1, M
! 780: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 781: END DO
! 782: END DO
! 783: *
! 784: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
! 785: $ WORK, LDWORK )
! 786: *
! 787: DO J = 1, K
! 788: DO I = 1, M
! 789: A( I, J ) = A( I, J ) - WORK( I, J )
! 790: END DO
! 791: END DO
! 792: *
! 793: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
! 794: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
! 795: CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
! 796: $ V, LDV, ONE, B, LDB )
! 797: CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
! 798: $ WORK( 1, KP ), LDWORK )
! 799: DO J = 1, L
! 800: DO I = 1, M
! 801: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
! 802: END DO
! 803: END DO
! 804: *
! 805: END IF
! 806: *
! 807: RETURN
! 808: *
! 809: * End of DTPRFB
! 810: *
! 811: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtprfb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack