← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dlarfb.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dlarfb.f
 
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 September 2012
 
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.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
      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, LASTV, LASTC, lastv2
 
221
*     ..
 
222
*     .. External Functions ..
 
223
      LOGICAL            LSAME
 
224
      INTEGER            ILADLR, ILADLC
 
225
      EXTERNAL           LSAME, ILADLR, ILADLC
 
226
*     ..
 
227
*     .. External Subroutines ..
 
228
      EXTERNAL           DCOPY, DGEMM, DTRMM
 
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, ILADLR( M, K, V, LDV ) )
 
257
               LASTC = ILADLC( 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 DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 
265
   10          CONTINUE
 
266
*
 
267
*              W := W * V1
 
268
*
 
269
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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, ILADLR( N, K, V, LDV ) )
 
316
               LASTC = ILADLR( 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 DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
 
324
   40          CONTINUE
 
325
*
 
326
*              W := W * V1
 
327
*
 
328
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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 = ILADLC( 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 DCOPY( LASTC, C( M-K+J, 1 ), LDC,
 
390
     $                 WORK( 1, J ), 1 )
 
391
   70          CONTINUE
 
392
*
 
393
*              W := W * V2
 
394
*
 
395
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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 = ILADLR( 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 DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
 
449
  100          CONTINUE
 
450
*
 
451
*              W := W * V2
 
452
*
 
453
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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, ILADLC( K, M, V, LDV ) )
 
510
               LASTC = ILADLC( 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 DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 
518
  130          CONTINUE
 
519
*
 
520
*              W := W * V1**T
 
521
*
 
522
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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, ILADLC( K, N, V, LDV ) )
 
569
               LASTC = ILADLR( 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 DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
 
577
  160          CONTINUE
 
578
*
 
579
*              W := W * V1**T
 
580
*
 
581
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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 = ILADLC( 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 DCOPY( LASTC, C( M-K+J, 1 ), LDC,
 
643
     $                 WORK( 1, J ), 1 )
 
644
  190          CONTINUE
 
645
*
 
646
*              W := W * V2**T
 
647
*
 
648
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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 = ILADLR( 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 DCOPY( LASTC, C( 1, N-K+J ), 1,
 
702
     $                 WORK( 1, J ), 1 )
 
703
  220          CONTINUE
 
704
*
 
705
*              W := W * V2**T
 
706
*
 
707
               CALL DTRMM( '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 DGEMM( '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 DTRMM( '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 DGEMM( '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 DTRMM( '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 DLARFB
 
757
*
 
758
      END