Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 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 |