SRC/slarzb.f

   1 *> \brief \b SLARZB applies a block reflector or its 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 SLARZB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarzb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarzb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarzb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLARZB( 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 *       REAL               C( LDC, * ), T( LDT, * ), V( LDV, * ),
  30 *      $                   WORK( LDWORK, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLARZB applies a real block reflector H or its transpose H**T to
  40 *> a real 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**T from the Left
  52 *>          = 'R': apply H or H**T 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**T (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 REAL 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 REAL 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 REAL array, dimension (LDC,N)
 136 *>          On entry, the M-by-N matrix C.
 137 *>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
 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 REAL 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 realOTHERcomputational
 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 SLARZB( 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       REAL               C( LDC, * ), T( LDT, * ), V( LDV, * ),
 197      $                   WORK( LDWORK, * )
 198 *     ..
 199 *
 200 *  =====================================================================
 201 *
 202 *     .. Parameters ..
 203       REAL               ONE
 204       PARAMETER          ( ONE = 1.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           SCOPY, SGEMM, STRMM, 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( 'SLARZB', -INFO )
 234          RETURN
 235       END IF
 236 *
 237       IF( LSAME( TRANS, 'N' ) ) THEN
 238          TRANST = 'T'
 239       ELSE
 240          TRANST = 'N'
 241       END IF
 242 *
 243       IF( LSAME( SIDE, 'L' ) ) THEN
 244 *
 245 *        Form  H * C  or  H**T * C
 246 *
 247 *        W( 1:n, 1:k ) = C( 1:k, 1:n )**T
 248 *
 249          DO 10 J = 1, K
 250             CALL SCOPY( 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 )**T * V( 1:k, 1:l )**T
 255 *
 256          IF( L.GT.0 )
 257      $      CALL SGEMM( 'Transpose', 'Transpose', N, K, L, ONE,
 258      $                  C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, LDWORK )
 259 *
 260 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
 261 *
 262          CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
 263      $               LDT, WORK, LDWORK )
 264 *
 265 *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**T
 266 *
 267          DO 30 J = 1, N
 268             DO 20 I = 1, K
 269                C( I, J ) = C( I, J ) - WORK( J, I )
 270    20       CONTINUE
 271    30    CONTINUE
 272 *
 273 *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
 274 *                            V( 1:k, 1:l )**T * W( 1:n, 1:k )**T
 275 *
 276          IF( L.GT.0 )
 277      $      CALL SGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
 278      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
 279 *
 280       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 281 *
 282 *        Form  C * H  or  C * H**T
 283 *
 284 *        W( 1:m, 1:k ) = C( 1:m, 1:k )
 285 *
 286          DO 40 J = 1, K
 287             CALL SCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 288    40    CONTINUE
 289 *
 290 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
 291 *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**T
 292 *
 293          IF( L.GT.0 )
 294      $      CALL SGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
 295      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
 296 *
 297 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T**T
 298 *
 299          CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
 300      $               LDT, WORK, LDWORK )
 301 *
 302 *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
 303 *
 304          DO 60 J = 1, K
 305             DO 50 I = 1, M
 306                C( I, J ) = C( I, J ) - WORK( I, J )
 307    50       CONTINUE
 308    60    CONTINUE
 309 *
 310 *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
 311 *                            W( 1:m, 1:k ) * V( 1:k, 1:l )
 312 *
 313          IF( L.GT.0 )
 314      $      CALL SGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
 315      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
 316 *
 317       END IF
 318 *
 319       RETURN
 320 *
 321 *     End of SLARZB
 322 *
 323       END