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