LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
clarfb.f
Go to the documentation of this file.
1 *> \brief \b CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
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 *> If SIDE = 'L', M >= K >= 0;
96 *> if SIDE = 'R', N >= K >= 0.
97 *> \endverbatim
98 *>
99 *> \param[in] V
100 *> \verbatim
101 *> V is COMPLEX array, dimension
102 *> (LDV,K) if STOREV = 'C'
103 *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
104 *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
105 *> The matrix V. See Further Details.
106 *> \endverbatim
107 *>
108 *> \param[in] LDV
109 *> \verbatim
110 *> LDV is INTEGER
111 *> The leading dimension of the array V.
112 *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
113 *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
114 *> if STOREV = 'R', LDV >= K.
115 *> \endverbatim
116 *>
117 *> \param[in] T
118 *> \verbatim
119 *> T is COMPLEX array, dimension (LDT,K)
120 *> The triangular K-by-K matrix T in the representation of the
121 *> block reflector.
122 *> \endverbatim
123 *>
124 *> \param[in] LDT
125 *> \verbatim
126 *> LDT is INTEGER
127 *> The leading dimension of the array T. LDT >= K.
128 *> \endverbatim
129 *>
130 *> \param[in,out] C
131 *> \verbatim
132 *> C is COMPLEX array, dimension (LDC,N)
133 *> On entry, the M-by-N matrix C.
134 *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
135 *> \endverbatim
136 *>
137 *> \param[in] LDC
138 *> \verbatim
139 *> LDC is INTEGER
140 *> The leading dimension of the array C. LDC >= max(1,M).
141 *> \endverbatim
142 *>
143 *> \param[out] WORK
144 *> \verbatim
145 *> WORK is COMPLEX array, dimension (LDWORK,K)
146 *> \endverbatim
147 *>
148 *> \param[in] LDWORK
149 *> \verbatim
150 *> LDWORK is INTEGER
151 *> The leading dimension of the array WORK.
152 *> If SIDE = 'L', LDWORK >= max(1,N);
153 *> if SIDE = 'R', LDWORK >= max(1,M).
154 *> \endverbatim
155 *
156 * Authors:
157 * ========
158 *
159 *> \author Univ. of Tennessee
160 *> \author Univ. of California Berkeley
161 *> \author Univ. of Colorado Denver
162 *> \author NAG Ltd.
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 --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 *
202 * .. Scalar Arguments ..
203  CHARACTER DIRECT, SIDE, STOREV, TRANS
204  INTEGER K, LDC, LDT, LDV, LDWORK, M, N
205 * ..
206 * .. Array Arguments ..
207  COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
208  $ work( ldwork, * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  COMPLEX ONE
215  parameter( one = ( 1.0e+0, 0.0e+0 ) )
216 * ..
217 * .. Local Scalars ..
218  CHARACTER TRANST
219  INTEGER I, J
220 * ..
221 * .. External Functions ..
222  LOGICAL LSAME
223  EXTERNAL lsame
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL ccopy, cgemm, clacgv, ctrmm
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC conjg
230 * ..
231 * .. Executable Statements ..
232 *
233 * Quick return if possible
234 *
235  IF( m.LE.0 .OR. n.LE.0 )
236  $ RETURN
237 *
238  IF( lsame( trans, 'N' ) ) THEN
239  transt = 'C'
240  ELSE
241  transt = 'N'
242  END IF
243 *
244  IF( lsame( storev, 'C' ) ) THEN
245 *
246  IF( lsame( direct, 'F' ) ) THEN
247 *
248 * Let V = ( V1 ) (first K rows)
249 * ( V2 )
250 * where V1 is unit lower triangular.
251 *
252  IF( lsame( side, 'L' ) ) THEN
253 *
254 * Form H * C or H**H * C where C = ( C1 )
255 * ( C2 )
256 *
257 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
258 *
259 * W := C1**H
260 *
261  DO 10 j = 1, k
262  CALL ccopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
263  CALL clacgv( n, work( 1, j ), 1 )
264  10 CONTINUE
265 *
266 * W := W * V1
267 *
268  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
269  $ k, one, v, ldv, work, ldwork )
270  IF( m.GT.k ) THEN
271 *
272 * W := W + C2**H *V2
273 *
274  CALL cgemm( 'Conjugate transpose', 'No transpose', n,
275  $ k, m-k, one, c( k+1, 1 ), ldc,
276  $ v( k+1, 1 ), ldv, one, work, ldwork )
277  END IF
278 *
279 * W := W * T**H or W * T
280 *
281  CALL ctrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
282  $ one, t, ldt, work, ldwork )
283 *
284 * C := C - V * W**H
285 *
286  IF( m.GT.k ) THEN
287 *
288 * C2 := C2 - V2 * W**H
289 *
290  CALL cgemm( 'No transpose', 'Conjugate transpose',
291  $ m-k, n, k, -one, v( k+1, 1 ), ldv, work,
292  $ ldwork, one, c( k+1, 1 ), ldc )
293  END IF
294 *
295 * W := W * V1**H
296 *
297  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
298  $ 'Unit', n, k, one, v, ldv, work, ldwork )
299 *
300 * C1 := C1 - W**H
301 *
302  DO 30 j = 1, k
303  DO 20 i = 1, n
304  c( j, i ) = c( j, i ) - conjg( work( i, j ) )
305  20 CONTINUE
306  30 CONTINUE
307 *
308  ELSE IF( lsame( side, 'R' ) ) THEN
309 *
310 * Form C * H or C * H**H where C = ( C1 C2 )
311 *
312 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
313 *
314 * W := C1
315 *
316  DO 40 j = 1, k
317  CALL ccopy( m, c( 1, j ), 1, work( 1, j ), 1 )
318  40 CONTINUE
319 *
320 * W := W * V1
321 *
322  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
323  $ k, one, v, ldv, work, ldwork )
324  IF( n.GT.k ) THEN
325 *
326 * W := W + C2 * V2
327 *
328  CALL cgemm( 'No transpose', 'No transpose', m, k, n-k,
329  $ one, c( 1, k+1 ), ldc, v( k+1, 1 ), ldv,
330  $ one, work, ldwork )
331  END IF
332 *
333 * W := W * T or W * T**H
334 *
335  CALL ctrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
336  $ one, t, ldt, work, ldwork )
337 *
338 * C := C - W * V**H
339 *
340  IF( n.GT.k ) THEN
341 *
342 * C2 := C2 - W * V2**H
343 *
344  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
345  $ n-k, k, -one, work, ldwork, v( k+1, 1 ),
346  $ ldv, one, c( 1, k+1 ), ldc )
347  END IF
348 *
349 * W := W * V1**H
350 *
351  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
352  $ 'Unit', m, k, one, v, ldv, work, ldwork )
353 *
354 * C1 := C1 - W
355 *
356  DO 60 j = 1, k
357  DO 50 i = 1, m
358  c( i, j ) = c( i, j ) - work( i, j )
359  50 CONTINUE
360  60 CONTINUE
361  END IF
362 *
363  ELSE
364 *
365 * Let V = ( V1 )
366 * ( V2 ) (last K rows)
367 * where V2 is unit upper triangular.
368 *
369  IF( lsame( side, 'L' ) ) THEN
370 *
371 * Form H * C or H**H * C where C = ( C1 )
372 * ( C2 )
373 *
374 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
375 *
376 * W := C2**H
377 *
378  DO 70 j = 1, k
379  CALL ccopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
380  CALL clacgv( n, work( 1, j ), 1 )
381  70 CONTINUE
382 *
383 * W := W * V2
384 *
385  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
386  $ k, one, v( m-k+1, 1 ), ldv, work, ldwork )
387  IF( m.GT.k ) THEN
388 *
389 * W := W + C1**H * V1
390 *
391  CALL cgemm( 'Conjugate transpose', 'No transpose', n,
392  $ k, m-k, one, c, ldc, v, ldv, one, work,
393  $ ldwork )
394  END IF
395 *
396 * W := W * T**H or W * T
397 *
398  CALL ctrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
399  $ one, t, ldt, work, ldwork )
400 *
401 * C := C - V * W**H
402 *
403  IF( m.GT.k ) THEN
404 *
405 * C1 := C1 - V1 * W**H
406 *
407  CALL cgemm( 'No transpose', 'Conjugate transpose',
408  $ m-k, n, k, -one, v, ldv, work, ldwork,
409  $ one, c, ldc )
410  END IF
411 *
412 * W := W * V2**H
413 *
414  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
415  $ 'Unit', n, k, one, v( m-k+1, 1 ), ldv, work,
416  $ ldwork )
417 *
418 * C2 := C2 - W**H
419 *
420  DO 90 j = 1, k
421  DO 80 i = 1, n
422  c( m-k+j, i ) = c( m-k+j, i ) -
423  $ conjg( work( i, j ) )
424  80 CONTINUE
425  90 CONTINUE
426 *
427  ELSE IF( lsame( side, 'R' ) ) THEN
428 *
429 * Form C * H or C * H**H where C = ( C1 C2 )
430 *
431 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
432 *
433 * W := C2
434 *
435  DO 100 j = 1, k
436  CALL ccopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
437  100 CONTINUE
438 *
439 * W := W * V2
440 *
441  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
442  $ k, one, v( n-k+1, 1 ), ldv, work, ldwork )
443  IF( n.GT.k ) THEN
444 *
445 * W := W + C1 * V1
446 *
447  CALL cgemm( 'No transpose', 'No transpose', m, k, n-k,
448  $ one, c, ldc, v, ldv, one, work, ldwork )
449  END IF
450 *
451 * W := W * T or W * T**H
452 *
453  CALL ctrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
454  $ one, t, ldt, work, ldwork )
455 *
456 * C := C - W * V**H
457 *
458  IF( n.GT.k ) THEN
459 *
460 * C1 := C1 - W * V1**H
461 *
462  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
463  $ n-k, k, -one, work, ldwork, v, ldv, one,
464  $ c, ldc )
465  END IF
466 *
467 * W := W * V2**H
468 *
469  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
470  $ 'Unit', m, k, one, v( n-k+1, 1 ), ldv, work,
471  $ ldwork )
472 *
473 * C2 := C2 - W
474 *
475  DO 120 j = 1, k
476  DO 110 i = 1, m
477  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
478  110 CONTINUE
479  120 CONTINUE
480  END IF
481  END IF
482 *
483  ELSE IF( lsame( storev, 'R' ) ) THEN
484 *
485  IF( lsame( direct, 'F' ) ) THEN
486 *
487 * Let V = ( V1 V2 ) (V1: first K columns)
488 * where V1 is unit upper triangular.
489 *
490  IF( lsame( side, 'L' ) ) THEN
491 *
492 * Form H * C or H**H * C where C = ( C1 )
493 * ( C2 )
494 *
495 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
496 *
497 * W := C1**H
498 *
499  DO 130 j = 1, k
500  CALL ccopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
501  CALL clacgv( n, work( 1, j ), 1 )
502  130 CONTINUE
503 *
504 * W := W * V1**H
505 *
506  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
507  $ 'Unit', n, k, one, v, ldv, work, ldwork )
508  IF( m.GT.k ) THEN
509 *
510 * W := W + C2**H * V2**H
511 *
512  CALL cgemm( 'Conjugate transpose',
513  $ 'Conjugate transpose', n, k, m-k, one,
514  $ c( k+1, 1 ), ldc, v( 1, k+1 ), ldv, one,
515  $ work, ldwork )
516  END IF
517 *
518 * W := W * T**H or W * T
519 *
520  CALL ctrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
521  $ one, t, ldt, work, ldwork )
522 *
523 * C := C - V**H * W**H
524 *
525  IF( m.GT.k ) THEN
526 *
527 * C2 := C2 - V2**H * W**H
528 *
529  CALL cgemm( 'Conjugate transpose',
530  $ 'Conjugate transpose', m-k, n, k, -one,
531  $ v( 1, k+1 ), ldv, work, ldwork, one,
532  $ c( k+1, 1 ), ldc )
533  END IF
534 *
535 * W := W * V1
536 *
537  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
538  $ k, one, v, ldv, work, ldwork )
539 *
540 * C1 := C1 - W**H
541 *
542  DO 150 j = 1, k
543  DO 140 i = 1, n
544  c( j, i ) = c( j, i ) - conjg( work( i, j ) )
545  140 CONTINUE
546  150 CONTINUE
547 *
548  ELSE IF( lsame( side, 'R' ) ) THEN
549 *
550 * Form C * H or C * H**H where C = ( C1 C2 )
551 *
552 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
553 *
554 * W := C1
555 *
556  DO 160 j = 1, k
557  CALL ccopy( m, c( 1, j ), 1, work( 1, j ), 1 )
558  160 CONTINUE
559 *
560 * W := W * V1**H
561 *
562  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
563  $ 'Unit', m, k, one, v, ldv, work, ldwork )
564  IF( n.GT.k ) THEN
565 *
566 * W := W + C2 * V2**H
567 *
568  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
569  $ k, n-k, one, c( 1, k+1 ), ldc,
570  $ v( 1, k+1 ), ldv, one, work, ldwork )
571  END IF
572 *
573 * W := W * T or W * T**H
574 *
575  CALL ctrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
576  $ one, t, ldt, work, ldwork )
577 *
578 * C := C - W * V
579 *
580  IF( n.GT.k ) THEN
581 *
582 * C2 := C2 - W * V2
583 *
584  CALL cgemm( 'No transpose', 'No transpose', m, n-k, k,
585  $ -one, work, ldwork, v( 1, k+1 ), ldv, one,
586  $ c( 1, k+1 ), ldc )
587  END IF
588 *
589 * W := W * V1
590 *
591  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
592  $ k, one, v, ldv, work, ldwork )
593 *
594 * C1 := C1 - W
595 *
596  DO 180 j = 1, k
597  DO 170 i = 1, m
598  c( i, j ) = c( i, j ) - work( i, j )
599  170 CONTINUE
600  180 CONTINUE
601 *
602  END IF
603 *
604  ELSE
605 *
606 * Let V = ( V1 V2 ) (V2: last K columns)
607 * where V2 is unit lower triangular.
608 *
609  IF( lsame( side, 'L' ) ) THEN
610 *
611 * Form H * C or H**H * C where C = ( C1 )
612 * ( C2 )
613 *
614 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
615 *
616 * W := C2**H
617 *
618  DO 190 j = 1, k
619  CALL ccopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
620  CALL clacgv( n, work( 1, j ), 1 )
621  190 CONTINUE
622 *
623 * W := W * V2**H
624 *
625  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
626  $ 'Unit', n, k, one, v( 1, m-k+1 ), ldv, work,
627  $ ldwork )
628  IF( m.GT.k ) THEN
629 *
630 * W := W + C1**H * V1**H
631 *
632  CALL cgemm( 'Conjugate transpose',
633  $ 'Conjugate transpose', n, k, m-k, one, c,
634  $ ldc, v, ldv, one, work, ldwork )
635  END IF
636 *
637 * W := W * T**H or W * T
638 *
639  CALL ctrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
640  $ one, t, ldt, work, ldwork )
641 *
642 * C := C - V**H * W**H
643 *
644  IF( m.GT.k ) THEN
645 *
646 * C1 := C1 - V1**H * W**H
647 *
648  CALL cgemm( 'Conjugate transpose',
649  $ 'Conjugate transpose', m-k, n, k, -one, v,
650  $ ldv, work, ldwork, one, c, ldc )
651  END IF
652 *
653 * W := W * V2
654 *
655  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
656  $ k, one, v( 1, m-k+1 ), ldv, work, ldwork )
657 *
658 * C2 := C2 - W**H
659 *
660  DO 210 j = 1, k
661  DO 200 i = 1, n
662  c( m-k+j, i ) = c( m-k+j, i ) -
663  $ conjg( work( i, j ) )
664  200 CONTINUE
665  210 CONTINUE
666 *
667  ELSE IF( lsame( side, 'R' ) ) THEN
668 *
669 * Form C * H or C * H**H where C = ( C1 C2 )
670 *
671 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
672 *
673 * W := C2
674 *
675  DO 220 j = 1, k
676  CALL ccopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
677  220 CONTINUE
678 *
679 * W := W * V2**H
680 *
681  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
682  $ 'Unit', m, k, one, v( 1, n-k+1 ), ldv, work,
683  $ ldwork )
684  IF( n.GT.k ) THEN
685 *
686 * W := W + C1 * V1**H
687 *
688  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
689  $ k, n-k, one, c, ldc, v, ldv, one, work,
690  $ ldwork )
691  END IF
692 *
693 * W := W * T or W * T**H
694 *
695  CALL ctrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
696  $ one, t, ldt, work, ldwork )
697 *
698 * C := C - W * V
699 *
700  IF( n.GT.k ) THEN
701 *
702 * C1 := C1 - W * V1
703 *
704  CALL cgemm( 'No transpose', 'No transpose', m, n-k, k,
705  $ -one, work, ldwork, v, ldv, one, c, ldc )
706  END IF
707 *
708 * W := W * V2
709 *
710  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
711  $ k, one, v( 1, n-k+1 ), ldv, work, ldwork )
712 *
713 * C1 := C1 - W
714 *
715  DO 240 j = 1, k
716  DO 230 i = 1, m
717  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
718  230 CONTINUE
719  240 CONTINUE
720 *
721  END IF
722 *
723  END IF
724  END IF
725 *
726  RETURN
727 *
728 * End of CLARFB
729 *
730  END
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:177
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
Definition: clarfb.f:197