[136] | 1 | *> \brief \b CLARFB
|
---|
| 2 | *
|
---|
| 3 | * =========== DOCUMENTATION ===========
|
---|
| 4 | *
|
---|
| 5 | * Online html documentation available at
|
---|
| 6 | * http://www.netlib.org/lapack/explore-html/
|
---|
| 7 | *
|
---|
| 8 | *> \htmlonly
|
---|
| 9 | *> Download CLARFB + dependencies
|
---|
| 10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb.f">
|
---|
| 11 | *> [TGZ]</a>
|
---|
| 12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb.f">
|
---|
| 13 | *> [ZIP]</a>
|
---|
| 14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb.f">
|
---|
| 15 | *> [TXT]</a>
|
---|
| 16 | *> \endhtmlonly
|
---|
| 17 | *
|
---|
| 18 | * Definition:
|
---|
| 19 | * ===========
|
---|
| 20 | *
|
---|
| 21 | * SUBROUTINE CLARFB( 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 | * COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
|
---|
| 30 | * $ WORK( LDWORK, * )
|
---|
| 31 | * ..
|
---|
| 32 | *
|
---|
| 33 | *
|
---|
| 34 | *> \par Purpose:
|
---|
| 35 | * =============
|
---|
| 36 | *>
|
---|
| 37 | *> \verbatim
|
---|
| 38 | *>
|
---|
| 39 | *> CLARFB applies a complex block reflector H or its transpose H**H to a
|
---|
| 40 | *> complex 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**H from the Left
|
---|
| 50 | *> = 'R': apply H or H**H from the Right
|
---|
| 51 | *> \endverbatim
|
---|
| 52 | *>
|
---|
| 53 | *> \param[in] TRANS
|
---|
| 54 | *> \verbatim
|
---|
| 55 | *> TRANS is CHARACTER*1
|
---|
| 56 | *> = 'N': apply H (No transpose)
|
---|
| 57 | *> = 'C': apply H**H (Conjugate 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 COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDC,N)
|
---|
| 131 | *> On entry, the M-by-N matrix C.
|
---|
| 132 | *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
|
---|
| 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 COMPLEX 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 complexOTHERauxiliary
|
---|
| 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 CLARFB( 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 | COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
|
---|
| 209 | $ WORK( LDWORK, * )
|
---|
| 210 | * ..
|
---|
| 211 | *
|
---|
| 212 | * =====================================================================
|
---|
| 213 | *
|
---|
| 214 | * .. Parameters ..
|
---|
| 215 | COMPLEX ONE
|
---|
| 216 | PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
|
---|
| 217 | * ..
|
---|
| 218 | * .. Local Scalars ..
|
---|
| 219 | CHARACTER TRANST
|
---|
| 220 | INTEGER I, J, LASTV, LASTC
|
---|
| 221 | * ..
|
---|
| 222 | * .. External Functions ..
|
---|
| 223 | LOGICAL LSAME
|
---|
| 224 | INTEGER ILACLR, ILACLC
|
---|
| 225 | EXTERNAL LSAME, ILACLR, ILACLC
|
---|
| 226 | * ..
|
---|
| 227 | * .. External Subroutines ..
|
---|
| 228 | EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM
|
---|
| 229 | * ..
|
---|
| 230 | * .. Intrinsic Functions ..
|
---|
| 231 | INTRINSIC CONJG
|
---|
| 232 | * ..
|
---|
| 233 | * .. Executable Statements ..
|
---|
| 234 | *
|
---|
| 235 | * Quick return if possible
|
---|
| 236 | *
|
---|
| 237 | IF( M.LE.0 .OR. N.LE.0 )
|
---|
| 238 | $ RETURN
|
---|
| 239 | *
|
---|
| 240 | IF( LSAME( TRANS, 'N' ) ) THEN
|
---|
| 241 | TRANST = 'C'
|
---|
| 242 | ELSE
|
---|
| 243 | TRANST = 'N'
|
---|
| 244 | END IF
|
---|
| 245 | *
|
---|
| 246 | IF( LSAME( STOREV, 'C' ) ) THEN
|
---|
| 247 | *
|
---|
| 248 | IF( LSAME( DIRECT, 'F' ) ) THEN
|
---|
| 249 | *
|
---|
| 250 | * Let V = ( V1 ) (first K rows)
|
---|
| 251 | * ( V2 )
|
---|
| 252 | * where V1 is unit lower triangular.
|
---|
| 253 | *
|
---|
| 254 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 255 | *
|
---|
| 256 | * Form H * C or H**H * C where C = ( C1 )
|
---|
| 257 | * ( C2 )
|
---|
| 258 | *
|
---|
| 259 | LASTV = MAX( K, ILACLR( M, K, V, LDV ) )
|
---|
| 260 | LASTC = ILACLC( LASTV, N, C, LDC )
|
---|
| 261 | *
|
---|
| 262 | * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
|
---|
| 263 | *
|
---|
| 264 | * W := C1**H
|
---|
| 265 | *
|
---|
| 266 | DO 10 J = 1, K
|
---|
| 267 | CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
|
---|
| 268 | CALL CLACGV( LASTC, WORK( 1, J ), 1 )
|
---|
| 269 | 10 CONTINUE
|
---|
| 270 | *
|
---|
| 271 | * W := W * V1
|
---|
| 272 | *
|
---|
| 273 | CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 274 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 275 | IF( LASTV.GT.K ) THEN
|
---|
| 276 | *
|
---|
| 277 | * W := W + C2**H *V2
|
---|
| 278 | *
|
---|
| 279 | CALL CGEMM( 'Conjugate transpose', 'No transpose',
|
---|
| 280 | $ LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC,
|
---|
| 281 | $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
|
---|
| 282 | END IF
|
---|
| 283 | *
|
---|
| 284 | * W := W * T**H or W * T
|
---|
| 285 | *
|
---|
| 286 | CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
|
---|
| 287 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 288 | *
|
---|
| 289 | * C := C - V * W**H
|
---|
| 290 | *
|
---|
| 291 | IF( M.GT.K ) THEN
|
---|
| 292 | *
|
---|
| 293 | * C2 := C2 - V2 * W**H
|
---|
| 294 | *
|
---|
| 295 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 296 | $ LASTV-K, LASTC, K, -ONE, V( K+1, 1 ), LDV,
|
---|
| 297 | $ WORK, LDWORK, ONE, C( K+1, 1 ), LDC )
|
---|
| 298 | END IF
|
---|
| 299 | *
|
---|
| 300 | * W := W * V1**H
|
---|
| 301 | *
|
---|
| 302 | CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
|
---|
| 303 | $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 304 | *
|
---|
| 305 | * C1 := C1 - W**H
|
---|
| 306 | *
|
---|
| 307 | DO 30 J = 1, K
|
---|
| 308 | DO 20 I = 1, LASTC
|
---|
| 309 | C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
|
---|
| 310 | 20 CONTINUE
|
---|
| 311 | 30 CONTINUE
|
---|
| 312 | *
|
---|
| 313 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 314 | *
|
---|
| 315 | * Form C * H or C * H**H where C = ( C1 C2 )
|
---|
| 316 | *
|
---|
| 317 | LASTV = MAX( K, ILACLR( N, K, V, LDV ) )
|
---|
| 318 | LASTC = ILACLR( M, LASTV, C, LDC )
|
---|
| 319 | *
|
---|
| 320 | * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
|
---|
| 321 | *
|
---|
| 322 | * W := C1
|
---|
| 323 | *
|
---|
| 324 | DO 40 J = 1, K
|
---|
| 325 | CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
|
---|
| 326 | 40 CONTINUE
|
---|
| 327 | *
|
---|
| 328 | * W := W * V1
|
---|
| 329 | *
|
---|
| 330 | CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 331 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 332 | IF( LASTV.GT.K ) THEN
|
---|
| 333 | *
|
---|
| 334 | * W := W + C2 * V2
|
---|
| 335 | *
|
---|
| 336 | CALL CGEMM( 'No transpose', 'No transpose',
|
---|
| 337 | $ LASTC, K, LASTV-K,
|
---|
| 338 | $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
|
---|
| 339 | $ ONE, WORK, LDWORK )
|
---|
| 340 | END IF
|
---|
| 341 | *
|
---|
| 342 | * W := W * T or W * T**H
|
---|
| 343 | *
|
---|
| 344 | CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
|
---|
| 345 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 346 | *
|
---|
| 347 | * C := C - W * V**H
|
---|
| 348 | *
|
---|
| 349 | IF( LASTV.GT.K ) THEN
|
---|
| 350 | *
|
---|
| 351 | * C2 := C2 - W * V2**H
|
---|
| 352 | *
|
---|
| 353 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 354 | $ LASTC, LASTV-K, K,
|
---|
| 355 | $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV,
|
---|
| 356 | $ ONE, C( 1, K+1 ), LDC )
|
---|
| 357 | END IF
|
---|
| 358 | *
|
---|
| 359 | * W := W * V1**H
|
---|
| 360 | *
|
---|
| 361 | CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
|
---|
| 362 | $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 363 | *
|
---|
| 364 | * C1 := C1 - W
|
---|
| 365 | *
|
---|
| 366 | DO 60 J = 1, K
|
---|
| 367 | DO 50 I = 1, LASTC
|
---|
| 368 | C( I, J ) = C( I, J ) - WORK( I, J )
|
---|
| 369 | 50 CONTINUE
|
---|
| 370 | 60 CONTINUE
|
---|
| 371 | END IF
|
---|
| 372 | *
|
---|
| 373 | ELSE
|
---|
| 374 | *
|
---|
| 375 | * Let V = ( V1 )
|
---|
| 376 | * ( V2 ) (last K rows)
|
---|
| 377 | * where V2 is unit upper triangular.
|
---|
| 378 | *
|
---|
| 379 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 380 | *
|
---|
| 381 | * Form H * C or H**H * C where C = ( C1 )
|
---|
| 382 | * ( C2 )
|
---|
| 383 | *
|
---|
| 384 | LASTV = MAX( K, ILACLR( M, K, V, LDV ) )
|
---|
| 385 | LASTC = ILACLC( LASTV, N, C, LDC )
|
---|
| 386 | *
|
---|
| 387 | * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
|
---|
| 388 | *
|
---|
| 389 | * W := C2**H
|
---|
| 390 | *
|
---|
| 391 | DO 70 J = 1, K
|
---|
| 392 | CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
|
---|
| 393 | $ WORK( 1, J ), 1 )
|
---|
| 394 | CALL CLACGV( LASTC, WORK( 1, J ), 1 )
|
---|
| 395 | 70 CONTINUE
|
---|
| 396 | *
|
---|
| 397 | * W := W * V2
|
---|
| 398 | *
|
---|
| 399 | CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 400 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 401 | $ WORK, LDWORK )
|
---|
| 402 | IF( LASTV.GT.K ) THEN
|
---|
| 403 | *
|
---|
| 404 | * W := W + C1**H*V1
|
---|
| 405 | *
|
---|
| 406 | CALL CGEMM( 'Conjugate transpose', 'No transpose',
|
---|
| 407 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
|
---|
| 408 | $ ONE, WORK, LDWORK )
|
---|
| 409 | END IF
|
---|
| 410 | *
|
---|
| 411 | * W := W * T**H or W * T
|
---|
| 412 | *
|
---|
| 413 | CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
|
---|
| 414 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 415 | *
|
---|
| 416 | * C := C - V * W**H
|
---|
| 417 | *
|
---|
| 418 | IF( LASTV.GT.K ) THEN
|
---|
| 419 | *
|
---|
| 420 | * C1 := C1 - V1 * W**H
|
---|
| 421 | *
|
---|
| 422 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 423 | $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
|
---|
| 424 | $ ONE, C, LDC )
|
---|
| 425 | END IF
|
---|
| 426 | *
|
---|
| 427 | * W := W * V2**H
|
---|
| 428 | *
|
---|
| 429 | CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
|
---|
| 430 | $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 431 | $ WORK, LDWORK )
|
---|
| 432 | *
|
---|
| 433 | * C2 := C2 - W**H
|
---|
| 434 | *
|
---|
| 435 | DO 90 J = 1, K
|
---|
| 436 | DO 80 I = 1, LASTC
|
---|
| 437 | C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
|
---|
| 438 | $ CONJG( WORK( I, J ) )
|
---|
| 439 | 80 CONTINUE
|
---|
| 440 | 90 CONTINUE
|
---|
| 441 | *
|
---|
| 442 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 443 | *
|
---|
| 444 | * Form C * H or C * H**H where C = ( C1 C2 )
|
---|
| 445 | *
|
---|
| 446 | LASTV = MAX( K, ILACLR( N, K, V, LDV ) )
|
---|
| 447 | LASTC = ILACLR( M, LASTV, C, LDC )
|
---|
| 448 | *
|
---|
| 449 | * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
|
---|
| 450 | *
|
---|
| 451 | * W := C2
|
---|
| 452 | *
|
---|
| 453 | DO 100 J = 1, K
|
---|
| 454 | CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1,
|
---|
| 455 | $ WORK( 1, J ), 1 )
|
---|
| 456 | 100 CONTINUE
|
---|
| 457 | *
|
---|
| 458 | * W := W * V2
|
---|
| 459 | *
|
---|
| 460 | CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 461 | $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 462 | $ WORK, LDWORK )
|
---|
| 463 | IF( LASTV.GT.K ) THEN
|
---|
| 464 | *
|
---|
| 465 | * W := W + C1 * V1
|
---|
| 466 | *
|
---|
| 467 | CALL CGEMM( 'No transpose', 'No transpose',
|
---|
| 468 | $ LASTC, K, LASTV-K,
|
---|
| 469 | $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
|
---|
| 470 | END IF
|
---|
| 471 | *
|
---|
| 472 | * W := W * T or W * T**H
|
---|
| 473 | *
|
---|
| 474 | CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
|
---|
| 475 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 476 | *
|
---|
| 477 | * C := C - W * V**H
|
---|
| 478 | *
|
---|
| 479 | IF( LASTV.GT.K ) THEN
|
---|
| 480 | *
|
---|
| 481 | * C1 := C1 - W * V1**H
|
---|
| 482 | *
|
---|
| 483 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 484 | $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
|
---|
| 485 | $ ONE, C, LDC )
|
---|
| 486 | END IF
|
---|
| 487 | *
|
---|
| 488 | * W := W * V2**H
|
---|
| 489 | *
|
---|
| 490 | CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
|
---|
| 491 | $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
|
---|
| 492 | $ WORK, LDWORK )
|
---|
| 493 | *
|
---|
| 494 | * C2 := C2 - W
|
---|
| 495 | *
|
---|
| 496 | DO 120 J = 1, K
|
---|
| 497 | DO 110 I = 1, LASTC
|
---|
| 498 | C( I, LASTV-K+J ) = C( I, LASTV-K+J )
|
---|
| 499 | $ - WORK( I, J )
|
---|
| 500 | 110 CONTINUE
|
---|
| 501 | 120 CONTINUE
|
---|
| 502 | END IF
|
---|
| 503 | END IF
|
---|
| 504 | *
|
---|
| 505 | ELSE IF( LSAME( STOREV, 'R' ) ) THEN
|
---|
| 506 | *
|
---|
| 507 | IF( LSAME( DIRECT, 'F' ) ) THEN
|
---|
| 508 | *
|
---|
| 509 | * Let V = ( V1 V2 ) (V1: first K columns)
|
---|
| 510 | * where V1 is unit upper triangular.
|
---|
| 511 | *
|
---|
| 512 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 513 | *
|
---|
| 514 | * Form H * C or H**H * C where C = ( C1 )
|
---|
| 515 | * ( C2 )
|
---|
| 516 | *
|
---|
| 517 | LASTV = MAX( K, ILACLC( K, M, V, LDV ) )
|
---|
| 518 | LASTC = ILACLC( LASTV, N, C, LDC )
|
---|
| 519 | *
|
---|
| 520 | * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
|
---|
| 521 | *
|
---|
| 522 | * W := C1**H
|
---|
| 523 | *
|
---|
| 524 | DO 130 J = 1, K
|
---|
| 525 | CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
|
---|
| 526 | CALL CLACGV( LASTC, WORK( 1, J ), 1 )
|
---|
| 527 | 130 CONTINUE
|
---|
| 528 | *
|
---|
| 529 | * W := W * V1**H
|
---|
| 530 | *
|
---|
| 531 | CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
|
---|
| 532 | $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 533 | IF( LASTV.GT.K ) THEN
|
---|
| 534 | *
|
---|
| 535 | * W := W + C2**H*V2**H
|
---|
| 536 | *
|
---|
| 537 | CALL CGEMM( 'Conjugate transpose',
|
---|
| 538 | $ 'Conjugate transpose', LASTC, K, LASTV-K,
|
---|
| 539 | $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
|
---|
| 540 | $ ONE, WORK, LDWORK )
|
---|
| 541 | END IF
|
---|
| 542 | *
|
---|
| 543 | * W := W * T**H or W * T
|
---|
| 544 | *
|
---|
| 545 | CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
|
---|
| 546 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 547 | *
|
---|
| 548 | * C := C - V**H * W**H
|
---|
| 549 | *
|
---|
| 550 | IF( LASTV.GT.K ) THEN
|
---|
| 551 | *
|
---|
| 552 | * C2 := C2 - V2**H * W**H
|
---|
| 553 | *
|
---|
| 554 | CALL CGEMM( 'Conjugate transpose',
|
---|
| 555 | $ 'Conjugate transpose', LASTV-K, LASTC, K,
|
---|
| 556 | $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
|
---|
| 557 | $ ONE, C( K+1, 1 ), LDC )
|
---|
| 558 | END IF
|
---|
| 559 | *
|
---|
| 560 | * W := W * V1
|
---|
| 561 | *
|
---|
| 562 | CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 563 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 564 | *
|
---|
| 565 | * C1 := C1 - W**H
|
---|
| 566 | *
|
---|
| 567 | DO 150 J = 1, K
|
---|
| 568 | DO 140 I = 1, LASTC
|
---|
| 569 | C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
|
---|
| 570 | 140 CONTINUE
|
---|
| 571 | 150 CONTINUE
|
---|
| 572 | *
|
---|
| 573 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 574 | *
|
---|
| 575 | * Form C * H or C * H**H where C = ( C1 C2 )
|
---|
| 576 | *
|
---|
| 577 | LASTV = MAX( K, ILACLC( K, N, V, LDV ) )
|
---|
| 578 | LASTC = ILACLR( M, LASTV, C, LDC )
|
---|
| 579 | *
|
---|
| 580 | * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
|
---|
| 581 | *
|
---|
| 582 | * W := C1
|
---|
| 583 | *
|
---|
| 584 | DO 160 J = 1, K
|
---|
| 585 | CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
|
---|
| 586 | 160 CONTINUE
|
---|
| 587 | *
|
---|
| 588 | * W := W * V1**H
|
---|
| 589 | *
|
---|
| 590 | CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
|
---|
| 591 | $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 592 | IF( LASTV.GT.K ) THEN
|
---|
| 593 | *
|
---|
| 594 | * W := W + C2 * V2**H
|
---|
| 595 | *
|
---|
| 596 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 597 | $ LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC,
|
---|
| 598 | $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
|
---|
| 599 | END IF
|
---|
| 600 | *
|
---|
| 601 | * W := W * T or W * T**H
|
---|
| 602 | *
|
---|
| 603 | CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
|
---|
| 604 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 605 | *
|
---|
| 606 | * C := C - W * V
|
---|
| 607 | *
|
---|
| 608 | IF( LASTV.GT.K ) THEN
|
---|
| 609 | *
|
---|
| 610 | * C2 := C2 - W * V2
|
---|
| 611 | *
|
---|
| 612 | CALL CGEMM( 'No transpose', 'No transpose',
|
---|
| 613 | $ LASTC, LASTV-K, K,
|
---|
| 614 | $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
|
---|
| 615 | $ ONE, C( 1, K+1 ), LDC )
|
---|
| 616 | END IF
|
---|
| 617 | *
|
---|
| 618 | * W := W * V1
|
---|
| 619 | *
|
---|
| 620 | CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
|
---|
| 621 | $ LASTC, K, ONE, V, LDV, WORK, LDWORK )
|
---|
| 622 | *
|
---|
| 623 | * C1 := C1 - W
|
---|
| 624 | *
|
---|
| 625 | DO 180 J = 1, K
|
---|
| 626 | DO 170 I = 1, LASTC
|
---|
| 627 | C( I, J ) = C( I, J ) - WORK( I, J )
|
---|
| 628 | 170 CONTINUE
|
---|
| 629 | 180 CONTINUE
|
---|
| 630 | *
|
---|
| 631 | END IF
|
---|
| 632 | *
|
---|
| 633 | ELSE
|
---|
| 634 | *
|
---|
| 635 | * Let V = ( V1 V2 ) (V2: last K columns)
|
---|
| 636 | * where V2 is unit lower triangular.
|
---|
| 637 | *
|
---|
| 638 | IF( LSAME( SIDE, 'L' ) ) THEN
|
---|
| 639 | *
|
---|
| 640 | * Form H * C or H**H * C where C = ( C1 )
|
---|
| 641 | * ( C2 )
|
---|
| 642 | *
|
---|
| 643 | LASTV = MAX( K, ILACLC( K, M, V, LDV ) )
|
---|
| 644 | LASTC = ILACLC( LASTV, N, C, LDC )
|
---|
| 645 | *
|
---|
| 646 | * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
|
---|
| 647 | *
|
---|
| 648 | * W := C2**H
|
---|
| 649 | *
|
---|
| 650 | DO 190 J = 1, K
|
---|
| 651 | CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
|
---|
| 652 | $ WORK( 1, J ), 1 )
|
---|
| 653 | CALL CLACGV( LASTC, WORK( 1, J ), 1 )
|
---|
| 654 | 190 CONTINUE
|
---|
| 655 | *
|
---|
| 656 | * W := W * V2**H
|
---|
| 657 | *
|
---|
| 658 | CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
|
---|
| 659 | $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 660 | $ WORK, LDWORK )
|
---|
| 661 | IF( LASTV.GT.K ) THEN
|
---|
| 662 | *
|
---|
| 663 | * W := W + C1**H * V1**H
|
---|
| 664 | *
|
---|
| 665 | CALL CGEMM( 'Conjugate transpose',
|
---|
| 666 | $ 'Conjugate transpose', LASTC, K, LASTV-K,
|
---|
| 667 | $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
|
---|
| 668 | END IF
|
---|
| 669 | *
|
---|
| 670 | * W := W * T**H or W * T
|
---|
| 671 | *
|
---|
| 672 | CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
|
---|
| 673 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 674 | *
|
---|
| 675 | * C := C - V**H * W**H
|
---|
| 676 | *
|
---|
| 677 | IF( LASTV.GT.K ) THEN
|
---|
| 678 | *
|
---|
| 679 | * C1 := C1 - V1**H * W**H
|
---|
| 680 | *
|
---|
| 681 | CALL CGEMM( 'Conjugate transpose',
|
---|
| 682 | $ 'Conjugate transpose', LASTV-K, LASTC, K,
|
---|
| 683 | $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
|
---|
| 684 | END IF
|
---|
| 685 | *
|
---|
| 686 | * W := W * V2
|
---|
| 687 | *
|
---|
| 688 | CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 689 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 690 | $ WORK, LDWORK )
|
---|
| 691 | *
|
---|
| 692 | * C2 := C2 - W**H
|
---|
| 693 | *
|
---|
| 694 | DO 210 J = 1, K
|
---|
| 695 | DO 200 I = 1, LASTC
|
---|
| 696 | C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
|
---|
| 697 | $ CONJG( WORK( I, J ) )
|
---|
| 698 | 200 CONTINUE
|
---|
| 699 | 210 CONTINUE
|
---|
| 700 | *
|
---|
| 701 | ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
---|
| 702 | *
|
---|
| 703 | * Form C * H or C * H**H where C = ( C1 C2 )
|
---|
| 704 | *
|
---|
| 705 | LASTV = MAX( K, ILACLC( K, N, V, LDV ) )
|
---|
| 706 | LASTC = ILACLR( M, LASTV, C, LDC )
|
---|
| 707 | *
|
---|
| 708 | * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
|
---|
| 709 | *
|
---|
| 710 | * W := C2
|
---|
| 711 | *
|
---|
| 712 | DO 220 J = 1, K
|
---|
| 713 | CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1,
|
---|
| 714 | $ WORK( 1, J ), 1 )
|
---|
| 715 | 220 CONTINUE
|
---|
| 716 | *
|
---|
| 717 | * W := W * V2**H
|
---|
| 718 | *
|
---|
| 719 | CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
|
---|
| 720 | $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 721 | $ WORK, LDWORK )
|
---|
| 722 | IF( LASTV.GT.K ) THEN
|
---|
| 723 | *
|
---|
| 724 | * W := W + C1 * V1**H
|
---|
| 725 | *
|
---|
| 726 | CALL CGEMM( 'No transpose', 'Conjugate transpose',
|
---|
| 727 | $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE,
|
---|
| 728 | $ WORK, LDWORK )
|
---|
| 729 | END IF
|
---|
| 730 | *
|
---|
| 731 | * W := W * T or W * T**H
|
---|
| 732 | *
|
---|
| 733 | CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
|
---|
| 734 | $ LASTC, K, ONE, T, LDT, WORK, LDWORK )
|
---|
| 735 | *
|
---|
| 736 | * C := C - W * V
|
---|
| 737 | *
|
---|
| 738 | IF( LASTV.GT.K ) THEN
|
---|
| 739 | *
|
---|
| 740 | * C1 := C1 - W * V1
|
---|
| 741 | *
|
---|
| 742 | CALL CGEMM( 'No transpose', 'No transpose',
|
---|
| 743 | $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
|
---|
| 744 | $ ONE, C, LDC )
|
---|
| 745 | END IF
|
---|
| 746 | *
|
---|
| 747 | * W := W * V2
|
---|
| 748 | *
|
---|
| 749 | CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
|
---|
| 750 | $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
|
---|
| 751 | $ WORK, LDWORK )
|
---|
| 752 | *
|
---|
| 753 | * C1 := C1 - W
|
---|
| 754 | *
|
---|
| 755 | DO 240 J = 1, K
|
---|
| 756 | DO 230 I = 1, LASTC
|
---|
| 757 | C( I, LASTV-K+J ) = C( I, LASTV-K+J )
|
---|
| 758 | $ - WORK( I, J )
|
---|
| 759 | 230 CONTINUE
|
---|
| 760 | 240 CONTINUE
|
---|
| 761 | *
|
---|
| 762 | END IF
|
---|
| 763 | *
|
---|
| 764 | END IF
|
---|
| 765 | END IF
|
---|
| 766 | *
|
---|
| 767 | RETURN
|
---|
| 768 | *
|
---|
| 769 | * End of CLARFB
|
---|
| 770 | *
|
---|
| 771 | END
|
---|