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