LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slarfb.f
Go to the documentation of this file.
1 *> \brief \b SLARFB applies a block reflector or its 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 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 September 2012
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.2) --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 * September 2012
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  lastc = ilaslc( m, n, c, ldc )
383 *
384 * W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK)
385 *
386 * W := C2**T
387 *
388  DO 70 j = 1, k
389  CALL scopy( lastc, c( m-k+j, 1 ), ldc,
390  $ work( 1, j ), 1 )
391  70 continue
392 *
393 * W := W * V2
394 *
395  CALL strmm( 'Right', 'Upper', 'No transpose', 'Unit',
396  $ lastc, k, one, v( m-k+1, 1 ), ldv,
397  $ work, ldwork )
398  IF( m.GT.k ) THEN
399 *
400 * W := W + C1**T*V1
401 *
402  CALL sgemm( 'Transpose', 'No transpose',
403  $ lastc, k, m-k, one, c, ldc, v, ldv,
404  $ one, work, ldwork )
405  END IF
406 *
407 * W := W * T**T or W * T
408 *
409  CALL strmm( 'Right', 'Lower', transt, 'Non-unit',
410  $ lastc, k, one, t, ldt, work, ldwork )
411 *
412 * C := C - V * W**T
413 *
414  IF( m.GT.k ) THEN
415 *
416 * C1 := C1 - V1 * W**T
417 *
418  CALL sgemm( 'No transpose', 'Transpose',
419  $ m-k, lastc, k, -one, v, ldv, work, ldwork,
420  $ one, c, ldc )
421  END IF
422 *
423 * W := W * V2**T
424 *
425  CALL strmm( 'Right', 'Upper', 'Transpose', 'Unit',
426  $ lastc, k, one, v( m-k+1, 1 ), ldv,
427  $ work, ldwork )
428 *
429 * C2 := C2 - W**T
430 *
431  DO 90 j = 1, k
432  DO 80 i = 1, lastc
433  c( m-k+j, i ) = c( m-k+j, i ) - work(i, j)
434  80 continue
435  90 continue
436 *
437  ELSE IF( lsame( side, 'R' ) ) THEN
438 *
439 * Form C * H or C * H**T where C = ( C1 C2 )
440 *
441  lastc = ilaslr( m, n, c, ldc )
442 *
443 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
444 *
445 * W := C2
446 *
447  DO 100 j = 1, k
448  CALL scopy( lastc, c( 1, n-k+j ), 1, work( 1, j ), 1 )
449  100 continue
450 *
451 * W := W * V2
452 *
453  CALL strmm( 'Right', 'Upper', 'No transpose', 'Unit',
454  $ lastc, k, one, v( n-k+1, 1 ), ldv,
455  $ work, ldwork )
456  IF( n.GT.k ) THEN
457 *
458 * W := W + C1 * V1
459 *
460  CALL sgemm( 'No transpose', 'No transpose',
461  $ lastc, k, n-k, one, c, ldc, v, ldv,
462  $ one, work, ldwork )
463  END IF
464 *
465 * W := W * T or W * T**T
466 *
467  CALL strmm( 'Right', 'Lower', trans, 'Non-unit',
468  $ lastc, k, one, t, ldt, work, ldwork )
469 *
470 * C := C - W * V**T
471 *
472  IF( n.GT.k ) THEN
473 *
474 * C1 := C1 - W * V1**T
475 *
476  CALL sgemm( 'No transpose', 'Transpose',
477  $ lastc, n-k, k, -one, work, ldwork, v, ldv,
478  $ one, c, ldc )
479  END IF
480 *
481 * W := W * V2**T
482 *
483  CALL strmm( 'Right', 'Upper', 'Transpose', 'Unit',
484  $ lastc, k, one, v( n-k+1, 1 ), ldv,
485  $ work, ldwork )
486 *
487 * C2 := C2 - W
488 *
489  DO 120 j = 1, k
490  DO 110 i = 1, lastc
491  c( i, n-k+j ) = c( i, n-k+j ) - work(i, j)
492  110 continue
493  120 continue
494  END IF
495  END IF
496 *
497  ELSE IF( lsame( storev, 'R' ) ) THEN
498 *
499  IF( lsame( direct, 'F' ) ) THEN
500 *
501 * Let V = ( V1 V2 ) (V1: first K columns)
502 * where V1 is unit upper triangular.
503 *
504  IF( lsame( side, 'L' ) ) THEN
505 *
506 * Form H * C or H**T * C where C = ( C1 )
507 * ( C2 )
508 *
509  lastv = max( k, ilaslc( k, m, v, ldv ) )
510  lastc = ilaslc( lastv, n, c, ldc )
511 *
512 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
513 *
514 * W := C1**T
515 *
516  DO 130 j = 1, k
517  CALL scopy( lastc, c( j, 1 ), ldc, work( 1, j ), 1 )
518  130 continue
519 *
520 * W := W * V1**T
521 *
522  CALL strmm( 'Right', 'Upper', 'Transpose', 'Unit',
523  $ lastc, k, one, v, ldv, work, ldwork )
524  IF( lastv.GT.k ) THEN
525 *
526 * W := W + C2**T*V2**T
527 *
528  CALL sgemm( 'Transpose', 'Transpose',
529  $ lastc, k, lastv-k,
530  $ one, c( k+1, 1 ), ldc, v( 1, k+1 ), ldv,
531  $ one, work, ldwork )
532  END IF
533 *
534 * W := W * T**T or W * T
535 *
536  CALL strmm( 'Right', 'Upper', transt, 'Non-unit',
537  $ lastc, k, one, t, ldt, work, ldwork )
538 *
539 * C := C - V**T * W**T
540 *
541  IF( lastv.GT.k ) THEN
542 *
543 * C2 := C2 - V2**T * W**T
544 *
545  CALL sgemm( 'Transpose', 'Transpose',
546  $ lastv-k, lastc, k,
547  $ -one, v( 1, k+1 ), ldv, work, ldwork,
548  $ one, c( k+1, 1 ), ldc )
549  END IF
550 *
551 * W := W * V1
552 *
553  CALL strmm( 'Right', 'Upper', 'No transpose', 'Unit',
554  $ lastc, k, one, v, ldv, work, ldwork )
555 *
556 * C1 := C1 - W**T
557 *
558  DO 150 j = 1, k
559  DO 140 i = 1, lastc
560  c( j, i ) = c( j, i ) - work( i, j )
561  140 continue
562  150 continue
563 *
564  ELSE IF( lsame( side, 'R' ) ) THEN
565 *
566 * Form C * H or C * H**T where C = ( C1 C2 )
567 *
568  lastv = max( k, ilaslc( k, n, v, ldv ) )
569  lastc = ilaslr( m, lastv, c, ldc )
570 *
571 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
572 *
573 * W := C1
574 *
575  DO 160 j = 1, k
576  CALL scopy( lastc, c( 1, j ), 1, work( 1, j ), 1 )
577  160 continue
578 *
579 * W := W * V1**T
580 *
581  CALL strmm( 'Right', 'Upper', 'Transpose', 'Unit',
582  $ lastc, k, one, v, ldv, work, ldwork )
583  IF( lastv.GT.k ) THEN
584 *
585 * W := W + C2 * V2**T
586 *
587  CALL sgemm( 'No transpose', 'Transpose',
588  $ lastc, k, lastv-k,
589  $ one, c( 1, k+1 ), ldc, v( 1, k+1 ), ldv,
590  $ one, work, ldwork )
591  END IF
592 *
593 * W := W * T or W * T**T
594 *
595  CALL strmm( 'Right', 'Upper', trans, 'Non-unit',
596  $ lastc, k, one, t, ldt, work, ldwork )
597 *
598 * C := C - W * V
599 *
600  IF( lastv.GT.k ) THEN
601 *
602 * C2 := C2 - W * V2
603 *
604  CALL sgemm( 'No transpose', 'No transpose',
605  $ lastc, lastv-k, k,
606  $ -one, work, ldwork, v( 1, k+1 ), ldv,
607  $ one, c( 1, k+1 ), ldc )
608  END IF
609 *
610 * W := W * V1
611 *
612  CALL strmm( 'Right', 'Upper', 'No transpose', 'Unit',
613  $ lastc, k, one, v, ldv, work, ldwork )
614 *
615 * C1 := C1 - W
616 *
617  DO 180 j = 1, k
618  DO 170 i = 1, lastc
619  c( i, j ) = c( i, j ) - work( i, j )
620  170 continue
621  180 continue
622 *
623  END IF
624 *
625  ELSE
626 *
627 * Let V = ( V1 V2 ) (V2: last K columns)
628 * where V2 is unit lower triangular.
629 *
630  IF( lsame( side, 'L' ) ) THEN
631 *
632 * Form H * C or H**T * C where C = ( C1 )
633 * ( C2 )
634 *
635  lastc = ilaslc( m, n, c, ldc )
636 *
637 * W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
638 *
639 * W := C2**T
640 *
641  DO 190 j = 1, k
642  CALL scopy( lastc, c( m-k+j, 1 ), ldc,
643  $ work( 1, j ), 1 )
644  190 continue
645 *
646 * W := W * V2**T
647 *
648  CALL strmm( 'Right', 'Lower', 'Transpose', 'Unit',
649  $ lastc, k, one, v( 1, m-k+1 ), ldv,
650  $ work, ldwork )
651  IF( m.GT.k ) THEN
652 *
653 * W := W + C1**T * V1**T
654 *
655  CALL sgemm( 'Transpose', 'Transpose',
656  $ lastc, k, m-k, one, c, ldc, v, ldv,
657  $ one, work, ldwork )
658  END IF
659 *
660 * W := W * T**T or W * T
661 *
662  CALL strmm( 'Right', 'Lower', transt, 'Non-unit',
663  $ lastc, k, one, t, ldt, work, ldwork )
664 *
665 * C := C - V**T * W**T
666 *
667  IF( m.GT.k ) THEN
668 *
669 * C1 := C1 - V1**T * W**T
670 *
671  CALL sgemm( 'Transpose', 'Transpose',
672  $ m-k, lastc, k, -one, v, ldv, work, ldwork,
673  $ one, c, ldc )
674  END IF
675 *
676 * W := W * V2
677 *
678  CALL strmm( 'Right', 'Lower', 'No transpose', 'Unit',
679  $ lastc, k, one, v( 1, m-k+1 ), ldv,
680  $ work, ldwork )
681 *
682 * C2 := C2 - W**T
683 *
684  DO 210 j = 1, k
685  DO 200 i = 1, lastc
686  c( m-k+j, i ) = c( m-k+j, i ) - work(i, j)
687  200 continue
688  210 continue
689 *
690  ELSE IF( lsame( side, 'R' ) ) THEN
691 *
692 * Form C * H or C * H**T where C = ( C1 C2 )
693 *
694  lastc = ilaslr( m, n, c, ldc )
695 *
696 * W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
697 *
698 * W := C2
699 *
700  DO 220 j = 1, k
701  CALL scopy( lastc, c( 1, n-k+j ), 1,
702  $ work( 1, j ), 1 )
703  220 continue
704 *
705 * W := W * V2**T
706 *
707  CALL strmm( 'Right', 'Lower', 'Transpose', 'Unit',
708  $ lastc, k, one, v( 1, n-k+1 ), ldv,
709  $ work, ldwork )
710  IF( n.GT.k ) THEN
711 *
712 * W := W + C1 * V1**T
713 *
714  CALL sgemm( 'No transpose', 'Transpose',
715  $ lastc, k, n-k, one, c, ldc, v, ldv,
716  $ one, work, ldwork )
717  END IF
718 *
719 * W := W * T or W * T**T
720 *
721  CALL strmm( 'Right', 'Lower', trans, 'Non-unit',
722  $ lastc, k, one, t, ldt, work, ldwork )
723 *
724 * C := C - W * V
725 *
726  IF( n.GT.k ) THEN
727 *
728 * C1 := C1 - W * V1
729 *
730  CALL sgemm( 'No transpose', 'No transpose',
731  $ lastc, n-k, k, -one, work, ldwork, v, ldv,
732  $ one, c, ldc )
733  END IF
734 *
735 * W := W * V2
736 *
737  CALL strmm( 'Right', 'Lower', 'No transpose', 'Unit',
738  $ lastc, k, one, v( 1, n-k+1 ), ldv,
739  $ work, ldwork )
740 *
741 * C1 := C1 - W
742 *
743  DO 240 j = 1, k
744  DO 230 i = 1, lastc
745  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
746  230 continue
747  240 continue
748 *
749  END IF
750 *
751  END IF
752  END IF
753 *
754  return
755 *
756 * End of SLARFB
757 *
758  END