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