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