[136] | 1 | *> \brief \b SLARFB
|
---|
| 2 | *
|
---|
| 3 | * =========== DOCUMENTATION ===========
|
---|
| 4 | *
|
---|
| 5 | * Online html documentation available at
|
---|
| 6 | * http://www.netlib.org/lapack/explore-html/
|
---|
| 7 | *
|
---|
| 8 | *> \htmlonly
|
---|
| 9 | *> Download SLARFB + dependencies
|
---|
| 10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfb.f">
|
---|
| 11 | *> [TGZ]</a>
|
---|
| 12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfb.f">
|
---|
| 13 | *> [ZIP]</a>
|
---|
| 14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfb.f">
|
---|
| 15 | *> [TXT]</a>
|
---|
| 16 | *> \endhtmlonly
|
---|
| 17 | *
|
---|
| 18 | * Definition:
|
---|
| 19 | * ===========
|
---|
| 20 | *
|
---|
| 21 | * SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
|
---|
| 22 | * T, LDT, C, LDC, WORK, LDWORK )
|
---|
| 23 | *
|
---|
| 24 | * .. Scalar Arguments ..
|
---|
| 25 | * CHARACTER DIRECT, SIDE, STOREV, TRANS
|
---|
| 26 | * INTEGER K, 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 | *> SLARFB applies a real block reflector H or its transpose H**T to a
|
---|
| 40 | *> real m by n matrix C, from either the left or the right.
|
---|
| 41 | *> \endverbatim
|
---|
| 42 | *
|
---|
| 43 | * Arguments:
|
---|
| 44 | * ==========
|
---|
| 45 | *
|
---|
| 46 | *> \param[in] SIDE
|
---|
| 47 | *> \verbatim
|
---|
| 48 | *> SIDE is CHARACTER*1
|
---|
| 49 | *> = 'L': apply H or H**T from the Left
|
---|
| 50 | *> = 'R': apply H or H**T from the Right
|
---|
| 51 | *> \endverbatim
|
---|
| 52 | *>
|
---|
| 53 | *> \param[in] TRANS
|
---|
| 54 | *> \verbatim
|
---|
| 55 | *> TRANS is CHARACTER*1
|
---|
| 56 | *> = 'N': apply H (No transpose)
|
---|
| 57 | *> = 'T': apply H**T (Transpose)
|
---|
| 58 | *> \endverbatim
|
---|
| 59 | *>
|
---|
| 60 | *> \param[in] DIRECT
|
---|
| 61 | *> \verbatim
|
---|
| 62 | *> DIRECT is CHARACTER*1
|
---|
| 63 | *> Indicates how H is formed from a product of elementary
|
---|
| 64 | *> reflectors
|
---|
| 65 | *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
|
---|
| 66 | *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
|
---|
| 67 | *> \endverbatim
|
---|
| 68 | *>
|
---|
| 69 | *> \param[in] STOREV
|
---|
| 70 | *> \verbatim
|
---|
| 71 | *> STOREV is CHARACTER*1
|
---|
| 72 | *> Indicates how the vectors which define the elementary
|
---|
| 73 | *> reflectors are stored:
|
---|
| 74 | *> = 'C': Columnwise
|
---|
| 75 | *> = 'R': Rowwise
|
---|
| 76 | *> \endverbatim
|
---|
| 77 | *>
|
---|
| 78 | *> \param[in] M
|
---|
| 79 | *> \verbatim
|
---|
| 80 | *> M is INTEGER
|
---|
| 81 | *> The number of rows of the matrix C.
|
---|
| 82 | *> \endverbatim
|
---|
| 83 | *>
|
---|
| 84 | *> \param[in] N
|
---|
| 85 | *> \verbatim
|
---|
| 86 | *> N is INTEGER
|
---|
| 87 | *> The number of columns of the matrix C.
|
---|
| 88 | *> \endverbatim
|
---|
| 89 | *>
|
---|
| 90 | *> \param[in] K
|
---|
| 91 | *> \verbatim
|
---|
| 92 | *> K is INTEGER
|
---|
| 93 | *> The order of the matrix T (= the number of elementary
|
---|
| 94 | *> reflectors whose product defines the block reflector).
|
---|
| 95 | *> \endverbatim
|
---|
| 96 | *>
|
---|
| 97 | *> \param[in] V
|
---|
| 98 | *> \verbatim
|
---|
| 99 | *> V is REAL array, dimension
|
---|
| 100 | *> (LDV,K) if STOREV = 'C'
|
---|
| 101 | *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
|
---|
| 102 | *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
|
---|
| 103 | *> The matrix V. See Further Details.
|
---|
| 104 | *> \endverbatim
|
---|
| 105 | *>
|
---|
| 106 | *> \param[in] LDV
|
---|
| 107 | *> \verbatim
|
---|
| 108 | *> LDV is INTEGER
|
---|
| 109 | *> The leading dimension of the array V.
|
---|
| 110 | *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
|
---|
| 111 | *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
|
---|
| 112 | *> if STOREV = 'R', LDV >= K.
|
---|
| 113 | *> \endverbatim
|
---|
| 114 | *>
|
---|
| 115 | *> \param[in] T
|
---|
| 116 | *> \verbatim
|
---|
| 117 | *> T is REAL array, dimension (LDT,K)
|
---|
| 118 | *> The triangular k by k matrix T in the representation of the
|
---|
| 119 | *> block reflector.
|
---|
| 120 | *> \endverbatim
|
---|
| 121 | *>
|
---|
| 122 | *> \param[in] LDT
|
---|
| 123 | *> \verbatim
|
---|
| 124 | *> LDT is INTEGER
|
---|
| 125 | *> The leading dimension of the array T. LDT >= K.
|
---|
| 126 | *> \endverbatim
|
---|
| 127 | *>
|
---|
| 128 | *> \param[in,out] C
|
---|
| 129 | *> \verbatim
|
---|
| 130 | *> C is REAL array, dimension (LDC,N)
|
---|
| 131 | *> On entry, the m by n matrix C.
|
---|
| 132 | *> On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
|
---|
| 133 | *> \endverbatim
|
---|
| 134 | *>
|
---|
| 135 | *> \param[in] LDC
|
---|
| 136 | *> \verbatim
|
---|
| 137 | *> LDC is INTEGER
|
---|
| 138 | *> The leading dimension of the array C. LDC >= max(1,M).
|
---|
| 139 | *> \endverbatim
|
---|
| 140 | *>
|
---|
| 141 | *> \param[out] WORK
|
---|
| 142 | *> \verbatim
|
---|
| 143 | *> WORK is REAL array, dimension (LDWORK,K)
|
---|
| 144 | *> \endverbatim
|
---|
| 145 | *>
|
---|
| 146 | *> \param[in] LDWORK
|
---|
| 147 | *> \verbatim
|
---|
| 148 | *> LDWORK is INTEGER
|
---|
| 149 | *> The leading dimension of the array WORK.
|
---|
| 150 | *> If SIDE = 'L', LDWORK >= max(1,N);
|
---|
| 151 | *> if SIDE = 'R', LDWORK >= max(1,M).
|
---|
| 152 | *> \endverbatim
|
---|
| 153 | *
|
---|
| 154 | * Authors:
|
---|
| 155 | * ========
|
---|
| 156 | *
|
---|
| 157 | *> \author Univ. of Tennessee
|
---|
| 158 | *> \author Univ. of California Berkeley
|
---|
| 159 | *> \author Univ. of Colorado Denver
|
---|
| 160 | *> \author NAG Ltd.
|
---|
| 161 | *
|
---|
| 162 | *> \date November 2011
|
---|
| 163 | *
|
---|
| 164 | *> \ingroup realOTHERauxiliary
|
---|
| 165 | *
|
---|
| 166 | *> \par Further Details:
|
---|
| 167 | * =====================
|
---|
| 168 | *>
|
---|
| 169 | *> \verbatim
|
---|
| 170 | *>
|
---|
| 171 | *> The shape of the matrix V and the storage of the vectors which define
|
---|
| 172 | *> the H(i) is best illustrated by the following example with n = 5 and
|
---|
| 173 | *> k = 3. The elements equal to 1 are not stored; the corresponding
|
---|
| 174 | *> array elements are modified but restored on exit. The rest of the
|
---|
| 175 | *> array is not used.
|
---|
| 176 | *>
|
---|
| 177 | *> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
|
---|
| 178 | *>
|
---|
| 179 | *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
|
---|
| 180 | *> ( v1 1 ) ( 1 v2 v2 v2 )
|
---|
| 181 | *> ( v1 v2 1 ) ( 1 v3 v3 )
|
---|
| 182 | *> ( v1 v2 v3 )
|
---|
| 183 | *> ( v1 v2 v3 )
|
---|
| 184 | *>
|
---|
| 185 | *> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
|
---|
| 186 | *>
|
---|
| 187 | *> V = ( v1 v2 v3 ) V = ( v1 v1 1 )
|
---|
| 188 | *> ( v1 v2 v3 ) ( v2 v2 v2 1 )
|
---|
| 189 | *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
|
---|
| 190 | *> ( 1 v3 )
|
---|
| 191 | *> ( 1 )
|
---|
| 192 | *> \endverbatim
|
---|
| 193 | *>
|
---|
| 194 | * =====================================================================
|
---|
| 195 | SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
|
---|
| 196 | $ T, LDT, C, LDC, WORK, LDWORK )
|
---|
| 197 | *
|
---|
| 198 | * -- LAPACK auxiliary routine (version 3.4.0) --
|
---|
| 199 | * -- LAPACK is a software package provided by Univ. of Tennessee, --
|
---|
| 200 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
---|
| 201 | * November 2011
|
---|
| 202 | *
|
---|
| 203 | * .. Scalar Arguments ..
|
---|
| 204 | CHARACTER DIRECT, SIDE, STOREV, TRANS
|
---|
| 205 | INTEGER K, LDC, LDT, LDV, LDWORK, M, N
|
---|
| 206 | * ..
|
---|
| 207 | * .. Array Arguments ..
|
---|
| 208 | REAL C( LDC, * ), T( LDT, * ), V( LDV, * ),
|
---|
| 209 | $ WORK( LDWORK, * )
|
---|
| 210 | * ..
|
---|
| 211 | *
|
---|
| 212 | * =====================================================================
|
---|
| 213 | *
|
---|
| 214 | * .. Parameters ..
|
---|
| 215 | REAL ONE
|
---|
| 216 | PARAMETER ( ONE = 1.0E+0 )
|
---|
| 217 | * ..
|
---|
| 218 | * .. Local Scalars ..
|
---|
| 219 | CHARACTER TRANST
|
---|
| 220 | INTEGER I, J, LASTV, LASTC
|
---|
| 221 | * ..
|
---|
| 222 | * .. External Functions ..
|
---|
| 223 | LOGICAL LSAME
|
---|
| 224 | INTEGER ILASLR, ILASLC
|
---|
| 225 | EXTERNAL LSAME, ILASLR, ILASLC
|
---|
| 226 | * ..
|
---|
| 227 | * .. External Subroutines ..
|
---|
| 228 | EXTERNAL SCOPY, SGEMM, STRMM
|
---|
| 229 | * ..
|
---|
| 230 | * .. Executable Statements ..
|
---|
| 231 | *
|
---|
| 232 | * Quick return if possible
|
---|
| 233 | *
|
---|
| 234 | IF( M.LE.0 .OR. N.LE.0 )
|
---|
| 235 | $ RETURN
|
---|
| 236 | *
|
---|
| 237 | IF( LSAME( TRANS, 'N' ) ) THEN
|
---|
| 238 | TRANST = 'T'
|
---|
| 239 | ELSE
|
---|
| 240 | TRANST = 'N'
|
---|
| 241 | END IF
|
---|
| 242 | *
|
---|
| 243 | IF( LSAME( STOREV, 'C' ) ) THEN
|
---|
| 244 | *
|
---|
| 245 | IF( LSAME( DIRECT, 'F' ) ) THEN
|
---|
| 246 | *
|
---|
| 247 | * Let V = ( V1 ) (first K rows)
|
---|
| 248 | * ( V2 )
|
---|
| 249 | * where V1 is unit lower triangular.
|
---|
| 250 | *
|
---|
| 251 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 252 | *
|
---|
| 253 | * Form H * C or H**T * C where C = ( C1 )
|
---|
| 254 | * ( C2 )
|
---|
| 255 | *
|
---|
| 256 | LASTV = MAX( K, ILASLR( M, K, V, LDV ) )
|
---|
| 257 | LASTC = ILASLC( LASTV, N, C, LDC )
|
---|
| 258 | *
|
---|
| 259 | * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK)
|
---|
| 260 | *
|
---|
| 261 | * W := C1**T
|
---|
| 262 | *
|
---|
| 263 | DO 10 J = 1, K
|
---|
| 264 | CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
|
---|
| 265 | 10 CONTINUE
|
---|
| 266 | *
|
---|
| 267 | * W := W * V1
|
---|
| 268 | *
|
---|
| 269 | CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 270 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 271 | IF( LASTV.GT.K ) THEN
|
---|
| 272 | *
|
---|
| 273 | * W := W + C2**T *V2
|
---|
| 274 | *
|
---|
| 275 | CALL SGEMM( 'Transpose', 'No transpose',
|
---|
| 276 | $ LASTC, K, LASTV-K,
|
---|
| 277 | $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
|
---|
| 278 | $ ONE, WORK, LDWORK )
|
---|
| 279 | END IF
|
---|
| 280 | *
|
---|
| 281 | * W := W * T**T or W * T
|
---|
| 282 | *
|
---|
| 283 | CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit',
|
---|
| 284 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 285 | *
|
---|
| 286 | * C := C - V * W**T
|
---|
| 287 | *
|
---|
| 288 | IF( LASTV.GT.K ) THEN
|
---|
| 289 | *
|
---|
| 290 | * C2 := C2 - V2 * W**T
|
---|
| 291 | *
|
---|
| 292 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 293 | $ LASTV-K, LASTC, K,
|
---|
| 294 | $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
|
---|
| 295 | $ C( K+1, 1 ), LDC )
|
---|
| 296 | END IF
|
---|
| 297 | *
|
---|
| 298 | * W := W * V1**T
|
---|
| 299 | *
|
---|
| 300 | CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
|
---|
| 301 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 302 | *
|
---|
| 303 | * C1 := C1 - W**T
|
---|
| 304 | *
|
---|
| 305 | DO 30 J = 1, K
|
---|
| 306 | DO 20 I = 1, LASTC
|
---|
| 307 | C( J, I ) = C( J, I ) - WORK( I, J )
|
---|
| 308 | 20 CONTINUE
|
---|
| 309 | 30 CONTINUE
|
---|
| 310 | *
|
---|
| 311 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 312 | *
|
---|
| 313 | * Form C * H or C * H**T where C = ( C1 C2 )
|
---|
| 314 | *
|
---|
| 315 | LASTV = MAX( K, ILASLR( N, K, V, LDV ) )
|
---|
| 316 | LASTC = ILASLR( M, LASTV, C, LDC )
|
---|
| 317 | *
|
---|
| 318 | * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
|
---|
| 319 | *
|
---|
| 320 | * W := C1
|
---|
| 321 | *
|
---|
| 322 | DO 40 J = 1, K
|
---|
| 323 | CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
|
---|
| 324 | 40 CONTINUE
|
---|
| 325 | *
|
---|
| 326 | * W := W * V1
|
---|
| 327 | *
|
---|
| 328 | CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 329 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 330 | IF( LASTV.GT.K ) THEN
|
---|
| 331 | *
|
---|
| 332 | * W := W + C2 * V2
|
---|
| 333 | *
|
---|
| 334 | CALL SGEMM( 'No transpose', 'No transpose',
|
---|
| 335 | $ LASTC, K, LASTV-K,
|
---|
| 336 | $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
|
---|
| 337 | $ ONE, WORK, LDWORK )
|
---|
| 338 | END IF
|
---|
| 339 | *
|
---|
| 340 | * W := W * T or W * T**T
|
---|
| 341 | *
|
---|
| 342 | CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit',
|
---|
| 343 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 344 | *
|
---|
| 345 | * C := C - W * V**T
|
---|
| 346 | *
|
---|
| 347 | IF( LASTV.GT.K ) THEN
|
---|
| 348 | *
|
---|
| 349 | * C2 := C2 - W * V2**T
|
---|
| 350 | *
|
---|
| 351 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 352 | $ LASTC, LASTV-K, K,
|
---|
| 353 | $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
|
---|
| 354 | $ C( 1, K+1 ), LDC )
|
---|
| 355 | END IF
|
---|
| 356 | *
|
---|
| 357 | * W := W * V1**T
|
---|
| 358 | *
|
---|
| 359 | CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
|
---|
| 360 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 361 | *
|
---|
| 362 | * C1 := C1 - W
|
---|
| 363 | *
|
---|
| 364 | DO 60 J = 1, K
|
---|
| 365 | DO 50 I = 1, LASTC
|
---|
| 366 | C( I, J ) = C( I, J ) - WORK( I, J )
|
---|
| 367 | 50 CONTINUE
|
---|
| 368 | 60 CONTINUE
|
---|
| 369 | END IF
|
---|
| 370 | *
|
---|
| 371 | ELSE
|
---|
| 372 | *
|
---|
| 373 | * Let V = ( V1 )
|
---|
| 374 | * ( V2 ) (last K rows)
|
---|
| 375 | * where V2 is unit upper triangular.
|
---|
| 376 | *
|
---|
| 377 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 378 | *
|
---|
| 379 | * Form H * C or H**T * C where C = ( C1 )
|
---|
| 380 | * ( C2 )
|
---|
| 381 | *
|
---|
| 382 | LASTV = MAX( K, ILASLR( M, K, V, LDV ) )
|
---|
| 383 | LASTC = ILASLC( LASTV, N, C, LDC )
|
---|
| 384 | *
|
---|
| 385 | * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK)
|
---|
| 386 | *
|
---|
| 387 | * W := C2**T
|
---|
| 388 | *
|
---|
| 389 | DO 70 J = 1, K
|
---|
| 390 | CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
|
---|
| 391 | $ WORK( 1, J ), 1 )
|
---|
| 392 | 70 CONTINUE
|
---|
| 393 | *
|
---|
| 394 | * W := W * V2
|
---|
| 395 | *
|
---|
| 396 | CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 397 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 398 | $ WORK, LDWORK )
|
---|
| 399 | IF( LASTV.GT.K ) THEN
|
---|
| 400 | *
|
---|
| 401 | * W := W + C1**T*V1
|
---|
| 402 | *
|
---|
| 403 | CALL SGEMM( 'Transpose', 'No transpose',
|
---|
| 404 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
|
---|
| 405 | $ ONE, WORK, LDWORK )
|
---|
| 406 | END IF
|
---|
| 407 | *
|
---|
| 408 | * W := W * T**T or W * T
|
---|
| 409 | *
|
---|
| 410 | CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit',
|
---|
| 411 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 412 | *
|
---|
| 413 | * C := C - V * W**T
|
---|
| 414 | *
|
---|
| 415 | IF( LASTV.GT.K ) THEN
|
---|
| 416 | *
|
---|
| 417 | * C1 := C1 - V1 * W**T
|
---|
| 418 | *
|
---|
| 419 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 420 | $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
|
---|
| 421 | $ ONE, C, LDC )
|
---|
| 422 | END IF
|
---|
| 423 | *
|
---|
| 424 | * W := W * V2**T
|
---|
| 425 | *
|
---|
| 426 | CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
|
---|
| 427 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 428 | $ WORK, LDWORK )
|
---|
| 429 | *
|
---|
| 430 | * C2 := C2 - W**T
|
---|
| 431 | *
|
---|
| 432 | DO 90 J = 1, K
|
---|
| 433 | DO 80 I = 1, LASTC
|
---|
| 434 | C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
|
---|
| 435 | 80 CONTINUE
|
---|
| 436 | 90 CONTINUE
|
---|
| 437 | *
|
---|
| 438 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 439 | *
|
---|
| 440 | * Form C * H or C * H**T where C = ( C1 C2 )
|
---|
| 441 | *
|
---|
| 442 | LASTV = MAX( K, ILASLR( N, K, V, LDV ) )
|
---|
| 443 | LASTC = ILASLR( M, LASTV, C, LDC )
|
---|
| 444 | *
|
---|
| 445 | * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
|
---|
| 446 | *
|
---|
| 447 | * W := C2
|
---|
| 448 | *
|
---|
| 449 | DO 100 J = 1, K
|
---|
| 450 | CALL SCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
|
---|
| 451 | 100 CONTINUE
|
---|
| 452 | *
|
---|
| 453 | * W := W * V2
|
---|
| 454 | *
|
---|
| 455 | CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 456 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 457 | $ WORK, LDWORK )
|
---|
| 458 | IF( LASTV.GT.K ) THEN
|
---|
| 459 | *
|
---|
| 460 | * W := W + C1 * V1
|
---|
| 461 | *
|
---|
| 462 | CALL SGEMM( 'No transpose', 'No transpose',
|
---|
| 463 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
|
---|
| 464 | $ ONE, WORK, LDWORK )
|
---|
| 465 | END IF
|
---|
| 466 | *
|
---|
| 467 | * W := W * T or W * T**T
|
---|
| 468 | *
|
---|
| 469 | CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit',
|
---|
| 470 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 471 | *
|
---|
| 472 | * C := C - W * V**T
|
---|
| 473 | *
|
---|
| 474 | IF( LASTV.GT.K ) THEN
|
---|
| 475 | *
|
---|
| 476 | * C1 := C1 - W * V1**T
|
---|
| 477 | *
|
---|
| 478 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 479 | $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
|
---|
| 480 | $ ONE, C, LDC )
|
---|
| 481 | END IF
|
---|
| 482 | *
|
---|
| 483 | * W := W * V2**T
|
---|
| 484 | *
|
---|
| 485 | CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
|
---|
| 486 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 487 | $ WORK, LDWORK )
|
---|
| 488 | *
|
---|
| 489 | * C2 := C2 - W
|
---|
| 490 | *
|
---|
| 491 | DO 120 J = 1, K
|
---|
| 492 | DO 110 I = 1, LASTC
|
---|
| 493 | C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
|
---|
| 494 | 110 CONTINUE
|
---|
| 495 | 120 CONTINUE
|
---|
| 496 | END IF
|
---|
| 497 | END IF
|
---|
| 498 | *
|
---|
| 499 | ELSE IF( LSAME( STOREV, 'R' ) ) THEN
|
---|
| 500 | *
|
---|
| 501 | IF( LSAME( DIRECT, 'F' ) ) THEN
|
---|
| 502 | *
|
---|
| 503 | * Let V = ( V1 V2 ) (V1: first K columns)
|
---|
| 504 | * where V1 is unit upper triangular.
|
---|
| 505 | *
|
---|
| 506 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 507 | *
|
---|
| 508 | * Form H * C or H**T * C where C = ( C1 )
|
---|
| 509 | * ( C2 )
|
---|
| 510 | *
|
---|
| 511 | LASTV = MAX( K, ILASLC( K, M, V, LDV ) )
|
---|
| 512 | LASTC = ILASLC( LASTV, N, C, LDC )
|
---|
| 513 | *
|
---|
| 514 | * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
|
---|
| 515 | *
|
---|
| 516 | * W := C1**T
|
---|
| 517 | *
|
---|
| 518 | DO 130 J = 1, K
|
---|
| 519 | CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
|
---|
| 520 | 130 CONTINUE
|
---|
| 521 | *
|
---|
| 522 | * W := W * V1**T
|
---|
| 523 | *
|
---|
| 524 | CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
|
---|
| 525 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 526 | IF( LASTV.GT.K ) THEN
|
---|
| 527 | *
|
---|
| 528 | * W := W + C2**T*V2**T
|
---|
| 529 | *
|
---|
| 530 | CALL SGEMM( 'Transpose', 'Transpose',
|
---|
| 531 | $ LASTC, K, LASTV-K,
|
---|
| 532 | $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
|
---|
| 533 | $ ONE, WORK, LDWORK )
|
---|
| 534 | END IF
|
---|
| 535 | *
|
---|
| 536 | * W := W * T**T or W * T
|
---|
| 537 | *
|
---|
| 538 | CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit',
|
---|
| 539 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 540 | *
|
---|
| 541 | * C := C - V**T * W**T
|
---|
| 542 | *
|
---|
| 543 | IF( LASTV.GT.K ) THEN
|
---|
| 544 | *
|
---|
| 545 | * C2 := C2 - V2**T * W**T
|
---|
| 546 | *
|
---|
| 547 | CALL SGEMM( 'Transpose', 'Transpose',
|
---|
| 548 | $ LASTV-K, LASTC, K,
|
---|
| 549 | $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
|
---|
| 550 | $ ONE, C( K+1, 1 ), LDC )
|
---|
| 551 | END IF
|
---|
| 552 | *
|
---|
| 553 | * W := W * V1
|
---|
| 554 | *
|
---|
| 555 | CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 556 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 557 | *
|
---|
| 558 | * C1 := C1 - W**T
|
---|
| 559 | *
|
---|
| 560 | DO 150 J = 1, K
|
---|
| 561 | DO 140 I = 1, LASTC
|
---|
| 562 | C( J, I ) = C( J, I ) - WORK( I, J )
|
---|
| 563 | 140 CONTINUE
|
---|
| 564 | 150 CONTINUE
|
---|
| 565 | *
|
---|
| 566 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 567 | *
|
---|
| 568 | * Form C * H or C * H**T where C = ( C1 C2 )
|
---|
| 569 | *
|
---|
| 570 | LASTV = MAX( K, ILASLC( K, N, V, LDV ) )
|
---|
| 571 | LASTC = ILASLR( M, LASTV, C, LDC )
|
---|
| 572 | *
|
---|
| 573 | * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
|
---|
| 574 | *
|
---|
| 575 | * W := C1
|
---|
| 576 | *
|
---|
| 577 | DO 160 J = 1, K
|
---|
| 578 | CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
|
---|
| 579 | 160 CONTINUE
|
---|
| 580 | *
|
---|
| 581 | * W := W * V1**T
|
---|
| 582 | *
|
---|
| 583 | CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
|
---|
| 584 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 585 | IF( LASTV.GT.K ) THEN
|
---|
| 586 | *
|
---|
| 587 | * W := W + C2 * V2**T
|
---|
| 588 | *
|
---|
| 589 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 590 | $ LASTC, K, LASTV-K,
|
---|
| 591 | $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
|
---|
| 592 | $ ONE, WORK, LDWORK )
|
---|
| 593 | END IF
|
---|
| 594 | *
|
---|
| 595 | * W := W * T or W * T**T
|
---|
| 596 | *
|
---|
| 597 | CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit',
|
---|
| 598 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 599 | *
|
---|
| 600 | * C := C - W * V
|
---|
| 601 | *
|
---|
| 602 | IF( LASTV.GT.K ) THEN
|
---|
| 603 | *
|
---|
| 604 | * C2 := C2 - W * V2
|
---|
| 605 | *
|
---|
| 606 | CALL SGEMM( 'No transpose', 'No transpose',
|
---|
| 607 | $ LASTC, LASTV-K, K,
|
---|
| 608 | $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
|
---|
| 609 | $ ONE, C( 1, K+1 ), LDC )
|
---|
| 610 | END IF
|
---|
| 611 | *
|
---|
| 612 | * W := W * V1
|
---|
| 613 | *
|
---|
| 614 | CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 615 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 616 | *
|
---|
| 617 | * C1 := C1 - W
|
---|
| 618 | *
|
---|
| 619 | DO 180 J = 1, K
|
---|
| 620 | DO 170 I = 1, LASTC
|
---|
| 621 | C( I, J ) = C( I, J ) - WORK( I, J )
|
---|
| 622 | 170 CONTINUE
|
---|
| 623 | 180 CONTINUE
|
---|
| 624 | *
|
---|
| 625 | END IF
|
---|
| 626 | *
|
---|
| 627 | ELSE
|
---|
| 628 | *
|
---|
| 629 | * Let V = ( V1 V2 ) (V2: last K columns)
|
---|
| 630 | * where V2 is unit lower triangular.
|
---|
| 631 | *
|
---|
| 632 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 633 | *
|
---|
| 634 | * Form H * C or H**T * C where C = ( C1 )
|
---|
| 635 | * ( C2 )
|
---|
| 636 | *
|
---|
| 637 | LASTV = MAX( K, ILASLC( K, M, V, LDV ) )
|
---|
| 638 | LASTC = ILASLC( LASTV, N, C, LDC )
|
---|
| 639 | *
|
---|
| 640 | * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
|
---|
| 641 | *
|
---|
| 642 | * W := C2**T
|
---|
| 643 | *
|
---|
| 644 | DO 190 J = 1, K
|
---|
| 645 | CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
|
---|
| 646 | $ WORK( 1, J ), 1 )
|
---|
| 647 | 190 CONTINUE
|
---|
| 648 | *
|
---|
| 649 | * W := W * V2**T
|
---|
| 650 | *
|
---|
| 651 | CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
|
---|
| 652 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 653 | $ WORK, LDWORK )
|
---|
| 654 | IF( LASTV.GT.K ) THEN
|
---|
| 655 | *
|
---|
| 656 | * W := W + C1**T * V1**T
|
---|
| 657 | *
|
---|
| 658 | CALL SGEMM( 'Transpose', 'Transpose',
|
---|
| 659 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
|
---|
| 660 | $ ONE, WORK, LDWORK )
|
---|
| 661 | END IF
|
---|
| 662 | *
|
---|
| 663 | * W := W * T**T or W * T
|
---|
| 664 | *
|
---|
| 665 | CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit',
|
---|
| 666 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 667 | *
|
---|
| 668 | * C := C - V**T * W**T
|
---|
| 669 | *
|
---|
| 670 | IF( LASTV.GT.K ) THEN
|
---|
| 671 | *
|
---|
| 672 | * C1 := C1 - V1**T * W**T
|
---|
| 673 | *
|
---|
| 674 | CALL SGEMM( 'Transpose', 'Transpose',
|
---|
| 675 | $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
|
---|
| 676 | $ ONE, C, LDC )
|
---|
| 677 | END IF
|
---|
| 678 | *
|
---|
| 679 | * W := W * V2
|
---|
| 680 | *
|
---|
| 681 | CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 682 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 683 | $ WORK, LDWORK )
|
---|
| 684 | *
|
---|
| 685 | * C2 := C2 - W**T
|
---|
| 686 | *
|
---|
| 687 | DO 210 J = 1, K
|
---|
| 688 | DO 200 I = 1, LASTC
|
---|
| 689 | C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
|
---|
| 690 | 200 CONTINUE
|
---|
| 691 | 210 CONTINUE
|
---|
| 692 | *
|
---|
| 693 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 694 | *
|
---|
| 695 | * Form C * H or C * H**T where C = ( C1 C2 )
|
---|
| 696 | *
|
---|
| 697 | LASTV = MAX( K, ILASLC( K, N, V, LDV ) )
|
---|
| 698 | LASTC = ILASLR( M, LASTV, C, LDC )
|
---|
| 699 | *
|
---|
| 700 | * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
|
---|
| 701 | *
|
---|
| 702 | * W := C2
|
---|
| 703 | *
|
---|
| 704 | DO 220 J = 1, K
|
---|
| 705 | CALL SCOPY( LASTC, C( 1, LASTV-K+J ), 1,
|
---|
| 706 | $ WORK( 1, J ), 1 )
|
---|
| 707 | 220 CONTINUE
|
---|
| 708 | *
|
---|
| 709 | * W := W * V2**T
|
---|
| 710 | *
|
---|
| 711 | CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
|
---|
| 712 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 713 | $ WORK, LDWORK )
|
---|
| 714 | IF( LASTV.GT.K ) THEN
|
---|
| 715 | *
|
---|
| 716 | * W := W + C1 * V1**T
|
---|
| 717 | *
|
---|
| 718 | CALL SGEMM( 'No transpose', 'Transpose',
|
---|
| 719 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
|
---|
| 720 | $ ONE, WORK, LDWORK )
|
---|
| 721 | END IF
|
---|
| 722 | *
|
---|
| 723 | * W := W * T or W * T**T
|
---|
| 724 | *
|
---|
| 725 | CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit',
|
---|
| 726 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 727 | *
|
---|
| 728 | * C := C - W * V
|
---|
| 729 | *
|
---|
| 730 | IF( LASTV.GT.K ) THEN
|
---|
| 731 | *
|
---|
| 732 | * C1 := C1 - W * V1
|
---|
| 733 | *
|
---|
| 734 | CALL SGEMM( 'No transpose', 'No transpose',
|
---|
| 735 | $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
|
---|
| 736 | $ ONE, C, LDC )
|
---|
| 737 | END IF
|
---|
| 738 | *
|
---|
| 739 | * W := W * V2
|
---|
| 740 | *
|
---|
| 741 | CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 742 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 743 | $ WORK, LDWORK )
|
---|
| 744 | *
|
---|
| 745 | * C1 := C1 - W
|
---|
| 746 | *
|
---|
| 747 | DO 240 J = 1, K
|
---|
| 748 | DO 230 I = 1, LASTC
|
---|
| 749 | C( I, LASTV-K+J ) = C( I, LASTV-K+J )
|
---|
| 750 | $ - WORK( I, J )
|
---|
| 751 | 230 CONTINUE
|
---|
| 752 | 240 CONTINUE
|
---|
| 753 | *
|
---|
| 754 | END IF
|
---|
| 755 | *
|
---|
| 756 | END IF
|
---|
| 757 | END IF
|
---|
| 758 | *
|
---|
| 759 | RETURN
|
---|
| 760 | *
|
---|
| 761 | * End of SLARFB
|
---|
| 762 | *
|
---|
| 763 | END
|
---|