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