Annotation of rpl/lapack/lapack/dlarfb_gett.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DLARFB_GETT
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLARFB_GETT + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb_gett.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb_gett.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb_gett.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
! 22: * $ WORK, LDWORK )
! 23: * IMPLICIT NONE
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER IDENT
! 27: * INTEGER K, LDA, LDB, LDT, LDWORK, M, N
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 31: * $ WORK( LDWORK, * )
! 32: * ..
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLARFB_GETT applies a real Householder block reflector H from the
! 40: *> left to a real (K+M)-by-N "triangular-pentagonal" matrix
! 41: *> composed of two block matrices: an upper trapezoidal K-by-N matrix A
! 42: *> stored in the array A, and a rectangular M-by-(N-K) matrix B, stored
! 43: *> in the array B. The block reflector H is stored in a compact
! 44: *> WY-representation, where the elementary reflectors are in the
! 45: *> arrays A, B and T. See Further Details section.
! 46: *> \endverbatim
! 47: *
! 48: * Arguments:
! 49: * ==========
! 50: *
! 51: *> \param[in] IDENT
! 52: *> \verbatim
! 53: *> IDENT is CHARACTER*1
! 54: *> If IDENT = not 'I', or not 'i', then V1 is unit
! 55: *> lower-triangular and stored in the left K-by-K block of
! 56: *> the input matrix A,
! 57: *> If IDENT = 'I' or 'i', then V1 is an identity matrix and
! 58: *> not stored.
! 59: *> See Further Details section.
! 60: *> \endverbatim
! 61: *>
! 62: *> \param[in] M
! 63: *> \verbatim
! 64: *> M is INTEGER
! 65: *> The number of rows of the matrix B.
! 66: *> M >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in] N
! 70: *> \verbatim
! 71: *> N is INTEGER
! 72: *> The number of columns of the matrices A and B.
! 73: *> N >= 0.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[in] K
! 77: *> \verbatim
! 78: *> K is INTEGER
! 79: *> The number or rows of the matrix A.
! 80: *> K is also order of the matrix T, i.e. the number of
! 81: *> elementary reflectors whose product defines the block
! 82: *> reflector. 0 <= K <= N.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] T
! 86: *> \verbatim
! 87: *> T is DOUBLE PRECISION array, dimension (LDT,K)
! 88: *> The upper-triangular K-by-K matrix T in the representation
! 89: *> of the block reflector.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] LDT
! 93: *> \verbatim
! 94: *> LDT is INTEGER
! 95: *> The leading dimension of the array T. LDT >= K.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in,out] A
! 99: *> \verbatim
! 100: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 101: *>
! 102: *> On entry:
! 103: *> a) In the K-by-N upper-trapezoidal part A: input matrix A.
! 104: *> b) In the columns below the diagonal: columns of V1
! 105: *> (ones are not stored on the diagonal).
! 106: *>
! 107: *> On exit:
! 108: *> A is overwritten by rectangular K-by-N product H*A.
! 109: *>
! 110: *> See Further Details section.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] LDA
! 114: *> \verbatim
! 115: *> LDB is INTEGER
! 116: *> The leading dimension of the array A. LDA >= max(1,K).
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in,out] B
! 120: *> \verbatim
! 121: *> B is DOUBLE PRECISION array, dimension (LDB,N)
! 122: *>
! 123: *> On entry:
! 124: *> a) In the M-by-(N-K) right block: input matrix B.
! 125: *> b) In the M-by-N left block: columns of V2.
! 126: *>
! 127: *> On exit:
! 128: *> B is overwritten by rectangular M-by-N product H*B.
! 129: *>
! 130: *> See Further Details section.
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in] LDB
! 134: *> \verbatim
! 135: *> LDB is INTEGER
! 136: *> The leading dimension of the array B. LDB >= max(1,M).
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[out] WORK
! 140: *> \verbatim
! 141: *> WORK is DOUBLE PRECISION array,
! 142: *> dimension (LDWORK,max(K,N-K))
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[in] LDWORK
! 146: *> \verbatim
! 147: *> LDWORK is INTEGER
! 148: *> The leading dimension of the array WORK. LDWORK>=max(1,K).
! 149: *>
! 150: *> \endverbatim
! 151: *
! 152: * Authors:
! 153: * ========
! 154: *
! 155: *> \author Univ. of Tennessee
! 156: *> \author Univ. of California Berkeley
! 157: *> \author Univ. of Colorado Denver
! 158: *> \author NAG Ltd.
! 159: *
! 160: *> \ingroup doubleOTHERauxiliary
! 161: *
! 162: *> \par Contributors:
! 163: * ==================
! 164: *>
! 165: *> \verbatim
! 166: *>
! 167: *> November 2020, Igor Kozachenko,
! 168: *> Computer Science Division,
! 169: *> University of California, Berkeley
! 170: *>
! 171: *> \endverbatim
! 172: *
! 173: *> \par Further Details:
! 174: * =====================
! 175: *>
! 176: *> \verbatim
! 177: *>
! 178: *> (1) Description of the Algebraic Operation.
! 179: *>
! 180: *> The matrix A is a K-by-N matrix composed of two column block
! 181: *> matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
! 182: *> A = ( A1, A2 ).
! 183: *> The matrix B is an M-by-N matrix composed of two column block
! 184: *> matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
! 185: *> B = ( B1, B2 ).
! 186: *>
! 187: *> Perform the operation:
! 188: *>
! 189: *> ( A_out ) := H * ( A_in ) = ( I - V * T * V**T ) * ( A_in ) =
! 190: *> ( B_out ) ( B_in ) ( B_in )
! 191: *> = ( I - ( V1 ) * T * ( V1**T, V2**T ) ) * ( A_in )
! 192: *> ( V2 ) ( B_in )
! 193: *> On input:
! 194: *>
! 195: *> a) ( A_in ) consists of two block columns:
! 196: *> ( B_in )
! 197: *>
! 198: *> ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
! 199: *> ( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )),
! 200: *>
! 201: *> where the column blocks are:
! 202: *>
! 203: *> ( A1_in ) is a K-by-K upper-triangular matrix stored in the
! 204: *> upper triangular part of the array A(1:K,1:K).
! 205: *> ( B1_in ) is an M-by-K rectangular ZERO matrix and not stored.
! 206: *>
! 207: *> ( A2_in ) is a K-by-(N-K) rectangular matrix stored
! 208: *> in the array A(1:K,K+1:N).
! 209: *> ( B2_in ) is an M-by-(N-K) rectangular matrix stored
! 210: *> in the array B(1:M,K+1:N).
! 211: *>
! 212: *> b) V = ( V1 )
! 213: *> ( V2 )
! 214: *>
! 215: *> where:
! 216: *> 1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
! 217: *> 2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
! 218: *> stored in the lower-triangular part of the array
! 219: *> A(1:K,1:K) (ones are not stored),
! 220: *> and V2 is an M-by-K rectangular stored the array B(1:M,1:K),
! 221: *> (because on input B1_in is a rectangular zero
! 222: *> matrix that is not stored and the space is
! 223: *> used to store V2).
! 224: *>
! 225: *> c) T is a K-by-K upper-triangular matrix stored
! 226: *> in the array T(1:K,1:K).
! 227: *>
! 228: *> On output:
! 229: *>
! 230: *> a) ( A_out ) consists of two block columns:
! 231: *> ( B_out )
! 232: *>
! 233: *> ( A_out ) = (( A1_out ) ( A2_out ))
! 234: *> ( B_out ) (( B1_out ) ( B2_out )),
! 235: *>
! 236: *> where the column blocks are:
! 237: *>
! 238: *> ( A1_out ) is a K-by-K square matrix, or a K-by-K
! 239: *> upper-triangular matrix, if V1 is an
! 240: *> identity matrix. AiOut is stored in
! 241: *> the array A(1:K,1:K).
! 242: *> ( B1_out ) is an M-by-K rectangular matrix stored
! 243: *> in the array B(1:M,K:N).
! 244: *>
! 245: *> ( A2_out ) is a K-by-(N-K) rectangular matrix stored
! 246: *> in the array A(1:K,K+1:N).
! 247: *> ( B2_out ) is an M-by-(N-K) rectangular matrix stored
! 248: *> in the array B(1:M,K+1:N).
! 249: *>
! 250: *>
! 251: *> The operation above can be represented as the same operation
! 252: *> on each block column:
! 253: *>
! 254: *> ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**T ) * ( A1_in )
! 255: *> ( B1_out ) ( 0 ) ( 0 )
! 256: *>
! 257: *> ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**T ) * ( A2_in )
! 258: *> ( B2_out ) ( B2_in ) ( B2_in )
! 259: *>
! 260: *> If IDENT != 'I':
! 261: *>
! 262: *> The computation for column block 1:
! 263: *>
! 264: *> A1_out: = A1_in - V1*T*(V1**T)*A1_in
! 265: *>
! 266: *> B1_out: = - V2*T*(V1**T)*A1_in
! 267: *>
! 268: *> The computation for column block 2, which exists if N > K:
! 269: *>
! 270: *> A2_out: = A2_in - V1*T*( (V1**T)*A2_in + (V2**T)*B2_in )
! 271: *>
! 272: *> B2_out: = B2_in - V2*T*( (V1**T)*A2_in + (V2**T)*B2_in )
! 273: *>
! 274: *> If IDENT == 'I':
! 275: *>
! 276: *> The operation for column block 1:
! 277: *>
! 278: *> A1_out: = A1_in - V1*T**A1_in
! 279: *>
! 280: *> B1_out: = - V2*T**A1_in
! 281: *>
! 282: *> The computation for column block 2, which exists if N > K:
! 283: *>
! 284: *> A2_out: = A2_in - T*( A2_in + (V2**T)*B2_in )
! 285: *>
! 286: *> B2_out: = B2_in - V2*T*( A2_in + (V2**T)*B2_in )
! 287: *>
! 288: *> (2) Description of the Algorithmic Computation.
! 289: *>
! 290: *> In the first step, we compute column block 2, i.e. A2 and B2.
! 291: *> Here, we need to use the K-by-(N-K) rectangular workspace
! 292: *> matrix W2 that is of the same size as the matrix A2.
! 293: *> W2 is stored in the array WORK(1:K,1:(N-K)).
! 294: *>
! 295: *> In the second step, we compute column block 1, i.e. A1 and B1.
! 296: *> Here, we need to use the K-by-K square workspace matrix W1
! 297: *> that is of the same size as the as the matrix A1.
! 298: *> W1 is stored in the array WORK(1:K,1:K).
! 299: *>
! 300: *> NOTE: Hence, in this routine, we need the workspace array WORK
! 301: *> only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from
! 302: *> the first step and W1 from the second step.
! 303: *>
! 304: *> Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',
! 305: *> more computations than in the Case (B).
! 306: *>
! 307: *> if( IDENT != 'I' ) then
! 308: *> if ( N > K ) then
! 309: *> (First Step - column block 2)
! 310: *> col2_(1) W2: = A2
! 311: *> col2_(2) W2: = (V1**T) * W2 = (unit_lower_tr_of_(A1)**T) * W2
! 312: *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
! 313: *> col2_(4) W2: = T * W2
! 314: *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
! 315: *> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
! 316: *> col2_(7) A2: = A2 - W2
! 317: *> else
! 318: *> (Second Step - column block 1)
! 319: *> col1_(1) W1: = A1
! 320: *> col1_(2) W1: = (V1**T) * W1 = (unit_lower_tr_of_(A1)**T) * W1
! 321: *> col1_(3) W1: = T * W1
! 322: *> col1_(4) B1: = - V2 * W1 = - B1 * W1
! 323: *> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
! 324: *> col1_(6) square A1: = A1 - W1
! 325: *> end if
! 326: *> end if
! 327: *>
! 328: *> Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',
! 329: *> less computations than in the Case (A)
! 330: *>
! 331: *> if( IDENT == 'I' ) then
! 332: *> if ( N > K ) then
! 333: *> (First Step - column block 2)
! 334: *> col2_(1) W2: = A2
! 335: *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
! 336: *> col2_(4) W2: = T * W2
! 337: *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
! 338: *> col2_(7) A2: = A2 - W2
! 339: *> else
! 340: *> (Second Step - column block 1)
! 341: *> col1_(1) W1: = A1
! 342: *> col1_(3) W1: = T * W1
! 343: *> col1_(4) B1: = - V2 * W1 = - B1 * W1
! 344: *> col1_(6) upper-triangular_of_(A1): = A1 - W1
! 345: *> end if
! 346: *> end if
! 347: *>
! 348: *> Combine these cases (A) and (B) together, this is the resulting
! 349: *> algorithm:
! 350: *>
! 351: *> if ( N > K ) then
! 352: *>
! 353: *> (First Step - column block 2)
! 354: *>
! 355: *> col2_(1) W2: = A2
! 356: *> if( IDENT != 'I' ) then
! 357: *> col2_(2) W2: = (V1**T) * W2
! 358: *> = (unit_lower_tr_of_(A1)**T) * W2
! 359: *> end if
! 360: *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2]
! 361: *> col2_(4) W2: = T * W2
! 362: *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
! 363: *> if( IDENT != 'I' ) then
! 364: *> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
! 365: *> end if
! 366: *> col2_(7) A2: = A2 - W2
! 367: *>
! 368: *> else
! 369: *>
! 370: *> (Second Step - column block 1)
! 371: *>
! 372: *> col1_(1) W1: = A1
! 373: *> if( IDENT != 'I' ) then
! 374: *> col1_(2) W1: = (V1**T) * W1
! 375: *> = (unit_lower_tr_of_(A1)**T) * W1
! 376: *> end if
! 377: *> col1_(3) W1: = T * W1
! 378: *> col1_(4) B1: = - V2 * W1 = - B1 * W1
! 379: *> if( IDENT != 'I' ) then
! 380: *> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
! 381: *> col1_(6_a) below_diag_of_(A1): = - below_diag_of_(W1)
! 382: *> end if
! 383: *> col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)
! 384: *>
! 385: *> end if
! 386: *>
! 387: *> \endverbatim
! 388: *>
! 389: * =====================================================================
! 390: SUBROUTINE DLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
! 391: $ WORK, LDWORK )
! 392: IMPLICIT NONE
! 393: *
! 394: * -- LAPACK auxiliary routine --
! 395: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 396: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 397: *
! 398: * .. Scalar Arguments ..
! 399: CHARACTER IDENT
! 400: INTEGER K, LDA, LDB, LDT, LDWORK, M, N
! 401: * ..
! 402: * .. Array Arguments ..
! 403: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
! 404: $ WORK( LDWORK, * )
! 405: * ..
! 406: *
! 407: * =====================================================================
! 408: *
! 409: * .. Parameters ..
! 410: DOUBLE PRECISION ONE, ZERO
! 411: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 412: * ..
! 413: * .. Local Scalars ..
! 414: LOGICAL LNOTIDENT
! 415: INTEGER I, J
! 416: * ..
! 417: * .. EXTERNAL FUNCTIONS ..
! 418: LOGICAL LSAME
! 419: EXTERNAL LSAME
! 420: * ..
! 421: * .. External Subroutines ..
! 422: EXTERNAL DCOPY, DGEMM, DTRMM
! 423: * ..
! 424: * .. Executable Statements ..
! 425: *
! 426: * Quick return if possible
! 427: *
! 428: IF( M.LT.0 .OR. N.LE.0 .OR. K.EQ.0 .OR. K.GT.N )
! 429: $ RETURN
! 430: *
! 431: LNOTIDENT = .NOT.LSAME( IDENT, 'I' )
! 432: *
! 433: * ------------------------------------------------------------------
! 434: *
! 435: * First Step. Computation of the Column Block 2:
! 436: *
! 437: * ( A2 ) := H * ( A2 )
! 438: * ( B2 ) ( B2 )
! 439: *
! 440: * ------------------------------------------------------------------
! 441: *
! 442: IF( N.GT.K ) THEN
! 443: *
! 444: * col2_(1) Compute W2: = A2. Therefore, copy A2 = A(1:K, K+1:N)
! 445: * into W2=WORK(1:K, 1:N-K) column-by-column.
! 446: *
! 447: DO J = 1, N-K
! 448: CALL DCOPY( K, A( 1, K+J ), 1, WORK( 1, J ), 1 )
! 449: END DO
! 450:
! 451: IF( LNOTIDENT ) THEN
! 452: *
! 453: * col2_(2) Compute W2: = (V1**T) * W2 = (A1**T) * W2,
! 454: * V1 is not an identy matrix, but unit lower-triangular
! 455: * V1 stored in A1 (diagonal ones are not stored).
! 456: *
! 457: *
! 458: CALL DTRMM( 'L', 'L', 'T', 'U', K, N-K, ONE, A, LDA,
! 459: $ WORK, LDWORK )
! 460: END IF
! 461: *
! 462: * col2_(3) Compute W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
! 463: * V2 stored in B1.
! 464: *
! 465: IF( M.GT.0 ) THEN
! 466: CALL DGEMM( 'T', 'N', K, N-K, M, ONE, B, LDB,
! 467: $ B( 1, K+1 ), LDB, ONE, WORK, LDWORK )
! 468: END IF
! 469: *
! 470: * col2_(4) Compute W2: = T * W2,
! 471: * T is upper-triangular.
! 472: *
! 473: CALL DTRMM( 'L', 'U', 'N', 'N', K, N-K, ONE, T, LDT,
! 474: $ WORK, LDWORK )
! 475: *
! 476: * col2_(5) Compute B2: = B2 - V2 * W2 = B2 - B1 * W2,
! 477: * V2 stored in B1.
! 478: *
! 479: IF( M.GT.0 ) THEN
! 480: CALL DGEMM( 'N', 'N', M, N-K, K, -ONE, B, LDB,
! 481: $ WORK, LDWORK, ONE, B( 1, K+1 ), LDB )
! 482: END IF
! 483: *
! 484: IF( LNOTIDENT ) THEN
! 485: *
! 486: * col2_(6) Compute W2: = V1 * W2 = A1 * W2,
! 487: * V1 is not an identity matrix, but unit lower-triangular,
! 488: * V1 stored in A1 (diagonal ones are not stored).
! 489: *
! 490: CALL DTRMM( 'L', 'L', 'N', 'U', K, N-K, ONE, A, LDA,
! 491: $ WORK, LDWORK )
! 492: END IF
! 493: *
! 494: * col2_(7) Compute A2: = A2 - W2 =
! 495: * = A(1:K, K+1:N-K) - WORK(1:K, 1:N-K),
! 496: * column-by-column.
! 497: *
! 498: DO J = 1, N-K
! 499: DO I = 1, K
! 500: A( I, K+J ) = A( I, K+J ) - WORK( I, J )
! 501: END DO
! 502: END DO
! 503: *
! 504: END IF
! 505: *
! 506: * ------------------------------------------------------------------
! 507: *
! 508: * Second Step. Computation of the Column Block 1:
! 509: *
! 510: * ( A1 ) := H * ( A1 )
! 511: * ( B1 ) ( 0 )
! 512: *
! 513: * ------------------------------------------------------------------
! 514: *
! 515: * col1_(1) Compute W1: = A1. Copy the upper-triangular
! 516: * A1 = A(1:K, 1:K) into the upper-triangular
! 517: * W1 = WORK(1:K, 1:K) column-by-column.
! 518: *
! 519: DO J = 1, K
! 520: CALL DCOPY( J, A( 1, J ), 1, WORK( 1, J ), 1 )
! 521: END DO
! 522: *
! 523: * Set the subdiagonal elements of W1 to zero column-by-column.
! 524: *
! 525: DO J = 1, K - 1
! 526: DO I = J + 1, K
! 527: WORK( I, J ) = ZERO
! 528: END DO
! 529: END DO
! 530: *
! 531: IF( LNOTIDENT ) THEN
! 532: *
! 533: * col1_(2) Compute W1: = (V1**T) * W1 = (A1**T) * W1,
! 534: * V1 is not an identity matrix, but unit lower-triangular
! 535: * V1 stored in A1 (diagonal ones are not stored),
! 536: * W1 is upper-triangular with zeroes below the diagonal.
! 537: *
! 538: CALL DTRMM( 'L', 'L', 'T', 'U', K, K, ONE, A, LDA,
! 539: $ WORK, LDWORK )
! 540: END IF
! 541: *
! 542: * col1_(3) Compute W1: = T * W1,
! 543: * T is upper-triangular,
! 544: * W1 is upper-triangular with zeroes below the diagonal.
! 545: *
! 546: CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, T, LDT,
! 547: $ WORK, LDWORK )
! 548: *
! 549: * col1_(4) Compute B1: = - V2 * W1 = - B1 * W1,
! 550: * V2 = B1, W1 is upper-triangular with zeroes below the diagonal.
! 551: *
! 552: IF( M.GT.0 ) THEN
! 553: CALL DTRMM( 'R', 'U', 'N', 'N', M, K, -ONE, WORK, LDWORK,
! 554: $ B, LDB )
! 555: END IF
! 556: *
! 557: IF( LNOTIDENT ) THEN
! 558: *
! 559: * col1_(5) Compute W1: = V1 * W1 = A1 * W1,
! 560: * V1 is not an identity matrix, but unit lower-triangular
! 561: * V1 stored in A1 (diagonal ones are not stored),
! 562: * W1 is upper-triangular on input with zeroes below the diagonal,
! 563: * and square on output.
! 564: *
! 565: CALL DTRMM( 'L', 'L', 'N', 'U', K, K, ONE, A, LDA,
! 566: $ WORK, LDWORK )
! 567: *
! 568: * col1_(6) Compute A1: = A1 - W1 = A(1:K, 1:K) - WORK(1:K, 1:K)
! 569: * column-by-column. A1 is upper-triangular on input.
! 570: * If IDENT, A1 is square on output, and W1 is square,
! 571: * if NOT IDENT, A1 is upper-triangular on output,
! 572: * W1 is upper-triangular.
! 573: *
! 574: * col1_(6)_a Compute elements of A1 below the diagonal.
! 575: *
! 576: DO J = 1, K - 1
! 577: DO I = J + 1, K
! 578: A( I, J ) = - WORK( I, J )
! 579: END DO
! 580: END DO
! 581: *
! 582: END IF
! 583: *
! 584: * col1_(6)_b Compute elements of A1 on and above the diagonal.
! 585: *
! 586: DO J = 1, K
! 587: DO I = 1, J
! 588: A( I, J ) = A( I, J ) - WORK( I, J )
! 589: END DO
! 590: END DO
! 591: *
! 592: RETURN
! 593: *
! 594: * End of DLARFB_GETT
! 595: *
! 596: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarfb_gett.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack