 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, 00002 $ T, LDT, C, LDC, WORK, LDWORK ) 00003 IMPLICIT NONE 00004 * 00005 * -- LAPACK auxiliary routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER DIRECT, SIDE, STOREV, TRANS 00012 INTEGER K, LDC, LDT, LDV, LDWORK, M, N 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), 00016 $ WORK( LDWORK, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DLARFB applies a real block reflector H or its transpose H**T to a 00023 * real m by n matrix C, from either the left or the right. 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * SIDE (input) CHARACTER*1 00029 * = 'L': apply H or H**T from the Left 00030 * = 'R': apply H or H**T from the Right 00031 * 00032 * TRANS (input) CHARACTER*1 00033 * = 'N': apply H (No transpose) 00034 * = 'T': apply H**T (Transpose) 00035 * 00036 * DIRECT (input) CHARACTER*1 00037 * Indicates how H is formed from a product of elementary 00038 * reflectors 00039 * = 'F': H = H(1) H(2) . . . H(k) (Forward) 00040 * = 'B': H = H(k) . . . H(2) H(1) (Backward) 00041 * 00042 * STOREV (input) CHARACTER*1 00043 * Indicates how the vectors which define the elementary 00044 * reflectors are stored: 00045 * = 'C': Columnwise 00046 * = 'R': Rowwise 00047 * 00048 * M (input) INTEGER 00049 * The number of rows of the matrix C. 00050 * 00051 * N (input) INTEGER 00052 * The number of columns of the matrix C. 00053 * 00054 * K (input) INTEGER 00055 * The order of the matrix T (= the number of elementary 00056 * reflectors whose product defines the block reflector). 00057 * 00058 * V (input) DOUBLE PRECISION array, dimension 00059 * (LDV,K) if STOREV = 'C' 00060 * (LDV,M) if STOREV = 'R' and SIDE = 'L' 00061 * (LDV,N) if STOREV = 'R' and SIDE = 'R' 00062 * The matrix V. See Further Details. 00063 * 00064 * LDV (input) INTEGER 00065 * The leading dimension of the array V. 00066 * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); 00067 * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); 00068 * if STOREV = 'R', LDV >= K. 00069 * 00070 * T (input) DOUBLE PRECISION array, dimension (LDT,K) 00071 * The triangular k by k matrix T in the representation of the 00072 * block reflector. 00073 * 00074 * LDT (input) INTEGER 00075 * The leading dimension of the array T. LDT >= K. 00076 * 00077 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00078 * On entry, the m by n matrix C. 00079 * On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T. 00080 * 00081 * LDC (input) INTEGER 00082 * The leading dimension of the array C. LDC >= max(1,M). 00083 * 00084 * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) 00085 * 00086 * LDWORK (input) INTEGER 00087 * The leading dimension of the array WORK. 00088 * If SIDE = 'L', LDWORK >= max(1,N); 00089 * if SIDE = 'R', LDWORK >= max(1,M). 00090 * 00091 * Further Details 00092 * =============== 00093 * 00094 * The shape of the matrix V and the storage of the vectors which define 00095 * the H(i) is best illustrated by the following example with n = 5 and 00096 * k = 3. The elements equal to 1 are not stored; the corresponding 00097 * array elements are modified but restored on exit. The rest of the 00098 * array is not used. 00099 * 00100 * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 00101 * 00102 * V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 00103 * ( v1 1 ) ( 1 v2 v2 v2 ) 00104 * ( v1 v2 1 ) ( 1 v3 v3 ) 00105 * ( v1 v2 v3 ) 00106 * ( v1 v2 v3 ) 00107 * 00108 * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 00109 * 00110 * V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 00111 * ( v1 v2 v3 ) ( v2 v2 v2 1 ) 00112 * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 00113 * ( 1 v3 ) 00114 * ( 1 ) 00115 * 00116 * ===================================================================== 00117 * 00118 * .. Parameters .. 00119 DOUBLE PRECISION ONE 00120 PARAMETER ( ONE = 1.0D+0 ) 00121 * .. 00122 * .. Local Scalars .. 00123 CHARACTER TRANST 00124 INTEGER I, J, LASTV, LASTC 00125 * .. 00126 * .. External Functions .. 00127 LOGICAL LSAME 00128 INTEGER ILADLR, ILADLC 00129 EXTERNAL LSAME, ILADLR, ILADLC 00130 * .. 00131 * .. External Subroutines .. 00132 EXTERNAL DCOPY, DGEMM, DTRMM 00133 * .. 00134 * .. Executable Statements .. 00135 * 00136 * Quick return if possible 00137 * 00138 IF( M.LE.0 .OR. N.LE.0 ) 00139 $ RETURN 00140 * 00141 IF( LSAME( TRANS, 'N' ) ) THEN 00142 TRANST = 'T' 00143 ELSE 00144 TRANST = 'N' 00145 END IF 00146 * 00147 IF( LSAME( STOREV, 'C' ) ) THEN 00148 * 00149 IF( LSAME( DIRECT, 'F' ) ) THEN 00150 * 00151 * Let V = ( V1 ) (first K rows) 00152 * ( V2 ) 00153 * where V1 is unit lower triangular. 00154 * 00155 IF( LSAME( SIDE, 'L' ) ) THEN 00156 * 00157 * Form H * C or H**T * C where C = ( C1 ) 00158 * ( C2 ) 00159 * 00160 LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 00161 LASTC = ILADLC( LASTV, N, C, LDC ) 00162 * 00163 * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) 00164 * 00165 * W := C1**T 00166 * 00167 DO 10 J = 1, K 00168 CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 00169 10 CONTINUE 00170 * 00171 * W := W * V1 00172 * 00173 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00174 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00175 IF( LASTV.GT.K ) THEN 00176 * 00177 * W := W + C2**T *V2 00178 * 00179 CALL DGEMM( 'Transpose', 'No transpose', 00180 $ LASTC, K, LASTV-K, 00181 $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, 00182 $ ONE, WORK, LDWORK ) 00183 END IF 00184 * 00185 * W := W * T**T or W * T 00186 * 00187 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 00188 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00189 * 00190 * C := C - V * W**T 00191 * 00192 IF( LASTV.GT.K ) THEN 00193 * 00194 * C2 := C2 - V2 * W**T 00195 * 00196 CALL DGEMM( 'No transpose', 'Transpose', 00197 $ LASTV-K, LASTC, K, 00198 $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, 00199 $ C( K+1, 1 ), LDC ) 00200 END IF 00201 * 00202 * W := W * V1**T 00203 * 00204 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 00205 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00206 * 00207 * C1 := C1 - W**T 00208 * 00209 DO 30 J = 1, K 00210 DO 20 I = 1, LASTC 00211 C( J, I ) = C( J, I ) - WORK( I, J ) 00212 20 CONTINUE 00213 30 CONTINUE 00214 * 00215 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00216 * 00217 * Form C * H or C * H**T where C = ( C1 C2 ) 00218 * 00219 LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 00220 LASTC = ILADLR( M, LASTV, C, LDC ) 00221 * 00222 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00223 * 00224 * W := C1 00225 * 00226 DO 40 J = 1, K 00227 CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 00228 40 CONTINUE 00229 * 00230 * W := W * V1 00231 * 00232 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00233 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00234 IF( LASTV.GT.K ) THEN 00235 * 00236 * W := W + C2 * V2 00237 * 00238 CALL DGEMM( 'No transpose', 'No transpose', 00239 $ LASTC, K, LASTV-K, 00240 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 00241 $ ONE, WORK, LDWORK ) 00242 END IF 00243 * 00244 * W := W * T or W * T**T 00245 * 00246 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 00247 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00248 * 00249 * C := C - W * V**T 00250 * 00251 IF( LASTV.GT.K ) THEN 00252 * 00253 * C2 := C2 - W * V2**T 00254 * 00255 CALL DGEMM( 'No transpose', 'Transpose', 00256 $ LASTC, LASTV-K, K, 00257 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, 00258 $ C( 1, K+1 ), LDC ) 00259 END IF 00260 * 00261 * W := W * V1**T 00262 * 00263 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 00264 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00265 * 00266 * C1 := C1 - W 00267 * 00268 DO 60 J = 1, K 00269 DO 50 I = 1, LASTC 00270 C( I, J ) = C( I, J ) - WORK( I, J ) 00271 50 CONTINUE 00272 60 CONTINUE 00273 END IF 00274 * 00275 ELSE 00276 * 00277 * Let V = ( V1 ) 00278 * ( V2 ) (last K rows) 00279 * where V2 is unit upper triangular. 00280 * 00281 IF( LSAME( SIDE, 'L' ) ) THEN 00282 * 00283 * Form H * C or H**T * C where C = ( C1 ) 00284 * ( C2 ) 00285 * 00286 LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) 00287 LASTC = ILADLC( LASTV, N, C, LDC ) 00288 * 00289 * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) 00290 * 00291 * W := C2**T 00292 * 00293 DO 70 J = 1, K 00294 CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 00295 $ WORK( 1, J ), 1 ) 00296 70 CONTINUE 00297 * 00298 * W := W * V2 00299 * 00300 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00301 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00302 $ WORK, LDWORK ) 00303 IF( LASTV.GT.K ) THEN 00304 * 00305 * W := W + C1**T*V1 00306 * 00307 CALL DGEMM( 'Transpose', 'No transpose', 00308 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 00309 $ ONE, WORK, LDWORK ) 00310 END IF 00311 * 00312 * W := W * T**T or W * T 00313 * 00314 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 00315 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00316 * 00317 * C := C - V * W**T 00318 * 00319 IF( LASTV.GT.K ) THEN 00320 * 00321 * C1 := C1 - V1 * W**T 00322 * 00323 CALL DGEMM( 'No transpose', 'Transpose', 00324 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 00325 $ ONE, C, LDC ) 00326 END IF 00327 * 00328 * W := W * V2**T 00329 * 00330 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 00331 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00332 $ WORK, LDWORK ) 00333 * 00334 * C2 := C2 - W**T 00335 * 00336 DO 90 J = 1, K 00337 DO 80 I = 1, LASTC 00338 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 00339 80 CONTINUE 00340 90 CONTINUE 00341 * 00342 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00343 * 00344 * Form C * H or C * H**T where C = ( C1 C2 ) 00345 * 00346 LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) 00347 LASTC = ILADLR( M, LASTV, C, LDC ) 00348 * 00349 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00350 * 00351 * W := C2 00352 * 00353 DO 100 J = 1, K 00354 CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) 00355 100 CONTINUE 00356 * 00357 * W := W * V2 00358 * 00359 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00360 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00361 $ WORK, LDWORK ) 00362 IF( LASTV.GT.K ) THEN 00363 * 00364 * W := W + C1 * V1 00365 * 00366 CALL DGEMM( 'No transpose', 'No transpose', 00367 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 00368 $ ONE, WORK, LDWORK ) 00369 END IF 00370 * 00371 * W := W * T or W * T**T 00372 * 00373 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 00374 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00375 * 00376 * C := C - W * V**T 00377 * 00378 IF( LASTV.GT.K ) THEN 00379 * 00380 * C1 := C1 - W * V1**T 00381 * 00382 CALL DGEMM( 'No transpose', 'Transpose', 00383 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 00384 $ ONE, C, LDC ) 00385 END IF 00386 * 00387 * W := W * V2**T 00388 * 00389 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 00390 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00391 $ WORK, LDWORK ) 00392 * 00393 * C2 := C2 - W 00394 * 00395 DO 120 J = 1, K 00396 DO 110 I = 1, LASTC 00397 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 00398 110 CONTINUE 00399 120 CONTINUE 00400 END IF 00401 END IF 00402 * 00403 ELSE IF( LSAME( STOREV, 'R' ) ) THEN 00404 * 00405 IF( LSAME( DIRECT, 'F' ) ) THEN 00406 * 00407 * Let V = ( V1 V2 ) (V1: first K columns) 00408 * where V1 is unit upper triangular. 00409 * 00410 IF( LSAME( SIDE, 'L' ) ) THEN 00411 * 00412 * Form H * C or H**T * C where C = ( C1 ) 00413 * ( C2 ) 00414 * 00415 LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 00416 LASTC = ILADLC( LASTV, N, C, LDC ) 00417 * 00418 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) 00419 * 00420 * W := C1**T 00421 * 00422 DO 130 J = 1, K 00423 CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 00424 130 CONTINUE 00425 * 00426 * W := W * V1**T 00427 * 00428 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 00429 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00430 IF( LASTV.GT.K ) THEN 00431 * 00432 * W := W + C2**T*V2**T 00433 * 00434 CALL DGEMM( 'Transpose', 'Transpose', 00435 $ LASTC, K, LASTV-K, 00436 $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, 00437 $ ONE, WORK, LDWORK ) 00438 END IF 00439 * 00440 * W := W * T**T or W * T 00441 * 00442 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 00443 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00444 * 00445 * C := C - V**T * W**T 00446 * 00447 IF( LASTV.GT.K ) THEN 00448 * 00449 * C2 := C2 - V2**T * W**T 00450 * 00451 CALL DGEMM( 'Transpose', 'Transpose', 00452 $ LASTV-K, LASTC, K, 00453 $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, 00454 $ ONE, C( K+1, 1 ), LDC ) 00455 END IF 00456 * 00457 * W := W * V1 00458 * 00459 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00460 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00461 * 00462 * C1 := C1 - W**T 00463 * 00464 DO 150 J = 1, K 00465 DO 140 I = 1, LASTC 00466 C( J, I ) = C( J, I ) - WORK( I, J ) 00467 140 CONTINUE 00468 150 CONTINUE 00469 * 00470 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00471 * 00472 * Form C * H or C * H**T where C = ( C1 C2 ) 00473 * 00474 LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 00475 LASTC = ILADLR( M, LASTV, C, LDC ) 00476 * 00477 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) 00478 * 00479 * W := C1 00480 * 00481 DO 160 J = 1, K 00482 CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 00483 160 CONTINUE 00484 * 00485 * W := W * V1**T 00486 * 00487 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', 00488 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00489 IF( LASTV.GT.K ) THEN 00490 * 00491 * W := W + C2 * V2**T 00492 * 00493 CALL DGEMM( 'No transpose', 'Transpose', 00494 $ LASTC, K, LASTV-K, 00495 $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, 00496 $ ONE, WORK, LDWORK ) 00497 END IF 00498 * 00499 * W := W * T or W * T**T 00500 * 00501 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 00502 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00503 * 00504 * C := C - W * V 00505 * 00506 IF( LASTV.GT.K ) THEN 00507 * 00508 * C2 := C2 - W * V2 00509 * 00510 CALL DGEMM( 'No transpose', 'No transpose', 00511 $ LASTC, LASTV-K, K, 00512 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, 00513 $ ONE, C( 1, K+1 ), LDC ) 00514 END IF 00515 * 00516 * W := W * V1 00517 * 00518 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00519 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00520 * 00521 * C1 := C1 - W 00522 * 00523 DO 180 J = 1, K 00524 DO 170 I = 1, LASTC 00525 C( I, J ) = C( I, J ) - WORK( I, J ) 00526 170 CONTINUE 00527 180 CONTINUE 00528 * 00529 END IF 00530 * 00531 ELSE 00532 * 00533 * Let V = ( V1 V2 ) (V2: last K columns) 00534 * where V2 is unit lower triangular. 00535 * 00536 IF( LSAME( SIDE, 'L' ) ) THEN 00537 * 00538 * Form H * C or H**T * C where C = ( C1 ) 00539 * ( C2 ) 00540 * 00541 LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) 00542 LASTC = ILADLC( LASTV, N, C, LDC ) 00543 * 00544 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) 00545 * 00546 * W := C2**T 00547 * 00548 DO 190 J = 1, K 00549 CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 00550 $ WORK( 1, J ), 1 ) 00551 190 CONTINUE 00552 * 00553 * W := W * V2**T 00554 * 00555 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 00556 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00557 $ WORK, LDWORK ) 00558 IF( LASTV.GT.K ) THEN 00559 * 00560 * W := W + C1**T * V1**T 00561 * 00562 CALL DGEMM( 'Transpose', 'Transpose', 00563 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 00564 $ ONE, WORK, LDWORK ) 00565 END IF 00566 * 00567 * W := W * T**T or W * T 00568 * 00569 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 00570 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00571 * 00572 * C := C - V**T * W**T 00573 * 00574 IF( LASTV.GT.K ) THEN 00575 * 00576 * C1 := C1 - V1**T * W**T 00577 * 00578 CALL DGEMM( 'Transpose', 'Transpose', 00579 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 00580 $ ONE, C, LDC ) 00581 END IF 00582 * 00583 * W := W * V2 00584 * 00585 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00586 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00587 $ WORK, LDWORK ) 00588 * 00589 * C2 := C2 - W**T 00590 * 00591 DO 210 J = 1, K 00592 DO 200 I = 1, LASTC 00593 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) 00594 200 CONTINUE 00595 210 CONTINUE 00596 * 00597 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00598 * 00599 * Form C * H or C * H**T where C = ( C1 C2 ) 00600 * 00601 LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) 00602 LASTC = ILADLR( M, LASTV, C, LDC ) 00603 * 00604 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) 00605 * 00606 * W := C2 00607 * 00608 DO 220 J = 1, K 00609 CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1, 00610 $ WORK( 1, J ), 1 ) 00611 220 CONTINUE 00612 * 00613 * W := W * V2**T 00614 * 00615 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', 00616 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00617 $ WORK, LDWORK ) 00618 IF( LASTV.GT.K ) THEN 00619 * 00620 * W := W + C1 * V1**T 00621 * 00622 CALL DGEMM( 'No transpose', 'Transpose', 00623 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 00624 $ ONE, WORK, LDWORK ) 00625 END IF 00626 * 00627 * W := W * T or W * T**T 00628 * 00629 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 00630 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00631 * 00632 * C := C - W * V 00633 * 00634 IF( LASTV.GT.K ) THEN 00635 * 00636 * C1 := C1 - W * V1 00637 * 00638 CALL DGEMM( 'No transpose', 'No transpose', 00639 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 00640 $ ONE, C, LDC ) 00641 END IF 00642 * 00643 * W := W * V2 00644 * 00645 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00646 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00647 $ WORK, LDWORK ) 00648 * 00649 * C1 := C1 - W 00650 * 00651 DO 240 J = 1, K 00652 DO 230 I = 1, LASTC 00653 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) 00654 230 CONTINUE 00655 240 CONTINUE 00656 * 00657 END IF 00658 * 00659 END IF 00660 END IF 00661 * 00662 RETURN 00663 * 00664 * End of DLARFB 00665 * 00666 END
1.7.3 
