LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stprfb.f
Go to the documentation of this file.
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
 All Files Functions