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