SRC/clarzb.f

   1 *> \brief \b CLARZB applies a block reflector or its conjugate-transpose to a general matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CLARZB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarzb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarzb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarzb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
  22 *                          LDV, T, LDT, C, LDC, WORK, LDWORK )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          DIRECT, SIDE, STOREV, TRANS
  26 *       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX            C( LDC, * ), T( LDT, * ), V( LDV, * ),
  30 *      $                   WORK( LDWORK, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> CLARZB applies a complex block reflector H or its transpose H**H
  40 *> to a complex distributed M-by-N  C from the left or the right.
  41 *>
  42 *> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] SIDE
  49 *> \verbatim
  50 *>          SIDE is CHARACTER*1
  51 *>          = 'L': apply H or H**H from the Left
  52 *>          = 'R': apply H or H**H from the Right
  53 *> \endverbatim
  54 *>
  55 *> \param[in] TRANS
  56 *> \verbatim
  57 *>          TRANS is CHARACTER*1
  58 *>          = 'N': apply H (No transpose)
  59 *>          = 'C': apply H**H (Conjugate transpose)
  60 *> \endverbatim
  61 *>
  62 *> \param[in] DIRECT
  63 *> \verbatim
  64 *>          DIRECT is CHARACTER*1
  65 *>          Indicates how H is formed from a product of elementary
  66 *>          reflectors
  67 *>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
  68 *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
  69 *> \endverbatim
  70 *>
  71 *> \param[in] STOREV
  72 *> \verbatim
  73 *>          STOREV is CHARACTER*1
  74 *>          Indicates how the vectors which define the elementary
  75 *>          reflectors are stored:
  76 *>          = 'C': Columnwise                        (not supported yet)
  77 *>          = 'R': Rowwise
  78 *> \endverbatim
  79 *>
  80 *> \param[in] M
  81 *> \verbatim
  82 *>          M is INTEGER
  83 *>          The number of rows of the matrix C.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] N
  87 *> \verbatim
  88 *>          N is INTEGER
  89 *>          The number of columns of the matrix C.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] K
  93 *> \verbatim
  94 *>          K is INTEGER
  95 *>          The order of the matrix T (= the number of elementary
  96 *>          reflectors whose product defines the block reflector).
  97 *> \endverbatim
  98 *>
  99 *> \param[in] L
 100 *> \verbatim
 101 *>          L is INTEGER
 102 *>          The number of columns of the matrix V containing the
 103 *>          meaningful part of the Householder reflectors.
 104 *>          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] V
 108 *> \verbatim
 109 *>          V is COMPLEX array, dimension (LDV,NV).
 110 *>          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] LDV
 114 *> \verbatim
 115 *>          LDV is INTEGER
 116 *>          The leading dimension of the array V.
 117 *>          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] T
 121 *> \verbatim
 122 *>          T is COMPLEX array, dimension (LDT,K)
 123 *>          The triangular K-by-K matrix T in the representation of the
 124 *>          block reflector.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] LDT
 128 *> \verbatim
 129 *>          LDT is INTEGER
 130 *>          The leading dimension of the array T. LDT >= K.
 131 *> \endverbatim
 132 *>
 133 *> \param[in,out] C
 134 *> \verbatim
 135 *>          C is COMPLEX array, dimension (LDC,N)
 136 *>          On entry, the M-by-N matrix C.
 137 *>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] LDC
 141 *> \verbatim
 142 *>          LDC is INTEGER
 143 *>          The leading dimension of the array C. LDC >= max(1,M).
 144 *> \endverbatim
 145 *>
 146 *> \param[out] WORK
 147 *> \verbatim
 148 *>          WORK is COMPLEX array, dimension (LDWORK,K)
 149 *> \endverbatim
 150 *>
 151 *> \param[in] LDWORK
 152 *> \verbatim
 153 *>          LDWORK is INTEGER
 154 *>          The leading dimension of the array WORK.
 155 *>          If SIDE = 'L', LDWORK >= max(1,N);
 156 *>          if SIDE = 'R', LDWORK >= max(1,M).
 157 *> \endverbatim
 158 *
 159 *  Authors:
 160 *  ========
 161 *
 162 *> \author Univ. of Tennessee
 163 *> \author Univ. of California Berkeley
 164 *> \author Univ. of Colorado Denver
 165 *> \author NAG Ltd.
 166 *
 167 *> \date September 2012
 168 *
 169 *> \ingroup complexOTHERcomputational
 170 *
 171 *> \par Contributors:
 172 *  ==================
 173 *>
 174 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 175 *
 176 *> \par Further Details:
 177 *  =====================
 178 *>
 179 *> \verbatim
 180 *> \endverbatim
 181 *>
 182 *  =====================================================================
 183       SUBROUTINE CLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
 184      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
 185 *
 186 *  -- LAPACK computational routine (version 3.4.2) --
 187 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 188 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 189 *     September 2012
 190 *
 191 *     .. Scalar Arguments ..
 192       CHARACTER          DIRECT, SIDE, STOREV, TRANS
 193       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
 194 *     ..
 195 *     .. Array Arguments ..
 196       COMPLEX            C( LDC, * ), T( LDT, * ), V( LDV, * ),
 197      $                   WORK( LDWORK, * )
 198 *     ..
 199 *
 200 *  =====================================================================
 201 *
 202 *     .. Parameters ..
 203       COMPLEX            ONE
 204       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
 205 *     ..
 206 *     .. Local Scalars ..
 207       CHARACTER          TRANST
 208       INTEGER            I, INFO, J
 209 *     ..
 210 *     .. External Functions ..
 211       LOGICAL            LSAME
 212       EXTERNAL           LSAME
 213 *     ..
 214 *     .. External Subroutines ..
 215       EXTERNAL           CCOPY, CGEMM, CLACGV, CTRMM, XERBLA
 216 *     ..
 217 *     .. Executable Statements ..
 218 *
 219 *     Quick return if possible
 220 *
 221       IF( M.LE.0 .OR. N.LE.0 )
 222      $   RETURN
 223 *
 224 *     Check for currently supported options
 225 *
 226       INFO = 0
 227       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
 228          INFO = -3
 229       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
 230          INFO = -4
 231       END IF
 232       IF( INFO.NE.0 ) THEN
 233          CALL XERBLA( 'CLARZB', -INFO )
 234          RETURN
 235       END IF
 236 *
 237       IF( LSAME( TRANS, 'N' ) ) THEN
 238          TRANST = 'C'
 239       ELSE
 240          TRANST = 'N'
 241       END IF
 242 *
 243       IF( LSAME( SIDE, 'L' ) ) THEN
 244 *
 245 *        Form  H * C  or  H**H * C
 246 *
 247 *        W( 1:n, 1:k ) = C( 1:k, 1:n )**H
 248 *
 249          DO 10 J = 1, K
 250             CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 251    10    CONTINUE
 252 *
 253 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
 254 *                        C( m-l+1:m, 1:n )**H * V( 1:k, 1:l )**T
 255 *
 256          IF( L.GT.0 )
 257      $      CALL CGEMM( 'Transpose', 'Conjugate transpose', N, K, L,
 258      $                  ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK,
 259      $                  LDWORK )
 260 *
 261 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
 262 *
 263          CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
 264      $               LDT, WORK, LDWORK )
 265 *
 266 *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**H
 267 *
 268          DO 30 J = 1, N
 269             DO 20 I = 1, K
 270                C( I, J ) = C( I, J ) - WORK( J, I )
 271    20       CONTINUE
 272    30    CONTINUE
 273 *
 274 *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
 275 *                            V( 1:k, 1:l )**H * W( 1:n, 1:k )**H
 276 *
 277          IF( L.GT.0 )
 278      $      CALL CGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
 279      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
 280 *
 281       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 282 *
 283 *        Form  C * H  or  C * H**H
 284 *
 285 *        W( 1:m, 1:k ) = C( 1:m, 1:k )
 286 *
 287          DO 40 J = 1, K
 288             CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 289    40    CONTINUE
 290 *
 291 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
 292 *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**H
 293 *
 294          IF( L.GT.0 )
 295      $      CALL CGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
 296      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
 297 *
 298 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T )  or
 299 *                        W( 1:m, 1:k ) * T**H
 300 *
 301          DO 50 J = 1, K
 302             CALL CLACGV( K-J+1, T( J, J ), 1 )
 303    50    CONTINUE
 304          CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
 305      $               LDT, WORK, LDWORK )
 306          DO 60 J = 1, K
 307             CALL CLACGV( K-J+1, T( J, J ), 1 )
 308    60    CONTINUE
 309 *
 310 *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
 311 *
 312          DO 80 J = 1, K
 313             DO 70 I = 1, M
 314                C( I, J ) = C( I, J ) - WORK( I, J )
 315    70       CONTINUE
 316    80    CONTINUE
 317 *
 318 *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
 319 *                            W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) )
 320 *
 321          DO 90 J = 1, L
 322             CALL CLACGV( K, V( 1, J ), 1 )
 323    90    CONTINUE
 324          IF( L.GT.0 )
 325      $      CALL CGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
 326      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
 327          DO 100 J = 1, L
 328             CALL CLACGV( K, V( 1, J ), 1 )
 329   100    CONTINUE
 330 *
 331       END IF
 332 *
 333       RETURN
 334 *
 335 *     End of CLARZB
 336 *
 337       END