![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b STPRFB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download STPRFB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stprfb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stprfb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stprfb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE STPRFB( 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 * REAL A( LDA, * ), B( LDB, * ), T( LDT, * ), 00030 * $ V( LDV, * ), WORK( LDWORK, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> STPRFB applies a real "triangular-pentagonal" block reflector H or its 00040 *> conjugate transpose H^H to a real 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 realOTHERauxiliary 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 STPRFB( 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 REAL A( LDA, * ), B( LDB, * ), T( LDT, * ), 00265 $ V( LDV, * ), WORK( LDWORK, * ) 00266 * .. 00267 * 00268 * ========================================================================== 00269 * 00270 * .. Parameters .. 00271 REAL ONE, ZERO 00272 PARAMETER ( ONE = 1.0, ZERO = 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 SGEMM, STRMM 00284 * .. 00285 * .. Executable Statements .. 00286 * 00287 * Quick return if possible 00288 * 00289 IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN 00290 * 00291 IF( LSAME( STOREV, 'C' ) ) THEN 00292 COLUMN = .TRUE. 00293 ROW = .FALSE. 00294 ELSE IF ( LSAME( STOREV, 'R' ) ) THEN 00295 COLUMN = .FALSE. 00296 ROW = .TRUE. 00297 ELSE 00298 COLUMN = .FALSE. 00299 ROW = .FALSE. 00300 END IF 00301 * 00302 IF( LSAME( SIDE, 'L' ) ) THEN 00303 LEFT = .TRUE. 00304 RIGHT = .FALSE. 00305 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00306 LEFT = .FALSE. 00307 RIGHT = .TRUE. 00308 ELSE 00309 LEFT = .FALSE. 00310 RIGHT = .FALSE. 00311 END IF 00312 * 00313 IF( LSAME( DIRECT, 'F' ) ) THEN 00314 FORWARD = .TRUE. 00315 BACKWARD = .FALSE. 00316 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 00317 FORWARD = .FALSE. 00318 BACKWARD = .TRUE. 00319 ELSE 00320 FORWARD = .FALSE. 00321 BACKWARD = .FALSE. 00322 END IF 00323 * 00324 * --------------------------------------------------------------------------- 00325 * 00326 IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN 00327 * 00328 * --------------------------------------------------------------------------- 00329 * 00330 * Let W = [ I ] (K-by-K) 00331 * [ V ] (M-by-K) 00332 * 00333 * Form H C or H^H C where C = [ A ] (K-by-N) 00334 * [ B ] (M-by-N) 00335 * 00336 * H = I - W T W^H or H^H = I - W T^H W^H 00337 * 00338 * A = A - T (A + V^H B) or A = A - T^H (A + V^H B) 00339 * B = B - V T (A + V^H B) or B = B - V T^H (A + V^H B) 00340 * 00341 * --------------------------------------------------------------------------- 00342 * 00343 MP = MIN( M-L+1, M ) 00344 KP = MIN( L+1, K ) 00345 * 00346 DO J = 1, N 00347 DO I = 1, L 00348 WORK( I, J ) = B( M-L+I, J ) 00349 END DO 00350 END DO 00351 CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV, 00352 $ WORK, LDWORK ) 00353 CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, 00354 $ ONE, WORK, LDWORK ) 00355 CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, 00356 $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) 00357 * 00358 DO J = 1, N 00359 DO I = 1, K 00360 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00361 END DO 00362 END DO 00363 * 00364 CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, 00365 $ WORK, LDWORK ) 00366 * 00367 DO J = 1, N 00368 DO I = 1, K 00369 A( I, J ) = A( I, J ) - WORK( I, J ) 00370 END DO 00371 END DO 00372 * 00373 CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, 00374 $ ONE, B, LDB ) 00375 CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, 00376 $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) 00377 CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, 00378 $ WORK, LDWORK ) 00379 DO J = 1, N 00380 DO I = 1, L 00381 B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) 00382 END DO 00383 END DO 00384 * 00385 * --------------------------------------------------------------------------- 00386 * 00387 ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN 00388 * 00389 * --------------------------------------------------------------------------- 00390 * 00391 * Let W = [ I ] (K-by-K) 00392 * [ V ] (N-by-K) 00393 * 00394 * Form C H or C H^H where C = [ A B ] (A is M-by-K, B is M-by-N) 00395 * 00396 * H = I - W T W^H or H^H = I - W T^H W^H 00397 * 00398 * A = A - (A + B V) T or A = A - (A + B V) T^H 00399 * B = B - (A + B V) T V^H or B = B - (A + B V) T^H V^H 00400 * 00401 * --------------------------------------------------------------------------- 00402 * 00403 NP = MIN( N-L+1, N ) 00404 KP = MIN( L+1, K ) 00405 * 00406 DO J = 1, L 00407 DO I = 1, M 00408 WORK( I, J ) = B( I, N-L+J ) 00409 END DO 00410 END DO 00411 CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, 00412 $ WORK, LDWORK ) 00413 CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, 00414 $ V, LDV, ONE, WORK, LDWORK ) 00415 CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, 00416 $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) 00417 * 00418 DO J = 1, K 00419 DO I = 1, M 00420 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00421 END DO 00422 END DO 00423 * 00424 CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, 00425 $ WORK, LDWORK ) 00426 * 00427 DO J = 1, K 00428 DO I = 1, M 00429 A( I, J ) = A( I, J ) - WORK( I, J ) 00430 END DO 00431 END DO 00432 * 00433 CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, 00434 $ V, LDV, ONE, B, LDB ) 00435 CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, 00436 $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB ) 00437 CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV, 00438 $ WORK, LDWORK ) 00439 DO J = 1, L 00440 DO I = 1, M 00441 B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) 00442 END DO 00443 END DO 00444 * 00445 * --------------------------------------------------------------------------- 00446 * 00447 ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN 00448 * 00449 * --------------------------------------------------------------------------- 00450 * 00451 * Let W = [ V ] (M-by-K) 00452 * [ I ] (K-by-K) 00453 * 00454 * Form H C or H^H C where C = [ B ] (M-by-N) 00455 * [ A ] (K-by-N) 00456 * 00457 * H = I - W T W^H or H^H = I - W T^H W^H 00458 * 00459 * A = A - T (A + V^H B) or A = A - T^H (A + V^H B) 00460 * B = B - V T (A + V^H B) or B = B - V T^H (A + V^H B) 00461 * 00462 * --------------------------------------------------------------------------- 00463 * 00464 MP = MIN( L+1, M ) 00465 KP = MIN( K-L+1, K ) 00466 * 00467 DO J = 1, N 00468 DO I = 1, L 00469 WORK( K-L+I, J ) = B( I, J ) 00470 END DO 00471 END DO 00472 * 00473 CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV, 00474 $ WORK( KP, 1 ), LDWORK ) 00475 CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, 00476 $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) 00477 CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV, 00478 $ B, LDB, ZERO, WORK, LDWORK ) 00479 * 00480 DO J = 1, N 00481 DO I = 1, K 00482 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00483 END DO 00484 END DO 00485 * 00486 CALL STRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, 00487 $ WORK, LDWORK ) 00488 * 00489 DO J = 1, N 00490 DO I = 1, K 00491 A( I, J ) = A( I, J ) - WORK( I, J ) 00492 END DO 00493 END DO 00494 * 00495 CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, 00496 $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) 00497 CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, 00498 $ WORK, LDWORK, ONE, B, LDB ) 00499 CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV, 00500 $ WORK( KP, 1 ), LDWORK ) 00501 DO J = 1, N 00502 DO I = 1, L 00503 B( I, J ) = B( I, J ) - WORK( K-L+I, J ) 00504 END DO 00505 END DO 00506 * 00507 * --------------------------------------------------------------------------- 00508 * 00509 ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN 00510 * 00511 * --------------------------------------------------------------------------- 00512 * 00513 * Let W = [ V ] (N-by-K) 00514 * [ I ] (K-by-K) 00515 * 00516 * Form C H or C H^H where C = [ B A ] (B is M-by-N, A is M-by-K) 00517 * 00518 * H = I - W T W^H or H^H = I - W T^H W^H 00519 * 00520 * A = A - (A + B V) T or A = A - (A + B V) T^H 00521 * B = B - (A + B V) T V^H or B = B - (A + B V) T^H V^H 00522 * 00523 * --------------------------------------------------------------------------- 00524 * 00525 NP = MIN( L+1, N ) 00526 KP = MIN( K-L+1, K ) 00527 * 00528 DO J = 1, L 00529 DO I = 1, M 00530 WORK( I, K-L+J ) = B( I, J ) 00531 END DO 00532 END DO 00533 CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, 00534 $ WORK( 1, KP ), LDWORK ) 00535 CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, 00536 $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) 00537 CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, 00538 $ V, LDV, ZERO, WORK, LDWORK ) 00539 * 00540 DO J = 1, K 00541 DO I = 1, M 00542 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00543 END DO 00544 END DO 00545 * 00546 CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, 00547 $ WORK, LDWORK ) 00548 * 00549 DO J = 1, K 00550 DO I = 1, M 00551 A( I, J ) = A( I, J ) - WORK( I, J ) 00552 END DO 00553 END DO 00554 * 00555 CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, 00556 $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB ) 00557 CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK, 00558 $ V, LDV, ONE, B, LDB ) 00559 CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV, 00560 $ WORK( 1, KP ), LDWORK ) 00561 DO J = 1, L 00562 DO I = 1, M 00563 B( I, J ) = B( I, J ) - WORK( I, K-L+J ) 00564 END DO 00565 END DO 00566 * 00567 * --------------------------------------------------------------------------- 00568 * 00569 ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN 00570 * 00571 * --------------------------------------------------------------------------- 00572 * 00573 * Let W = [ I V ] ( I is K-by-K, V is K-by-M ) 00574 * 00575 * Form H C or H^H C where C = [ A ] (K-by-N) 00576 * [ B ] (M-by-N) 00577 * 00578 * H = I - W^H T W or H^H = I - W^H T^H W 00579 * 00580 * A = A - T (A + V B) or A = A - T^H (A + V B) 00581 * B = B - V^H T (A + V B) or B = B - V^H T^H (A + V B) 00582 * 00583 * --------------------------------------------------------------------------- 00584 * 00585 MP = MIN( M-L+1, M ) 00586 KP = MIN( L+1, K ) 00587 * 00588 DO J = 1, N 00589 DO I = 1, L 00590 WORK( I, J ) = B( M-L+I, J ) 00591 END DO 00592 END DO 00593 CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, 00594 $ WORK, LDB ) 00595 CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, 00596 $ ONE, WORK, LDWORK ) 00597 CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, 00598 $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) 00599 * 00600 DO J = 1, N 00601 DO I = 1, K 00602 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00603 END DO 00604 END DO 00605 * 00606 CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, 00607 $ WORK, LDWORK ) 00608 * 00609 DO J = 1, N 00610 DO I = 1, K 00611 A( I, J ) = A( I, J ) - WORK( I, J ) 00612 END DO 00613 END DO 00614 * 00615 CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, 00616 $ ONE, B, LDB ) 00617 CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, 00618 $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) 00619 CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV, 00620 $ WORK, LDWORK ) 00621 DO J = 1, N 00622 DO I = 1, L 00623 B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) 00624 END DO 00625 END DO 00626 * 00627 * --------------------------------------------------------------------------- 00628 * 00629 ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN 00630 * 00631 * --------------------------------------------------------------------------- 00632 * 00633 * Let W = [ I V ] ( I is K-by-K, V is K-by-N ) 00634 * 00635 * Form C H or C H^H where C = [ A B ] (A is M-by-K, B is M-by-N) 00636 * 00637 * H = I - W^H T W or H^H = I - W^H T^H W 00638 * 00639 * A = A - (A + B V^H) T or A = A - (A + B V^H) T^H 00640 * B = B - (A + B V^H) T V or B = B - (A + B V^H) T^H V 00641 * 00642 * --------------------------------------------------------------------------- 00643 * 00644 NP = MIN( N-L+1, N ) 00645 KP = MIN( L+1, K ) 00646 * 00647 DO J = 1, L 00648 DO I = 1, M 00649 WORK( I, J ) = B( I, N-L+J ) 00650 END DO 00651 END DO 00652 CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV, 00653 $ WORK, LDWORK ) 00654 CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV, 00655 $ ONE, WORK, LDWORK ) 00656 CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, 00657 $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) 00658 * 00659 DO J = 1, K 00660 DO I = 1, M 00661 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00662 END DO 00663 END DO 00664 * 00665 CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, 00666 $ WORK, LDWORK ) 00667 * 00668 DO J = 1, K 00669 DO I = 1, M 00670 A( I, J ) = A( I, J ) - WORK( I, J ) 00671 END DO 00672 END DO 00673 * 00674 CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 00675 $ V, LDV, ONE, B, LDB ) 00676 CALL SGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, 00677 $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) 00678 CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, 00679 $ WORK, LDWORK ) 00680 DO J = 1, L 00681 DO I = 1, M 00682 B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) 00683 END DO 00684 END DO 00685 * 00686 * --------------------------------------------------------------------------- 00687 * 00688 ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN 00689 * 00690 * --------------------------------------------------------------------------- 00691 * 00692 * Let W = [ V I ] ( I is K-by-K, V is K-by-M ) 00693 * 00694 * Form H C or H^H C where C = [ B ] (M-by-N) 00695 * [ A ] (K-by-N) 00696 * 00697 * H = I - W^H T W or H^H = I - W^H T^H W 00698 * 00699 * A = A - T (A + V B) or A = A - T^H (A + V B) 00700 * B = B - V^H T (A + V B) or B = B - V^H T^H (A + V B) 00701 * 00702 * --------------------------------------------------------------------------- 00703 * 00704 MP = MIN( L+1, M ) 00705 KP = MIN( K-L+1, K ) 00706 * 00707 DO J = 1, N 00708 DO I = 1, L 00709 WORK( K-L+I, J ) = B( I, J ) 00710 END DO 00711 END DO 00712 CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV, 00713 $ WORK( KP, 1 ), LDWORK ) 00714 CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV, 00715 $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) 00716 CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB, 00717 $ ZERO, WORK, LDWORK ) 00718 * 00719 DO J = 1, N 00720 DO I = 1, K 00721 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00722 END DO 00723 END DO 00724 * 00725 CALL STRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT, 00726 $ WORK, LDWORK ) 00727 * 00728 DO J = 1, N 00729 DO I = 1, K 00730 A( I, J ) = A( I, J ) - WORK( I, J ) 00731 END DO 00732 END DO 00733 * 00734 CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, 00735 $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) 00736 CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, 00737 $ WORK, LDWORK, ONE, B, LDB ) 00738 CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV, 00739 $ WORK( KP, 1 ), LDWORK ) 00740 DO J = 1, N 00741 DO I = 1, L 00742 B( I, J ) = B( I, J ) - WORK( K-L+I, J ) 00743 END DO 00744 END DO 00745 * 00746 * --------------------------------------------------------------------------- 00747 * 00748 ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN 00749 * 00750 * --------------------------------------------------------------------------- 00751 * 00752 * Let W = [ V I ] ( I is K-by-K, V is K-by-N ) 00753 * 00754 * Form C H or C H^H where C = [ B A ] (A is M-by-K, B is M-by-N) 00755 * 00756 * H = I - W^H T W or H^H = I - W^H T^H W 00757 * 00758 * A = A - (A + B V^H) T or A = A - (A + B V^H) T^H 00759 * B = B - (A + B V^H) T V or B = B - (A + B V^H) T^H V 00760 * 00761 * --------------------------------------------------------------------------- 00762 * 00763 NP = MIN( L+1, N ) 00764 KP = MIN( K-L+1, K ) 00765 * 00766 DO J = 1, L 00767 DO I = 1, M 00768 WORK( I, K-L+J ) = B( I, J ) 00769 END DO 00770 END DO 00771 CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV, 00772 $ WORK( 1, KP ), LDWORK ) 00773 CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB, 00774 $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) 00775 CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV, 00776 $ ZERO, WORK, LDWORK ) 00777 * 00778 DO J = 1, K 00779 DO I = 1, M 00780 WORK( I, J ) = WORK( I, J ) + A( I, J ) 00781 END DO 00782 END DO 00783 * 00784 CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, 00785 $ WORK, LDWORK ) 00786 * 00787 DO J = 1, K 00788 DO I = 1, M 00789 A( I, J ) = A( I, J ) - WORK( I, J ) 00790 END DO 00791 END DO 00792 * 00793 CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, 00794 $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) 00795 CALL SGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, 00796 $ V, LDV, ONE, B, LDB ) 00797 CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, 00798 $ WORK( 1, KP ), LDWORK ) 00799 DO J = 1, L 00800 DO I = 1, M 00801 B( I, J ) = B( I, J ) - WORK( I, K-L+J ) 00802 END DO 00803 END DO 00804 * 00805 END IF 00806 * 00807 RETURN 00808 * 00809 * End of STPRFB 00810 * 00811 END