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