LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlarzb.f
Go to the documentation of this file.
00001 *> \brief \b ZLARZB
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLARZB + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarzb.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarzb.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarzb.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
00022 *                          LDV, T, LDT, C, LDC, WORK, LDWORK )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          DIRECT, SIDE, STOREV, TRANS
00026 *       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
00030 *      $                   WORK( LDWORK, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> ZLARZB applies a complex block reflector H or its transpose H**H
00040 *> to a complex distributed M-by-N  C from the left or the right.
00041 *>
00042 *> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
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, not supported yet)
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': Columnwise                        (not supported yet)
00077 *>          = 'R': Rowwise
00078 *> \endverbatim
00079 *>
00080 *> \param[in] M
00081 *> \verbatim
00082 *>          M is INTEGER
00083 *>          The number of rows of the matrix C.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] N
00087 *> \verbatim
00088 *>          N is INTEGER
00089 *>          The number of columns of the matrix C.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] K
00093 *> \verbatim
00094 *>          K is INTEGER
00095 *>          The order of the matrix T (= the number of elementary
00096 *>          reflectors whose product defines the block reflector).
00097 *> \endverbatim
00098 *>
00099 *> \param[in] L
00100 *> \verbatim
00101 *>          L is INTEGER
00102 *>          The number of columns of the matrix V containing the
00103 *>          meaningful part of the Householder reflectors.
00104 *>          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] V
00108 *> \verbatim
00109 *>          V is COMPLEX*16 array, dimension (LDV,NV).
00110 *>          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
00111 *> \endverbatim
00112 *>
00113 *> \param[in] LDV
00114 *> \verbatim
00115 *>          LDV is INTEGER
00116 *>          The leading dimension of the array V.
00117 *>          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
00118 *> \endverbatim
00119 *>
00120 *> \param[in] T
00121 *> \verbatim
00122 *>          T is COMPLEX*16 array, dimension (LDT,K)
00123 *>          The triangular K-by-K matrix T in the representation of the
00124 *>          block reflector.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] LDT
00128 *> \verbatim
00129 *>          LDT is INTEGER
00130 *>          The leading dimension of the array T. LDT >= K.
00131 *> \endverbatim
00132 *>
00133 *> \param[in,out] C
00134 *> \verbatim
00135 *>          C is COMPLEX*16 array, dimension (LDC,N)
00136 *>          On entry, the M-by-N matrix C.
00137 *>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
00138 *> \endverbatim
00139 *>
00140 *> \param[in] LDC
00141 *> \verbatim
00142 *>          LDC is INTEGER
00143 *>          The leading dimension of the array C. LDC >= max(1,M).
00144 *> \endverbatim
00145 *>
00146 *> \param[out] WORK
00147 *> \verbatim
00148 *>          WORK is COMPLEX*16 array, dimension (LDWORK,K)
00149 *> \endverbatim
00150 *>
00151 *> \param[in] LDWORK
00152 *> \verbatim
00153 *>          LDWORK is INTEGER
00154 *>          The leading dimension of the array WORK.
00155 *>          If SIDE = 'L', LDWORK >= max(1,N);
00156 *>          if SIDE = 'R', LDWORK >= max(1,M).
00157 *> \endverbatim
00158 *
00159 *  Authors:
00160 *  ========
00161 *
00162 *> \author Univ. of Tennessee 
00163 *> \author Univ. of California Berkeley 
00164 *> \author Univ. of Colorado Denver 
00165 *> \author NAG Ltd. 
00166 *
00167 *> \date November 2011
00168 *
00169 *> \ingroup complex16OTHERcomputational
00170 *
00171 *> \par Contributors:
00172 *  ==================
00173 *>
00174 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00175 *
00176 *> \par Further Details:
00177 *  =====================
00178 *>
00179 *> \verbatim
00180 *> \endverbatim
00181 *>
00182 *  =====================================================================
00183       SUBROUTINE ZLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
00184      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
00185 *
00186 *  -- LAPACK computational routine (version 3.4.0) --
00187 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00188 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00189 *     November 2011
00190 *
00191 *     .. Scalar Arguments ..
00192       CHARACTER          DIRECT, SIDE, STOREV, TRANS
00193       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
00194 *     ..
00195 *     .. Array Arguments ..
00196       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
00197      $                   WORK( LDWORK, * )
00198 *     ..
00199 *
00200 *  =====================================================================
00201 *
00202 *     .. Parameters ..
00203       COMPLEX*16         ONE
00204       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
00205 *     ..
00206 *     .. Local Scalars ..
00207       CHARACTER          TRANST
00208       INTEGER            I, INFO, J
00209 *     ..
00210 *     .. External Functions ..
00211       LOGICAL            LSAME
00212       EXTERNAL           LSAME
00213 *     ..
00214 *     .. External Subroutines ..
00215       EXTERNAL           XERBLA, ZCOPY, ZGEMM, ZLACGV, ZTRMM
00216 *     ..
00217 *     .. Executable Statements ..
00218 *
00219 *     Quick return if possible
00220 *
00221       IF( M.LE.0 .OR. N.LE.0 )
00222      $   RETURN
00223 *
00224 *     Check for currently supported options
00225 *
00226       INFO = 0
00227       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
00228          INFO = -3
00229       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
00230          INFO = -4
00231       END IF
00232       IF( INFO.NE.0 ) THEN
00233          CALL XERBLA( 'ZLARZB', -INFO )
00234          RETURN
00235       END IF
00236 *
00237       IF( LSAME( TRANS, 'N' ) ) THEN
00238          TRANST = 'C'
00239       ELSE
00240          TRANST = 'N'
00241       END IF
00242 *
00243       IF( LSAME( SIDE, 'L' ) ) THEN
00244 *
00245 *        Form  H * C  or  H**H * C
00246 *
00247 *        W( 1:n, 1:k ) = C( 1:k, 1:n )**H
00248 *
00249          DO 10 J = 1, K
00250             CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
00251    10    CONTINUE
00252 *
00253 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
00254 *                        C( m-l+1:m, 1:n )**H * V( 1:k, 1:l )**T
00255 *
00256          IF( L.GT.0 )
00257      $      CALL ZGEMM( 'Transpose', 'Conjugate transpose', N, K, L,
00258      $                  ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK,
00259      $                  LDWORK )
00260 *
00261 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
00262 *
00263          CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
00264      $               LDT, WORK, LDWORK )
00265 *
00266 *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**H
00267 *
00268          DO 30 J = 1, N
00269             DO 20 I = 1, K
00270                C( I, J ) = C( I, J ) - WORK( J, I )
00271    20       CONTINUE
00272    30    CONTINUE
00273 *
00274 *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
00275 *                            V( 1:k, 1:l )**H * W( 1:n, 1:k )**H
00276 *
00277          IF( L.GT.0 )
00278      $      CALL ZGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
00279      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
00280 *
00281       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
00282 *
00283 *        Form  C * H  or  C * H**H
00284 *
00285 *        W( 1:m, 1:k ) = C( 1:m, 1:k )
00286 *
00287          DO 40 J = 1, K
00288             CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
00289    40    CONTINUE
00290 *
00291 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
00292 *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**H
00293 *
00294          IF( L.GT.0 )
00295      $      CALL ZGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
00296      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
00297 *
00298 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T )  or
00299 *                        W( 1:m, 1:k ) * T**H
00300 *
00301          DO 50 J = 1, K
00302             CALL ZLACGV( K-J+1, T( J, J ), 1 )
00303    50    CONTINUE
00304          CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
00305      $               LDT, WORK, LDWORK )
00306          DO 60 J = 1, K
00307             CALL ZLACGV( K-J+1, T( J, J ), 1 )
00308    60    CONTINUE
00309 *
00310 *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
00311 *
00312          DO 80 J = 1, K
00313             DO 70 I = 1, M
00314                C( I, J ) = C( I, J ) - WORK( I, J )
00315    70       CONTINUE
00316    80    CONTINUE
00317 *
00318 *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
00319 *                            W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) )
00320 *
00321          DO 90 J = 1, L
00322             CALL ZLACGV( K, V( 1, J ), 1 )
00323    90    CONTINUE
00324          IF( L.GT.0 )
00325      $      CALL ZGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
00326      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
00327          DO 100 J = 1, L
00328             CALL ZLACGV( K, V( 1, J ), 1 )
00329   100    CONTINUE
00330 *
00331       END IF
00332 *
00333       RETURN
00334 *
00335 *     End of ZLARZB
00336 *
00337       END
 All Files Functions