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