![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLARFB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLARFB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLARFB( 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 * COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), 00030 * $ WORK( LDWORK, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> CLARFB applies a complex block reflector H or its transpose H**H to a 00040 *> complex 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**H from the Left 00050 *> = 'R': apply H or H**H from the Right 00051 *> \endverbatim 00052 *> 00053 *> \param[in] TRANS 00054 *> \verbatim 00055 *> TRANS is CHARACTER*1 00056 *> = 'N': apply H (No transpose) 00057 *> = 'C': apply H**H (Conjugate 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 COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDC,N) 00131 *> On entry, the M-by-N matrix C. 00132 *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H. 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 COMPLEX 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 complexOTHERauxiliary 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 CLARFB( 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 COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), 00209 $ WORK( LDWORK, * ) 00210 * .. 00211 * 00212 * ===================================================================== 00213 * 00214 * .. Parameters .. 00215 COMPLEX ONE 00216 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 00217 * .. 00218 * .. Local Scalars .. 00219 CHARACTER TRANST 00220 INTEGER I, J, LASTV, LASTC 00221 * .. 00222 * .. External Functions .. 00223 LOGICAL LSAME 00224 INTEGER ILACLR, ILACLC 00225 EXTERNAL LSAME, ILACLR, ILACLC 00226 * .. 00227 * .. External Subroutines .. 00228 EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM 00229 * .. 00230 * .. Intrinsic Functions .. 00231 INTRINSIC CONJG 00232 * .. 00233 * .. Executable Statements .. 00234 * 00235 * Quick return if possible 00236 * 00237 IF( M.LE.0 .OR. N.LE.0 ) 00238 $ RETURN 00239 * 00240 IF( LSAME( TRANS, 'N' ) ) THEN 00241 TRANST = 'C' 00242 ELSE 00243 TRANST = 'N' 00244 END IF 00245 * 00246 IF( LSAME( STOREV, 'C' ) ) THEN 00247 * 00248 IF( LSAME( DIRECT, 'F' ) ) THEN 00249 * 00250 * Let V = ( V1 ) (first K rows) 00251 * ( V2 ) 00252 * where V1 is unit lower triangular. 00253 * 00254 IF( LSAME( SIDE, 'L' ) ) THEN 00255 * 00256 * Form H * C or H**H * C where C = ( C1 ) 00257 * ( C2 ) 00258 * 00259 LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) 00260 LASTC = ILACLC( LASTV, N, C, LDC ) 00261 * 00262 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) 00263 * 00264 * W := C1**H 00265 * 00266 DO 10 J = 1, K 00267 CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 00268 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 00269 10 CONTINUE 00270 * 00271 * W := W * V1 00272 * 00273 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00274 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00275 IF( LASTV.GT.K ) THEN 00276 * 00277 * W := W + C2**H *V2 00278 * 00279 CALL CGEMM( 'Conjugate transpose', 'No transpose', 00280 $ LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC, 00281 $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK ) 00282 END IF 00283 * 00284 * W := W * T**H or W * T 00285 * 00286 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 00287 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00288 * 00289 * C := C - V * W**H 00290 * 00291 IF( M.GT.K ) THEN 00292 * 00293 * C2 := C2 - V2 * W**H 00294 * 00295 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00296 $ LASTV-K, LASTC, K, -ONE, V( K+1, 1 ), LDV, 00297 $ WORK, LDWORK, ONE, C( K+1, 1 ), LDC ) 00298 END IF 00299 * 00300 * W := W * V1**H 00301 * 00302 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 00303 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00304 * 00305 * C1 := C1 - W**H 00306 * 00307 DO 30 J = 1, K 00308 DO 20 I = 1, LASTC 00309 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) 00310 20 CONTINUE 00311 30 CONTINUE 00312 * 00313 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00314 * 00315 * Form C * H or C * H**H where C = ( C1 C2 ) 00316 * 00317 LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) 00318 LASTC = ILACLR( M, LASTV, C, LDC ) 00319 * 00320 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00321 * 00322 * W := C1 00323 * 00324 DO 40 J = 1, K 00325 CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 00326 40 CONTINUE 00327 * 00328 * W := W * V1 00329 * 00330 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00331 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00332 IF( LASTV.GT.K ) THEN 00333 * 00334 * W := W + C2 * V2 00335 * 00336 CALL CGEMM( 'No transpose', 'No transpose', 00337 $ LASTC, K, LASTV-K, 00338 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 00339 $ ONE, WORK, LDWORK ) 00340 END IF 00341 * 00342 * W := W * T or W * T**H 00343 * 00344 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 00345 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00346 * 00347 * C := C - W * V**H 00348 * 00349 IF( LASTV.GT.K ) THEN 00350 * 00351 * C2 := C2 - W * V2**H 00352 * 00353 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00354 $ LASTC, LASTV-K, K, 00355 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, 00356 $ ONE, C( 1, K+1 ), LDC ) 00357 END IF 00358 * 00359 * W := W * V1**H 00360 * 00361 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 00362 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00363 * 00364 * C1 := C1 - W 00365 * 00366 DO 60 J = 1, K 00367 DO 50 I = 1, LASTC 00368 C( I, J ) = C( I, J ) - WORK( I, J ) 00369 50 CONTINUE 00370 60 CONTINUE 00371 END IF 00372 * 00373 ELSE 00374 * 00375 * Let V = ( V1 ) 00376 * ( V2 ) (last K rows) 00377 * where V2 is unit upper triangular. 00378 * 00379 IF( LSAME( SIDE, 'L' ) ) THEN 00380 * 00381 * Form H * C or H**H * C where C = ( C1 ) 00382 * ( C2 ) 00383 * 00384 LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) 00385 LASTC = ILACLC( LASTV, N, C, LDC ) 00386 * 00387 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) 00388 * 00389 * W := C2**H 00390 * 00391 DO 70 J = 1, K 00392 CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 00393 $ WORK( 1, J ), 1 ) 00394 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 00395 70 CONTINUE 00396 * 00397 * W := W * V2 00398 * 00399 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00400 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00401 $ WORK, LDWORK ) 00402 IF( LASTV.GT.K ) THEN 00403 * 00404 * W := W + C1**H*V1 00405 * 00406 CALL CGEMM( 'Conjugate transpose', 'No transpose', 00407 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, 00408 $ ONE, WORK, LDWORK ) 00409 END IF 00410 * 00411 * W := W * T**H or W * T 00412 * 00413 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 00414 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00415 * 00416 * C := C - V * W**H 00417 * 00418 IF( LASTV.GT.K ) THEN 00419 * 00420 * C1 := C1 - V1 * W**H 00421 * 00422 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00423 $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, 00424 $ ONE, C, LDC ) 00425 END IF 00426 * 00427 * W := W * V2**H 00428 * 00429 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 00430 $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00431 $ WORK, LDWORK ) 00432 * 00433 * C2 := C2 - W**H 00434 * 00435 DO 90 J = 1, K 00436 DO 80 I = 1, LASTC 00437 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - 00438 $ CONJG( WORK( I, J ) ) 00439 80 CONTINUE 00440 90 CONTINUE 00441 * 00442 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00443 * 00444 * Form C * H or C * H**H where C = ( C1 C2 ) 00445 * 00446 LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) 00447 LASTC = ILACLR( M, LASTV, C, LDC ) 00448 * 00449 * W := C * V = (C1*V1 + C2*V2) (stored in WORK) 00450 * 00451 * W := C2 00452 * 00453 DO 100 J = 1, K 00454 CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, 00455 $ WORK( 1, J ), 1 ) 00456 100 CONTINUE 00457 * 00458 * W := W * V2 00459 * 00460 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00461 $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00462 $ WORK, LDWORK ) 00463 IF( LASTV.GT.K ) THEN 00464 * 00465 * W := W + C1 * V1 00466 * 00467 CALL CGEMM( 'No transpose', 'No transpose', 00468 $ LASTC, K, LASTV-K, 00469 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 00470 END IF 00471 * 00472 * W := W * T or W * T**H 00473 * 00474 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 00475 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00476 * 00477 * C := C - W * V**H 00478 * 00479 IF( LASTV.GT.K ) THEN 00480 * 00481 * C1 := C1 - W * V1**H 00482 * 00483 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00484 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 00485 $ ONE, C, LDC ) 00486 END IF 00487 * 00488 * W := W * V2**H 00489 * 00490 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 00491 $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, 00492 $ WORK, LDWORK ) 00493 * 00494 * C2 := C2 - W 00495 * 00496 DO 120 J = 1, K 00497 DO 110 I = 1, LASTC 00498 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) 00499 $ - WORK( I, J ) 00500 110 CONTINUE 00501 120 CONTINUE 00502 END IF 00503 END IF 00504 * 00505 ELSE IF( LSAME( STOREV, 'R' ) ) THEN 00506 * 00507 IF( LSAME( DIRECT, 'F' ) ) THEN 00508 * 00509 * Let V = ( V1 V2 ) (V1: first K columns) 00510 * where V1 is unit upper triangular. 00511 * 00512 IF( LSAME( SIDE, 'L' ) ) THEN 00513 * 00514 * Form H * C or H**H * C where C = ( C1 ) 00515 * ( C2 ) 00516 * 00517 LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) 00518 LASTC = ILACLC( LASTV, N, C, LDC ) 00519 * 00520 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) 00521 * 00522 * W := C1**H 00523 * 00524 DO 130 J = 1, K 00525 CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 00526 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 00527 130 CONTINUE 00528 * 00529 * W := W * V1**H 00530 * 00531 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 00532 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00533 IF( LASTV.GT.K ) THEN 00534 * 00535 * W := W + C2**H*V2**H 00536 * 00537 CALL CGEMM( 'Conjugate transpose', 00538 $ 'Conjugate transpose', LASTC, K, LASTV-K, 00539 $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, 00540 $ ONE, WORK, LDWORK ) 00541 END IF 00542 * 00543 * W := W * T**H or W * T 00544 * 00545 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', 00546 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00547 * 00548 * C := C - V**H * W**H 00549 * 00550 IF( LASTV.GT.K ) THEN 00551 * 00552 * C2 := C2 - V2**H * W**H 00553 * 00554 CALL CGEMM( 'Conjugate transpose', 00555 $ 'Conjugate transpose', LASTV-K, LASTC, K, 00556 $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, 00557 $ ONE, C( K+1, 1 ), LDC ) 00558 END IF 00559 * 00560 * W := W * V1 00561 * 00562 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00563 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00564 * 00565 * C1 := C1 - W**H 00566 * 00567 DO 150 J = 1, K 00568 DO 140 I = 1, LASTC 00569 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) 00570 140 CONTINUE 00571 150 CONTINUE 00572 * 00573 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00574 * 00575 * Form C * H or C * H**H where C = ( C1 C2 ) 00576 * 00577 LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) 00578 LASTC = ILACLR( M, LASTV, C, LDC ) 00579 * 00580 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) 00581 * 00582 * W := C1 00583 * 00584 DO 160 J = 1, K 00585 CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) 00586 160 CONTINUE 00587 * 00588 * W := W * V1**H 00589 * 00590 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', 00591 $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00592 IF( LASTV.GT.K ) THEN 00593 * 00594 * W := W + C2 * V2**H 00595 * 00596 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00597 $ LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC, 00598 $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK ) 00599 END IF 00600 * 00601 * W := W * T or W * T**H 00602 * 00603 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', 00604 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00605 * 00606 * C := C - W * V 00607 * 00608 IF( LASTV.GT.K ) THEN 00609 * 00610 * C2 := C2 - W * V2 00611 * 00612 CALL CGEMM( 'No transpose', 'No transpose', 00613 $ LASTC, LASTV-K, K, 00614 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, 00615 $ ONE, C( 1, K+1 ), LDC ) 00616 END IF 00617 * 00618 * W := W * V1 00619 * 00620 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', 00621 $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) 00622 * 00623 * C1 := C1 - W 00624 * 00625 DO 180 J = 1, K 00626 DO 170 I = 1, LASTC 00627 C( I, J ) = C( I, J ) - WORK( I, J ) 00628 170 CONTINUE 00629 180 CONTINUE 00630 * 00631 END IF 00632 * 00633 ELSE 00634 * 00635 * Let V = ( V1 V2 ) (V2: last K columns) 00636 * where V2 is unit lower triangular. 00637 * 00638 IF( LSAME( SIDE, 'L' ) ) THEN 00639 * 00640 * Form H * C or H**H * C where C = ( C1 ) 00641 * ( C2 ) 00642 * 00643 LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) 00644 LASTC = ILACLC( LASTV, N, C, LDC ) 00645 * 00646 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) 00647 * 00648 * W := C2**H 00649 * 00650 DO 190 J = 1, K 00651 CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, 00652 $ WORK( 1, J ), 1 ) 00653 CALL CLACGV( LASTC, WORK( 1, J ), 1 ) 00654 190 CONTINUE 00655 * 00656 * W := W * V2**H 00657 * 00658 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 00659 $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00660 $ WORK, LDWORK ) 00661 IF( LASTV.GT.K ) THEN 00662 * 00663 * W := W + C1**H * V1**H 00664 * 00665 CALL CGEMM( 'Conjugate transpose', 00666 $ 'Conjugate transpose', LASTC, K, LASTV-K, 00667 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 00668 END IF 00669 * 00670 * W := W * T**H or W * T 00671 * 00672 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', 00673 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00674 * 00675 * C := C - V**H * W**H 00676 * 00677 IF( LASTV.GT.K ) THEN 00678 * 00679 * C1 := C1 - V1**H * W**H 00680 * 00681 CALL CGEMM( 'Conjugate transpose', 00682 $ 'Conjugate transpose', LASTV-K, LASTC, K, 00683 $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC ) 00684 END IF 00685 * 00686 * W := W * V2 00687 * 00688 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00689 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00690 $ WORK, LDWORK ) 00691 * 00692 * C2 := C2 - W**H 00693 * 00694 DO 210 J = 1, K 00695 DO 200 I = 1, LASTC 00696 C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - 00697 $ CONJG( WORK( I, J ) ) 00698 200 CONTINUE 00699 210 CONTINUE 00700 * 00701 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00702 * 00703 * Form C * H or C * H**H where C = ( C1 C2 ) 00704 * 00705 LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) 00706 LASTC = ILACLR( M, LASTV, C, LDC ) 00707 * 00708 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) 00709 * 00710 * W := C2 00711 * 00712 DO 220 J = 1, K 00713 CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, 00714 $ WORK( 1, J ), 1 ) 00715 220 CONTINUE 00716 * 00717 * W := W * V2**H 00718 * 00719 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', 00720 $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00721 $ WORK, LDWORK ) 00722 IF( LASTV.GT.K ) THEN 00723 * 00724 * W := W + C1 * V1**H 00725 * 00726 CALL CGEMM( 'No transpose', 'Conjugate transpose', 00727 $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE, 00728 $ WORK, LDWORK ) 00729 END IF 00730 * 00731 * W := W * T or W * T**H 00732 * 00733 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', 00734 $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) 00735 * 00736 * C := C - W * V 00737 * 00738 IF( LASTV.GT.K ) THEN 00739 * 00740 * C1 := C1 - W * V1 00741 * 00742 CALL CGEMM( 'No transpose', 'No transpose', 00743 $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, 00744 $ ONE, C, LDC ) 00745 END IF 00746 * 00747 * W := W * V2 00748 * 00749 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', 00750 $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, 00751 $ WORK, LDWORK ) 00752 * 00753 * C1 := C1 - W 00754 * 00755 DO 240 J = 1, K 00756 DO 230 I = 1, LASTC 00757 C( I, LASTV-K+J ) = C( I, LASTV-K+J ) 00758 $ - WORK( I, J ) 00759 230 CONTINUE 00760 240 CONTINUE 00761 * 00762 END IF 00763 * 00764 END IF 00765 END IF 00766 * 00767 RETURN 00768 * 00769 * End of CLARFB 00770 * 00771 END