Annotation of rpl/lapack/lapack/ztprfb.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZTPRFB
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTPRFB + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztprfb.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztprfb.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztprfb.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTPRFB( 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: * COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 30: * $ V( LDV, * ), WORK( LDWORK, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> ZTPRFB applies a complex "triangular-pentagonal" block reflector H or its
! 40: *> conjugate transpose H**H to a complex 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**H from the Left
! 52: *> = 'R': apply H or H**H from the Right
! 53: *> \endverbatim
! 54: *>
! 55: *> \param[in] TRANS
! 56: *> \verbatim
! 57: *> TRANS is CHARACTER*1
! 58: *> = 'N': apply H (No transpose)
! 59: *> = 'C': apply H**H (Conjugate 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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**H*C or C*H or C*H**H. 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 COMPLEX*16 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**H*C or C*H or C*H**H. 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 COMPLEX*16 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 complex16OTHERauxiliary
! 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 ZTPRFB( 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: COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 265: $ V( LDV, * ), WORK( LDWORK, * )
! 266: * ..
! 267: *
! 268: * ==========================================================================
! 269: *
! 270: * .. Parameters ..
! 271: COMPLEX*16 ONE, ZERO
! 272: PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,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 ZGEMM, ZTRMM
! 284: * ..
! 285: * .. Intrinsic Functions ..
! 286: INTRINSIC CONJG
! 287: * ..
! 288: * .. Executable Statements ..
! 289: *
! 290: * Quick return if possible
! 291: *
! 292: IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
! 293: *
! 294: IF( LSAME( STOREV, 'C' ) ) THEN
! 295: COLUMN = .TRUE.
! 296: ROW = .FALSE.
! 297: ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
! 298: COLUMN = .FALSE.
! 299: ROW = .TRUE.
! 300: ELSE
! 301: COLUMN = .FALSE.
! 302: ROW = .FALSE.
! 303: END IF
! 304: *
! 305: IF( LSAME( SIDE, 'L' ) ) THEN
! 306: LEFT = .TRUE.
! 307: RIGHT = .FALSE.
! 308: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
! 309: LEFT = .FALSE.
! 310: RIGHT = .TRUE.
! 311: ELSE
! 312: LEFT = .FALSE.
! 313: RIGHT = .FALSE.
! 314: END IF
! 315: *
! 316: IF( LSAME( DIRECT, 'F' ) ) THEN
! 317: FORWARD = .TRUE.
! 318: BACKWARD = .FALSE.
! 319: ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
! 320: FORWARD = .FALSE.
! 321: BACKWARD = .TRUE.
! 322: ELSE
! 323: FORWARD = .FALSE.
! 324: BACKWARD = .FALSE.
! 325: END IF
! 326: *
! 327: * ---------------------------------------------------------------------------
! 328: *
! 329: IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
! 330: *
! 331: * ---------------------------------------------------------------------------
! 332: *
! 333: * Let W = [ I ] (K-by-K)
! 334: * [ V ] (M-by-K)
! 335: *
! 336: * Form H C or H**H C where C = [ A ] (K-by-N)
! 337: * [ B ] (M-by-N)
! 338: *
! 339: * H = I - W T W**H or H**H = I - W T**H W**H
! 340: *
! 341: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
! 342: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
! 343: *
! 344: * ---------------------------------------------------------------------------
! 345: *
! 346: MP = MIN( M-L+1, M )
! 347: KP = MIN( L+1, K )
! 348: *
! 349: DO J = 1, N
! 350: DO I = 1, L
! 351: WORK( I, J ) = B( M-L+I, J )
! 352: END DO
! 353: END DO
! 354: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV,
! 355: $ WORK, LDWORK )
! 356: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
! 357: $ ONE, WORK, LDWORK )
! 358: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
! 359: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
! 360: *
! 361: DO J = 1, N
! 362: DO I = 1, K
! 363: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 364: END DO
! 365: END DO
! 366: *
! 367: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
! 368: $ WORK, LDWORK )
! 369: *
! 370: DO J = 1, N
! 371: DO I = 1, K
! 372: A( I, J ) = A( I, J ) - WORK( I, J )
! 373: END DO
! 374: END DO
! 375: *
! 376: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
! 377: $ ONE, B, LDB )
! 378: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
! 379: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
! 380: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
! 381: $ WORK, LDWORK )
! 382: DO J = 1, N
! 383: DO I = 1, L
! 384: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
! 385: END DO
! 386: END DO
! 387: *
! 388: * ---------------------------------------------------------------------------
! 389: *
! 390: ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
! 391: *
! 392: * ---------------------------------------------------------------------------
! 393: *
! 394: * Let W = [ I ] (K-by-K)
! 395: * [ V ] (N-by-K)
! 396: *
! 397: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
! 398: *
! 399: * H = I - W T W**H or H**H = I - W T**H W**H
! 400: *
! 401: * A = A - (A + B V) T or A = A - (A + B V) T**H
! 402: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
! 403: *
! 404: * ---------------------------------------------------------------------------
! 405: *
! 406: NP = MIN( N-L+1, N )
! 407: KP = MIN( L+1, K )
! 408: *
! 409: DO J = 1, L
! 410: DO I = 1, M
! 411: WORK( I, J ) = B( I, N-L+J )
! 412: END DO
! 413: END DO
! 414: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
! 415: $ WORK, LDWORK )
! 416: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
! 417: $ V, LDV, ONE, WORK, LDWORK )
! 418: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
! 419: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
! 420: *
! 421: DO J = 1, K
! 422: DO I = 1, M
! 423: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 424: END DO
! 425: END DO
! 426: *
! 427: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
! 428: $ WORK, LDWORK )
! 429: *
! 430: DO J = 1, K
! 431: DO I = 1, M
! 432: A( I, J ) = A( I, J ) - WORK( I, J )
! 433: END DO
! 434: END DO
! 435: *
! 436: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
! 437: $ V, LDV, ONE, B, LDB )
! 438: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
! 439: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
! 440: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( NP, 1 ), LDV,
! 441: $ WORK, LDWORK )
! 442: DO J = 1, L
! 443: DO I = 1, M
! 444: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
! 445: END DO
! 446: END DO
! 447: *
! 448: * ---------------------------------------------------------------------------
! 449: *
! 450: ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
! 451: *
! 452: * ---------------------------------------------------------------------------
! 453: *
! 454: * Let W = [ V ] (M-by-K)
! 455: * [ I ] (K-by-K)
! 456: *
! 457: * Form H C or H**H C where C = [ B ] (M-by-N)
! 458: * [ A ] (K-by-N)
! 459: *
! 460: * H = I - W T W**H or H**H = I - W T**H W**H
! 461: *
! 462: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
! 463: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
! 464: *
! 465: * ---------------------------------------------------------------------------
! 466: *
! 467: MP = MIN( L+1, M )
! 468: KP = MIN( K-L+1, K )
! 469: *
! 470: DO J = 1, N
! 471: DO I = 1, L
! 472: WORK( K-L+I, J ) = B( I, J )
! 473: END DO
! 474: END DO
! 475: *
! 476: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV,
! 477: $ WORK( KP, 1 ), LDWORK )
! 478: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
! 479: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
! 480: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV,
! 481: $ B, LDB, ZERO, WORK, LDWORK )
! 482: *
! 483: DO J = 1, N
! 484: DO I = 1, K
! 485: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 486: END DO
! 487: END DO
! 488: *
! 489: CALL ZTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
! 490: $ WORK, LDWORK )
! 491: *
! 492: DO J = 1, N
! 493: DO I = 1, K
! 494: A( I, J ) = A( I, J ) - WORK( I, J )
! 495: END DO
! 496: END DO
! 497: *
! 498: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
! 499: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
! 500: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
! 501: $ WORK, LDWORK, ONE, B, LDB )
! 502: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
! 503: $ WORK( KP, 1 ), LDWORK )
! 504: DO J = 1, N
! 505: DO I = 1, L
! 506: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
! 507: END DO
! 508: END DO
! 509: *
! 510: * ---------------------------------------------------------------------------
! 511: *
! 512: ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
! 513: *
! 514: * ---------------------------------------------------------------------------
! 515: *
! 516: * Let W = [ V ] (N-by-K)
! 517: * [ I ] (K-by-K)
! 518: *
! 519: * Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
! 520: *
! 521: * H = I - W T W**H or H**H = I - W T**H W**H
! 522: *
! 523: * A = A - (A + B V) T or A = A - (A + B V) T**H
! 524: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
! 525: *
! 526: * ---------------------------------------------------------------------------
! 527: *
! 528: NP = MIN( L+1, N )
! 529: KP = MIN( K-L+1, K )
! 530: *
! 531: DO J = 1, L
! 532: DO I = 1, M
! 533: WORK( I, K-L+J ) = B( I, J )
! 534: END DO
! 535: END DO
! 536: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
! 537: $ WORK( 1, KP ), LDWORK )
! 538: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
! 539: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
! 540: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
! 541: $ V, LDV, ZERO, WORK, LDWORK )
! 542: *
! 543: DO J = 1, K
! 544: DO I = 1, M
! 545: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 546: END DO
! 547: END DO
! 548: *
! 549: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
! 550: $ WORK, LDWORK )
! 551: *
! 552: DO J = 1, K
! 553: DO I = 1, M
! 554: A( I, J ) = A( I, J ) - WORK( I, J )
! 555: END DO
! 556: END DO
! 557: *
! 558: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
! 559: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
! 560: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK,
! 561: $ V, LDV, ONE, B, LDB )
! 562: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, KP ), LDV,
! 563: $ WORK( 1, KP ), LDWORK )
! 564: DO J = 1, L
! 565: DO I = 1, M
! 566: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
! 567: END DO
! 568: END DO
! 569: *
! 570: * ---------------------------------------------------------------------------
! 571: *
! 572: ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
! 573: *
! 574: * ---------------------------------------------------------------------------
! 575: *
! 576: * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
! 577: *
! 578: * Form H C or H**H C where C = [ A ] (K-by-N)
! 579: * [ B ] (M-by-N)
! 580: *
! 581: * H = I - W**H T W or H**H = I - W**H T**H W
! 582: *
! 583: * A = A - T (A + V B) or A = A - T**H (A + V B)
! 584: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
! 585: *
! 586: * ---------------------------------------------------------------------------
! 587: *
! 588: MP = MIN( M-L+1, M )
! 589: KP = MIN( L+1, K )
! 590: *
! 591: DO J = 1, N
! 592: DO I = 1, L
! 593: WORK( I, J ) = B( M-L+I, J )
! 594: END DO
! 595: END DO
! 596: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
! 597: $ WORK, LDB )
! 598: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
! 599: $ ONE, WORK, LDWORK )
! 600: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
! 601: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
! 602: *
! 603: DO J = 1, N
! 604: DO I = 1, K
! 605: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 606: END DO
! 607: END DO
! 608: *
! 609: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
! 610: $ WORK, LDWORK )
! 611: *
! 612: DO J = 1, N
! 613: DO I = 1, K
! 614: A( I, J ) = A( I, J ) - WORK( I, J )
! 615: END DO
! 616: END DO
! 617: *
! 618: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
! 619: $ ONE, B, LDB )
! 620: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
! 621: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
! 622: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, MP ), LDV,
! 623: $ WORK, LDWORK )
! 624: DO J = 1, N
! 625: DO I = 1, L
! 626: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
! 627: END DO
! 628: END DO
! 629: *
! 630: * ---------------------------------------------------------------------------
! 631: *
! 632: ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
! 633: *
! 634: * ---------------------------------------------------------------------------
! 635: *
! 636: * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
! 637: *
! 638: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
! 639: *
! 640: * H = I - W**H T W or H**H = I - W**H T**H W
! 641: *
! 642: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
! 643: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
! 644: *
! 645: * ---------------------------------------------------------------------------
! 646: *
! 647: NP = MIN( N-L+1, N )
! 648: KP = MIN( L+1, K )
! 649: *
! 650: DO J = 1, L
! 651: DO I = 1, M
! 652: WORK( I, J ) = B( I, N-L+J )
! 653: END DO
! 654: END DO
! 655: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV,
! 656: $ WORK, LDWORK )
! 657: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV,
! 658: $ ONE, WORK, LDWORK )
! 659: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB,
! 660: $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
! 661: *
! 662: DO J = 1, K
! 663: DO I = 1, M
! 664: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 665: END DO
! 666: END DO
! 667: *
! 668: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
! 669: $ WORK, LDWORK )
! 670: *
! 671: DO J = 1, K
! 672: DO I = 1, M
! 673: A( I, J ) = A( I, J ) - WORK( I, J )
! 674: END DO
! 675: END DO
! 676: *
! 677: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
! 678: $ V, LDV, ONE, B, LDB )
! 679: CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
! 680: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
! 681: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
! 682: $ WORK, LDWORK )
! 683: DO J = 1, L
! 684: DO I = 1, M
! 685: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
! 686: END DO
! 687: END DO
! 688: *
! 689: * ---------------------------------------------------------------------------
! 690: *
! 691: ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
! 692: *
! 693: * ---------------------------------------------------------------------------
! 694: *
! 695: * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
! 696: *
! 697: * Form H C or H**H C where C = [ B ] (M-by-N)
! 698: * [ A ] (K-by-N)
! 699: *
! 700: * H = I - W**H T W or H**H = I - W**H T**H W
! 701: *
! 702: * A = A - T (A + V B) or A = A - T**H (A + V B)
! 703: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
! 704: *
! 705: * ---------------------------------------------------------------------------
! 706: *
! 707: MP = MIN( L+1, M )
! 708: KP = MIN( K-L+1, K )
! 709: *
! 710: DO J = 1, N
! 711: DO I = 1, L
! 712: WORK( K-L+I, J ) = B( I, J )
! 713: END DO
! 714: END DO
! 715: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
! 716: $ WORK( KP, 1 ), LDWORK )
! 717: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
! 718: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
! 719: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
! 720: $ ZERO, WORK, LDWORK )
! 721: *
! 722: DO J = 1, N
! 723: DO I = 1, K
! 724: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 725: END DO
! 726: END DO
! 727: *
! 728: CALL ZTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
! 729: $ WORK, LDWORK )
! 730: *
! 731: DO J = 1, N
! 732: DO I = 1, K
! 733: A( I, J ) = A( I, J ) - WORK( I, J )
! 734: END DO
! 735: END DO
! 736: *
! 737: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
! 738: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
! 739: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV,
! 740: $ WORK, LDWORK, ONE, B, LDB )
! 741: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV,
! 742: $ WORK( KP, 1 ), LDWORK )
! 743: DO J = 1, N
! 744: DO I = 1, L
! 745: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
! 746: END DO
! 747: END DO
! 748: *
! 749: * ---------------------------------------------------------------------------
! 750: *
! 751: ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
! 752: *
! 753: * ---------------------------------------------------------------------------
! 754: *
! 755: * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
! 756: *
! 757: * Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N)
! 758: *
! 759: * H = I - W**H T W or H**H = I - W**H T**H W
! 760: *
! 761: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
! 762: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
! 763: *
! 764: * ---------------------------------------------------------------------------
! 765: *
! 766: NP = MIN( L+1, N )
! 767: KP = MIN( K-L+1, K )
! 768: *
! 769: DO J = 1, L
! 770: DO I = 1, M
! 771: WORK( I, K-L+J ) = B( I, J )
! 772: END DO
! 773: END DO
! 774: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV,
! 775: $ WORK( 1, KP ), LDWORK )
! 776: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB,
! 777: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
! 778: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, V, LDV,
! 779: $ ZERO, WORK, LDWORK )
! 780: *
! 781: DO J = 1, K
! 782: DO I = 1, M
! 783: WORK( I, J ) = WORK( I, J ) + A( I, J )
! 784: END DO
! 785: END DO
! 786: *
! 787: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
! 788: $ WORK, LDWORK )
! 789: *
! 790: DO J = 1, K
! 791: DO I = 1, M
! 792: A( I, J ) = A( I, J ) - WORK( I, J )
! 793: END DO
! 794: END DO
! 795: *
! 796: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
! 797: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
! 798: CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
! 799: $ V, LDV, ONE, B, LDB )
! 800: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
! 801: $ WORK( 1, KP ), LDWORK )
! 802: DO J = 1, L
! 803: DO I = 1, M
! 804: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
! 805: END DO
! 806: END DO
! 807: *
! 808: END IF
! 809: *
! 810: RETURN
! 811: *
! 812: * End of ZTPRFB
! 813: *
! 814: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztprfb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack