KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlarfb.f
Go to the documentation of this file.
1  SUBROUTINE dlarfb( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
2  $ T, LDT, C, LDC, WORK, LDWORK )
3 *
4 * -- LAPACK auxiliary routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * February 29, 1992
8 *
9 * .. Scalar Arguments ..
10  CHARACTER DIRECT, SIDE, STOREV, TRANS
11  INTEGER K, LDC, LDT, LDV, LDWORK, M, N
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ),
15  $ work( ldwork, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLARFB applies a real block reflector H or its transpose H' to a
22 * real m by n matrix C, from either the left or the right.
23 *
24 * Arguments
25 * =========
26 *
27 * SIDE (input) CHARACTER*1
28 * = 'L': apply H or H' from the Left
29 * = 'R': apply H or H' from the Right
30 *
31 * TRANS (input) CHARACTER*1
32 * = 'N': apply H (No transpose)
33 * = 'T': apply H' (Transpose)
34 *
35 * DIRECT (input) CHARACTER*1
36 * Indicates how H is formed from a product of elementary
37 * reflectors
38 * = 'F': H = H(1) H(2) . . . H(k) (Forward)
39 * = 'B': H = H(k) . . . H(2) H(1) (Backward)
40 *
41 * STOREV (input) CHARACTER*1
42 * Indicates how the vectors which define the elementary
43 * reflectors are stored:
44 * = 'C': Columnwise
45 * = 'R': Rowwise
46 *
47 * M (input) INTEGER
48 * The number of rows of the matrix C.
49 *
50 * N (input) INTEGER
51 * The number of columns of the matrix C.
52 *
53 * K (input) INTEGER
54 * The order of the matrix T (= the number of elementary
55 * reflectors whose product defines the block reflector).
56 *
57 * V (input) DOUBLE PRECISION array, dimension
58 * (LDV,K) if STOREV = 'C'
59 * (LDV,M) if STOREV = 'R' and SIDE = 'L'
60 * (LDV,N) if STOREV = 'R' and SIDE = 'R'
61 * The matrix V. See further details.
62 *
63 * LDV (input) INTEGER
64 * The leading dimension of the array V.
65 * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
66 * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
67 * if STOREV = 'R', LDV >= K.
68 *
69 * T (input) DOUBLE PRECISION array, dimension (LDT,K)
70 * The triangular k by k matrix T in the representation of the
71 * block reflector.
72 *
73 * LDT (input) INTEGER
74 * The leading dimension of the array T. LDT >= K.
75 *
76 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
77 * On entry, the m by n matrix C.
78 * On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
79 *
80 * LDC (input) INTEGER
81 * The leading dimension of the array C. LDA >= max(1,M).
82 *
83 * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
84 *
85 * LDWORK (input) INTEGER
86 * The leading dimension of the array WORK.
87 * If SIDE = 'L', LDWORK >= max(1,N);
88 * if SIDE = 'R', LDWORK >= max(1,M).
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93  DOUBLE PRECISION ONE
94  parameter( one = 1.0d+0 )
95 * ..
96 * .. Local Scalars ..
97  CHARACTER TRANST
98  INTEGER I, J
99 * ..
100 * .. External Functions ..
101  LOGICAL LSAME
102  EXTERNAL lsame
103 * ..
104 * .. External Subroutines ..
105  EXTERNAL dcopy, dgemm, dtrmm
106 * ..
107 * .. Executable Statements ..
108 *
109 * Quick return if possible
110 *
111  IF( m.LE.0 .OR. n.LE.0 )
112  $ RETURN
113 *
114  IF( lsame( trans, 'N' ) ) THEN
115  transt = 'T'
116  ELSE
117  transt = 'N'
118  END IF
119 *
120  IF( lsame( storev, 'C' ) ) THEN
121 *
122  IF( lsame( direct, 'F' ) ) THEN
123 *
124 * Let V = ( V1 ) (first K rows)
125 * ( V2 )
126 * where V1 is unit lower triangular.
127 *
128  IF( lsame( side, 'L' ) ) THEN
129 *
130 * Form H * C or H' * C where C = ( C1 )
131 * ( C2 )
132 *
133 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
134 *
135 * W := C1'
136 *
137  DO 10 j = 1, k
138  CALL dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
139  10 CONTINUE
140 *
141 * W := W * V1
142 *
143  CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
144  $ k, one, v, ldv, work, ldwork )
145  IF( m.GT.k ) THEN
146 *
147 * W := W + C2'*V2
148 *
149  CALL dgemm( 'Transpose', 'No transpose', n, k, m-k,
150  $ one, c( k+1, 1 ), ldc, v( k+1, 1 ), ldv,
151  $ one, work, ldwork )
152  END IF
153 *
154 * W := W * T' or W * T
155 *
156  CALL dtrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
157  $ one, t, ldt, work, ldwork )
158 *
159 * C := C - V * W'
160 *
161  IF( m.GT.k ) THEN
162 *
163 * C2 := C2 - V2 * W'
164 *
165  CALL dgemm( 'No transpose', 'Transpose', m-k, n, k,
166  $ -one, v( k+1, 1 ), ldv, work, ldwork, one,
167  $ c( k+1, 1 ), ldc )
168  END IF
169 *
170 * W := W * V1'
171 *
172  CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', n, k,
173  $ one, v, ldv, work, ldwork )
174 *
175 * C1 := C1 - W'
176 *
177  DO 30 j = 1, k
178  DO 20 i = 1, n
179  c( j, i ) = c( j, i ) - work( i, j )
180  20 CONTINUE
181  30 CONTINUE
182 *
183  ELSE IF( lsame( side, 'R' ) ) THEN
184 *
185 * Form C * H or C * H' where C = ( C1 C2 )
186 *
187 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
188 *
189 * W := C1
190 *
191  DO 40 j = 1, k
192  CALL dcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
193  40 CONTINUE
194 *
195 * W := W * V1
196 *
197  CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
198  $ k, one, v, ldv, work, ldwork )
199  IF( n.GT.k ) THEN
200 *
201 * W := W + C2 * V2
202 *
203  CALL dgemm( 'No transpose', 'No transpose', m, k, n-k,
204  $ one, c( 1, k+1 ), ldc, v( k+1, 1 ), ldv,
205  $ one, work, ldwork )
206  END IF
207 *
208 * W := W * T or W * T'
209 *
210  CALL dtrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
211  $ one, t, ldt, work, ldwork )
212 *
213 * C := C - W * V'
214 *
215  IF( n.GT.k ) THEN
216 *
217 * C2 := C2 - W * V2'
218 *
219  CALL dgemm( 'No transpose', 'Transpose', m, n-k, k,
220  $ -one, work, ldwork, v( k+1, 1 ), ldv, one,
221  $ c( 1, k+1 ), ldc )
222  END IF
223 *
224 * W := W * V1'
225 *
226  CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', m, k,
227  $ one, v, ldv, work, ldwork )
228 *
229 * C1 := C1 - W
230 *
231  DO 60 j = 1, k
232  DO 50 i = 1, m
233  c( i, j ) = c( i, j ) - work( i, j )
234  50 CONTINUE
235  60 CONTINUE
236  END IF
237 *
238  ELSE
239 *
240 * Let V = ( V1 )
241 * ( V2 ) (last K rows)
242 * where V2 is unit upper triangular.
243 *
244  IF( lsame( side, 'L' ) ) THEN
245 *
246 * Form H * C or H' * C where C = ( C1 )
247 * ( C2 )
248 *
249 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
250 *
251 * W := C2'
252 *
253  DO 70 j = 1, k
254  CALL dcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
255  70 CONTINUE
256 *
257 * W := W * V2
258 *
259  CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
260  $ k, one, v( m-k+1, 1 ), ldv, work, ldwork )
261  IF( m.GT.k ) THEN
262 *
263 * W := W + C1'*V1
264 *
265  CALL dgemm( 'Transpose', 'No transpose', n, k, m-k,
266  $ one, c, ldc, v, ldv, one, work, ldwork )
267  END IF
268 *
269 * W := W * T' or W * T
270 *
271  CALL dtrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
272  $ one, t, ldt, work, ldwork )
273 *
274 * C := C - V * W'
275 *
276  IF( m.GT.k ) THEN
277 *
278 * C1 := C1 - V1 * W'
279 *
280  CALL dgemm( 'No transpose', 'Transpose', m-k, n, k,
281  $ -one, v, ldv, work, ldwork, one, c, ldc )
282  END IF
283 *
284 * W := W * V2'
285 *
286  CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', n, k,
287  $ one, v( m-k+1, 1 ), ldv, work, ldwork )
288 *
289 * C2 := C2 - W'
290 *
291  DO 90 j = 1, k
292  DO 80 i = 1, n
293  c( m-k+j, i ) = c( m-k+j, i ) - work( i, j )
294  80 CONTINUE
295  90 CONTINUE
296 *
297  ELSE IF( lsame( side, 'R' ) ) THEN
298 *
299 * Form C * H or C * H' where C = ( C1 C2 )
300 *
301 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
302 *
303 * W := C2
304 *
305  DO 100 j = 1, k
306  CALL dcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
307  100 CONTINUE
308 *
309 * W := W * V2
310 *
311  CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
312  $ k, one, v( n-k+1, 1 ), ldv, work, ldwork )
313  IF( n.GT.k ) THEN
314 *
315 * W := W + C1 * V1
316 *
317  CALL dgemm( 'No transpose', 'No transpose', m, k, n-k,
318  $ one, c, ldc, v, ldv, one, work, ldwork )
319  END IF
320 *
321 * W := W * T or W * T'
322 *
323  CALL dtrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
324  $ one, t, ldt, work, ldwork )
325 *
326 * C := C - W * V'
327 *
328  IF( n.GT.k ) THEN
329 *
330 * C1 := C1 - W * V1'
331 *
332  CALL dgemm( 'No transpose', 'Transpose', m, n-k, k,
333  $ -one, work, ldwork, v, ldv, one, c, ldc )
334  END IF
335 *
336 * W := W * V2'
337 *
338  CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', m, k,
339  $ one, v( n-k+1, 1 ), ldv, work, ldwork )
340 *
341 * C2 := C2 - W
342 *
343  DO 120 j = 1, k
344  DO 110 i = 1, m
345  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
346  110 CONTINUE
347  120 CONTINUE
348  END IF
349  END IF
350 *
351  ELSE IF( lsame( storev, 'R' ) ) THEN
352 *
353  IF( lsame( direct, 'F' ) ) THEN
354 *
355 * Let V = ( V1 V2 ) (V1: first K columns)
356 * where V1 is unit upper triangular.
357 *
358  IF( lsame( side, 'L' ) ) THEN
359 *
360 * Form H * C or H' * C where C = ( C1 )
361 * ( C2 )
362 *
363 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
364 *
365 * W := C1'
366 *
367  DO 130 j = 1, k
368  CALL dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
369  130 CONTINUE
370 *
371 * W := W * V1'
372 *
373  CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', n, k,
374  $ one, v, ldv, work, ldwork )
375  IF( m.GT.k ) THEN
376 *
377 * W := W + C2'*V2'
378 *
379  CALL dgemm( 'Transpose', 'Transpose', n, k, m-k, one,
380  $ c( k+1, 1 ), ldc, v( 1, k+1 ), ldv, one,
381  $ work, ldwork )
382  END IF
383 *
384 * W := W * T' or W * T
385 *
386  CALL dtrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
387  $ one, t, ldt, work, ldwork )
388 *
389 * C := C - V' * W'
390 *
391  IF( m.GT.k ) THEN
392 *
393 * C2 := C2 - V2' * W'
394 *
395  CALL dgemm( 'Transpose', 'Transpose', m-k, n, k, -one,
396  $ v( 1, k+1 ), ldv, work, ldwork, one,
397  $ c( k+1, 1 ), ldc )
398  END IF
399 *
400 * W := W * V1
401 *
402  CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
403  $ k, one, v, ldv, work, ldwork )
404 *
405 * C1 := C1 - W'
406 *
407  DO 150 j = 1, k
408  DO 140 i = 1, n
409  c( j, i ) = c( j, i ) - work( i, j )
410  140 CONTINUE
411  150 CONTINUE
412 *
413  ELSE IF( lsame( side, 'R' ) ) THEN
414 *
415 * Form C * H or C * H' where C = ( C1 C2 )
416 *
417 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
418 *
419 * W := C1
420 *
421  DO 160 j = 1, k
422  CALL dcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
423  160 CONTINUE
424 *
425 * W := W * V1'
426 *
427  CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', m, k,
428  $ one, v, ldv, work, ldwork )
429  IF( n.GT.k ) THEN
430 *
431 * W := W + C2 * V2'
432 *
433  CALL dgemm( 'No transpose', 'Transpose', m, k, n-k,
434  $ one, c( 1, k+1 ), ldc, v( 1, k+1 ), ldv,
435  $ one, work, ldwork )
436  END IF
437 *
438 * W := W * T or W * T'
439 *
440  CALL dtrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
441  $ one, t, ldt, work, ldwork )
442 *
443 * C := C - W * V
444 *
445  IF( n.GT.k ) THEN
446 *
447 * C2 := C2 - W * V2
448 *
449  CALL dgemm( 'No transpose', 'No transpose', m, n-k, k,
450  $ -one, work, ldwork, v( 1, k+1 ), ldv, one,
451  $ c( 1, k+1 ), ldc )
452  END IF
453 *
454 * W := W * V1
455 *
456  CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
457  $ k, one, v, ldv, work, ldwork )
458 *
459 * C1 := C1 - W
460 *
461  DO 180 j = 1, k
462  DO 170 i = 1, m
463  c( i, j ) = c( i, j ) - work( i, j )
464  170 CONTINUE
465  180 CONTINUE
466 *
467  END IF
468 *
469  ELSE
470 *
471 * Let V = ( V1 V2 ) (V2: last K columns)
472 * where V2 is unit lower triangular.
473 *
474  IF( lsame( side, 'L' ) ) THEN
475 *
476 * Form H * C or H' * C where C = ( C1 )
477 * ( C2 )
478 *
479 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
480 *
481 * W := C2'
482 *
483  DO 190 j = 1, k
484  CALL dcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
485  190 CONTINUE
486 *
487 * W := W * V2'
488 *
489  CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', n, k,
490  $ one, v( 1, m-k+1 ), ldv, work, ldwork )
491  IF( m.GT.k ) THEN
492 *
493 * W := W + C1'*V1'
494 *
495  CALL dgemm( 'Transpose', 'Transpose', n, k, m-k, one,
496  $ c, ldc, v, ldv, one, work, ldwork )
497  END IF
498 *
499 * W := W * T' or W * T
500 *
501  CALL dtrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
502  $ one, t, ldt, work, ldwork )
503 *
504 * C := C - V' * W'
505 *
506  IF( m.GT.k ) THEN
507 *
508 * C1 := C1 - V1' * W'
509 *
510  CALL dgemm( 'Transpose', 'Transpose', m-k, n, k, -one,
511  $ v, ldv, work, ldwork, one, c, ldc )
512  END IF
513 *
514 * W := W * V2
515 *
516  CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
517  $ k, one, v( 1, m-k+1 ), ldv, work, ldwork )
518 *
519 * C2 := C2 - W'
520 *
521  DO 210 j = 1, k
522  DO 200 i = 1, n
523  c( m-k+j, i ) = c( m-k+j, i ) - work( i, j )
524  200 CONTINUE
525  210 CONTINUE
526 *
527  ELSE IF( lsame( side, 'R' ) ) THEN
528 *
529 * Form C * H or C * H' where C = ( C1 C2 )
530 *
531 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
532 *
533 * W := C2
534 *
535  DO 220 j = 1, k
536  CALL dcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
537  220 CONTINUE
538 *
539 * W := W * V2'
540 *
541  CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', m, k,
542  $ one, v( 1, n-k+1 ), ldv, work, ldwork )
543  IF( n.GT.k ) THEN
544 *
545 * W := W + C1 * V1'
546 *
547  CALL dgemm( 'No transpose', 'Transpose', m, k, n-k,
548  $ one, c, ldc, v, ldv, one, work, ldwork )
549  END IF
550 *
551 * W := W * T or W * T'
552 *
553  CALL dtrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
554  $ one, t, ldt, work, ldwork )
555 *
556 * C := C - W * V
557 *
558  IF( n.GT.k ) THEN
559 *
560 * C1 := C1 - W * V1
561 *
562  CALL dgemm( 'No transpose', 'No transpose', m, n-k, k,
563  $ -one, work, ldwork, v, ldv, one, c, ldc )
564  END IF
565 *
566 * W := W * V2
567 *
568  CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
569  $ k, one, v( 1, n-k+1 ), ldv, work, ldwork )
570 *
571 * C1 := C1 - W
572 *
573  DO 240 j = 1, k
574  DO 230 i = 1, m
575  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
576  230 CONTINUE
577  240 CONTINUE
578 *
579  END IF
580 *
581  END IF
582  END IF
583 *
584  RETURN
585 *
586 * End of DLARFB
587 *
588  END
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
Definition: dcopy.f:2
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgemm.f:3
dlarfb
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
Definition: dlarfb.f:3
dtrmm
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrmm.f:3