Annotation of rpl/lapack/lapack/zlarfb.f, revision 1.17
1.12 bertrand 1: *> \brief \b ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 ! bertrand 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.17 ! bertrand 9: *> Download ZLARFB + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfb.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfb.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfb.f">
1.9 bertrand 15: *> [TXT]</a>
1.17 ! bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
22: * T, LDT, C, LDC, WORK, LDWORK )
1.17 ! bertrand 23: *
1.9 bertrand 24: * .. Scalar Arguments ..
25: * CHARACTER DIRECT, SIDE, STOREV, TRANS
26: * INTEGER K, LDC, LDT, LDV, LDWORK, M, N
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 C( LDC, * ), T( LDT, * ), V( LDV, * ),
30: * $ WORK( LDWORK, * )
31: * ..
1.17 ! bertrand 32: *
1.9 bertrand 33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZLARFB applies a complex block reflector H or its transpose H**H to a
40: *> complex 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**H from the Left
50: *> = 'R': apply H or H**H from the Right
51: *> \endverbatim
52: *>
53: *> \param[in] TRANS
54: *> \verbatim
55: *> TRANS is CHARACTER*1
56: *> = 'N': apply H (No transpose)
57: *> = 'C': apply H**H (Conjugate 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 COMPLEX*16 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: *> 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDC,N)
131: *> On entry, the M-by-N matrix C.
132: *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
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 COMPLEX*16 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: *
1.17 ! bertrand 157: *> \author Univ. of Tennessee
! 158: *> \author Univ. of California Berkeley
! 159: *> \author Univ. of Colorado Denver
! 160: *> \author NAG Ltd.
1.9 bertrand 161: *
1.14 bertrand 162: *> \date June 2013
1.9 bertrand 163: *
164: *> \ingroup complex16OTHERauxiliary
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: * =====================================================================
1.1 bertrand 195: SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
196: $ T, LDT, C, LDC, WORK, LDWORK )
197: *
1.17 ! bertrand 198: * -- LAPACK auxiliary routine (version 3.7.0) --
1.1 bertrand 199: * -- LAPACK is a software package provided by Univ. of Tennessee, --
200: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.14 bertrand 201: * June 2013
1.1 bertrand 202: *
203: * .. Scalar Arguments ..
204: CHARACTER DIRECT, SIDE, STOREV, TRANS
205: INTEGER K, LDC, LDT, LDV, LDWORK, M, N
206: * ..
207: * .. Array Arguments ..
208: COMPLEX*16 C( LDC, * ), T( LDT, * ), V( LDV, * ),
209: $ WORK( LDWORK, * )
210: * ..
211: *
212: * =====================================================================
213: *
214: * .. Parameters ..
215: COMPLEX*16 ONE
216: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
217: * ..
218: * .. Local Scalars ..
219: CHARACTER TRANST
1.14 bertrand 220: INTEGER I, J
1.1 bertrand 221: * ..
222: * .. External Functions ..
223: LOGICAL LSAME
1.14 bertrand 224: EXTERNAL LSAME
1.1 bertrand 225: * ..
226: * .. External Subroutines ..
227: EXTERNAL ZCOPY, ZGEMM, ZLACGV, ZTRMM
228: * ..
229: * .. Intrinsic Functions ..
230: INTRINSIC DCONJG
231: * ..
232: * .. Executable Statements ..
233: *
234: * Quick return if possible
235: *
236: IF( M.LE.0 .OR. N.LE.0 )
237: $ RETURN
238: *
239: IF( LSAME( TRANS, 'N' ) ) THEN
240: TRANST = 'C'
241: ELSE
242: TRANST = 'N'
243: END IF
244: *
245: IF( LSAME( STOREV, 'C' ) ) THEN
246: *
247: IF( LSAME( DIRECT, 'F' ) ) THEN
248: *
249: * Let V = ( V1 ) (first K rows)
250: * ( V2 )
251: * where V1 is unit lower triangular.
252: *
253: IF( LSAME( SIDE, 'L' ) ) THEN
254: *
1.8 bertrand 255: * Form H * C or H**H * C where C = ( C1 )
256: * ( C2 )
1.1 bertrand 257: *
1.8 bertrand 258: * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
1.1 bertrand 259: *
1.8 bertrand 260: * W := C1**H
1.1 bertrand 261: *
262: DO 10 J = 1, K
1.14 bertrand 263: CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
264: CALL ZLACGV( N, WORK( 1, J ), 1 )
1.1 bertrand 265: 10 CONTINUE
266: *
267: * W := W * V1
268: *
1.14 bertrand 269: CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
270: $ K, ONE, V, LDV, WORK, LDWORK )
271: IF( M.GT.K ) THEN
272: *
273: * W := W + C2**H * V2
274: *
275: CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
276: $ K, M-K, ONE, C( K+1, 1 ), LDC,
277: $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
1.1 bertrand 278: END IF
279: *
1.8 bertrand 280: * W := W * T**H or W * T
1.1 bertrand 281: *
1.14 bertrand 282: CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
283: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 284: *
1.8 bertrand 285: * C := C - V * W**H
1.1 bertrand 286: *
287: IF( M.GT.K ) THEN
288: *
1.8 bertrand 289: * C2 := C2 - V2 * W**H
1.1 bertrand 290: *
291: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
1.14 bertrand 292: $ M-K, N, K, -ONE, V( K+1, 1 ), LDV, WORK,
293: $ LDWORK, ONE, C( K+1, 1 ), LDC )
1.1 bertrand 294: END IF
295: *
1.8 bertrand 296: * W := W * V1**H
1.1 bertrand 297: *
298: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
1.14 bertrand 299: $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
1.1 bertrand 300: *
1.8 bertrand 301: * C1 := C1 - W**H
1.1 bertrand 302: *
303: DO 30 J = 1, K
1.14 bertrand 304: DO 20 I = 1, N
1.1 bertrand 305: C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
306: 20 CONTINUE
307: 30 CONTINUE
308: *
309: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
310: *
1.8 bertrand 311: * Form C * H or C * H**H where C = ( C1 C2 )
1.1 bertrand 312: *
313: * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
314: *
315: * W := C1
316: *
317: DO 40 J = 1, K
1.14 bertrand 318: CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
1.1 bertrand 319: 40 CONTINUE
320: *
321: * W := W * V1
322: *
1.14 bertrand 323: CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
324: $ K, ONE, V, LDV, WORK, LDWORK )
325: IF( N.GT.K ) THEN
1.1 bertrand 326: *
327: * W := W + C2 * V2
328: *
1.14 bertrand 329: CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
330: $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
331: $ ONE, WORK, LDWORK )
1.1 bertrand 332: END IF
333: *
1.8 bertrand 334: * W := W * T or W * T**H
1.1 bertrand 335: *
1.14 bertrand 336: CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
337: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 338: *
1.8 bertrand 339: * C := C - W * V**H
1.1 bertrand 340: *
1.14 bertrand 341: IF( N.GT.K ) THEN
1.1 bertrand 342: *
1.8 bertrand 343: * C2 := C2 - W * V2**H
1.1 bertrand 344: *
1.14 bertrand 345: CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
346: $ N-K, K, -ONE, WORK, LDWORK, V( K+1, 1 ),
347: $ LDV, ONE, C( 1, K+1 ), LDC )
1.1 bertrand 348: END IF
349: *
1.8 bertrand 350: * W := W * V1**H
1.1 bertrand 351: *
352: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
1.14 bertrand 353: $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
1.1 bertrand 354: *
355: * C1 := C1 - W
356: *
357: DO 60 J = 1, K
1.14 bertrand 358: DO 50 I = 1, M
1.1 bertrand 359: C( I, J ) = C( I, J ) - WORK( I, J )
360: 50 CONTINUE
361: 60 CONTINUE
362: END IF
363: *
364: ELSE
365: *
366: * Let V = ( V1 )
367: * ( V2 ) (last K rows)
368: * where V2 is unit upper triangular.
369: *
370: IF( LSAME( SIDE, 'L' ) ) THEN
371: *
1.8 bertrand 372: * Form H * C or H**H * C where C = ( C1 )
373: * ( C2 )
1.1 bertrand 374: *
1.8 bertrand 375: * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
1.1 bertrand 376: *
1.8 bertrand 377: * W := C2**H
1.1 bertrand 378: *
379: DO 70 J = 1, K
1.14 bertrand 380: CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
381: CALL ZLACGV( N, WORK( 1, J ), 1 )
1.1 bertrand 382: 70 CONTINUE
383: *
384: * W := W * V2
385: *
1.14 bertrand 386: CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
387: $ K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
1.12 bertrand 388: IF( M.GT.K ) THEN
1.1 bertrand 389: *
1.14 bertrand 390: * W := W + C1**H * V1
1.1 bertrand 391: *
1.14 bertrand 392: CALL ZGEMM( 'Conjugate transpose', 'No transpose', N,
393: $ K, M-K, ONE, C, LDC, V, LDV, ONE, WORK,
394: $ LDWORK )
1.1 bertrand 395: END IF
396: *
1.8 bertrand 397: * W := W * T**H or W * T
1.1 bertrand 398: *
1.14 bertrand 399: CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
400: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 401: *
1.8 bertrand 402: * C := C - V * W**H
1.1 bertrand 403: *
1.12 bertrand 404: IF( M.GT.K ) THEN
1.1 bertrand 405: *
1.8 bertrand 406: * C1 := C1 - V1 * W**H
1.1 bertrand 407: *
408: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
1.14 bertrand 409: $ M-K, N, K, -ONE, V, LDV, WORK, LDWORK,
410: $ ONE, C, LDC )
1.1 bertrand 411: END IF
412: *
1.8 bertrand 413: * W := W * V2**H
1.1 bertrand 414: *
415: CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
1.14 bertrand 416: $ 'Unit', N, K, ONE, V( M-K+1, 1 ), LDV, WORK,
417: $ LDWORK )
1.1 bertrand 418: *
1.8 bertrand 419: * C2 := C2 - W**H
1.1 bertrand 420: *
421: DO 90 J = 1, K
1.14 bertrand 422: DO 80 I = 1, N
1.12 bertrand 423: C( M-K+J, I ) = C( M-K+J, I ) -
1.1 bertrand 424: $ DCONJG( WORK( I, J ) )
425: 80 CONTINUE
426: 90 CONTINUE
427: *
428: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
429: *
1.8 bertrand 430: * Form C * H or C * H**H where C = ( C1 C2 )
1.1 bertrand 431: *
432: * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
433: *
434: * W := C2
435: *
436: DO 100 J = 1, K
1.14 bertrand 437: CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
1.1 bertrand 438: 100 CONTINUE
439: *
440: * W := W * V2
441: *
1.14 bertrand 442: CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
443: $ K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
1.12 bertrand 444: IF( N.GT.K ) THEN
1.1 bertrand 445: *
446: * W := W + C1 * V1
447: *
1.14 bertrand 448: CALL ZGEMM( 'No transpose', 'No transpose', M, K, N-K,
449: $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
1.1 bertrand 450: END IF
451: *
1.8 bertrand 452: * W := W * T or W * T**H
1.1 bertrand 453: *
1.14 bertrand 454: CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
455: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 456: *
1.8 bertrand 457: * C := C - W * V**H
1.1 bertrand 458: *
1.12 bertrand 459: IF( N.GT.K ) THEN
1.1 bertrand 460: *
1.8 bertrand 461: * C1 := C1 - W * V1**H
1.1 bertrand 462: *
1.14 bertrand 463: CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
464: $ N-K, K, -ONE, WORK, LDWORK, V, LDV, ONE,
465: $ C, LDC )
1.1 bertrand 466: END IF
467: *
1.8 bertrand 468: * W := W * V2**H
1.1 bertrand 469: *
470: CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
1.14 bertrand 471: $ 'Unit', M, K, ONE, V( N-K+1, 1 ), LDV, WORK,
472: $ LDWORK )
1.1 bertrand 473: *
474: * C2 := C2 - W
475: *
476: DO 120 J = 1, K
1.14 bertrand 477: DO 110 I = 1, M
478: C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
1.1 bertrand 479: 110 CONTINUE
480: 120 CONTINUE
481: END IF
482: END IF
483: *
484: ELSE IF( LSAME( STOREV, 'R' ) ) THEN
485: *
486: IF( LSAME( DIRECT, 'F' ) ) THEN
487: *
488: * Let V = ( V1 V2 ) (V1: first K columns)
489: * where V1 is unit upper triangular.
490: *
491: IF( LSAME( SIDE, 'L' ) ) THEN
492: *
1.8 bertrand 493: * Form H * C or H**H * C where C = ( C1 )
494: * ( C2 )
1.1 bertrand 495: *
1.8 bertrand 496: * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
1.1 bertrand 497: *
1.8 bertrand 498: * W := C1**H
1.1 bertrand 499: *
500: DO 130 J = 1, K
1.14 bertrand 501: CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
502: CALL ZLACGV( N, WORK( 1, J ), 1 )
1.1 bertrand 503: 130 CONTINUE
504: *
1.8 bertrand 505: * W := W * V1**H
1.1 bertrand 506: *
507: CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
1.14 bertrand 508: $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
509: IF( M.GT.K ) THEN
1.1 bertrand 510: *
1.14 bertrand 511: * W := W + C2**H * V2**H
1.1 bertrand 512: *
513: CALL ZGEMM( 'Conjugate transpose',
1.14 bertrand 514: $ 'Conjugate transpose', N, K, M-K, ONE,
515: $ C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
516: $ WORK, LDWORK )
1.1 bertrand 517: END IF
518: *
1.8 bertrand 519: * W := W * T**H or W * T
1.1 bertrand 520: *
1.14 bertrand 521: CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
522: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 523: *
1.8 bertrand 524: * C := C - V**H * W**H
1.1 bertrand 525: *
1.14 bertrand 526: IF( M.GT.K ) THEN
1.1 bertrand 527: *
1.8 bertrand 528: * C2 := C2 - V2**H * W**H
1.1 bertrand 529: *
530: CALL ZGEMM( 'Conjugate transpose',
1.14 bertrand 531: $ 'Conjugate transpose', M-K, N, K, -ONE,
532: $ V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
533: $ C( K+1, 1 ), LDC )
1.1 bertrand 534: END IF
535: *
536: * W := W * V1
537: *
1.14 bertrand 538: CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
539: $ K, ONE, V, LDV, WORK, LDWORK )
1.1 bertrand 540: *
1.8 bertrand 541: * C1 := C1 - W**H
1.1 bertrand 542: *
543: DO 150 J = 1, K
1.14 bertrand 544: DO 140 I = 1, N
1.1 bertrand 545: C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
546: 140 CONTINUE
547: 150 CONTINUE
548: *
549: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
550: *
1.8 bertrand 551: * Form C * H or C * H**H where C = ( C1 C2 )
1.1 bertrand 552: *
1.8 bertrand 553: * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
1.1 bertrand 554: *
555: * W := C1
556: *
557: DO 160 J = 1, K
1.14 bertrand 558: CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
1.1 bertrand 559: 160 CONTINUE
560: *
1.8 bertrand 561: * W := W * V1**H
1.1 bertrand 562: *
563: CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
1.14 bertrand 564: $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
565: IF( N.GT.K ) THEN
1.1 bertrand 566: *
1.8 bertrand 567: * W := W + C2 * V2**H
1.1 bertrand 568: *
1.14 bertrand 569: CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
570: $ K, N-K, ONE, C( 1, K+1 ), LDC,
571: $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
1.1 bertrand 572: END IF
573: *
1.8 bertrand 574: * W := W * T or W * T**H
1.1 bertrand 575: *
1.14 bertrand 576: CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
577: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 578: *
579: * C := C - W * V
580: *
1.14 bertrand 581: IF( N.GT.K ) THEN
1.1 bertrand 582: *
583: * C2 := C2 - W * V2
584: *
1.14 bertrand 585: CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
586: $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
587: $ C( 1, K+1 ), LDC )
1.1 bertrand 588: END IF
589: *
590: * W := W * V1
591: *
1.14 bertrand 592: CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
593: $ K, ONE, V, LDV, WORK, LDWORK )
1.1 bertrand 594: *
595: * C1 := C1 - W
596: *
597: DO 180 J = 1, K
1.14 bertrand 598: DO 170 I = 1, M
1.1 bertrand 599: C( I, J ) = C( I, J ) - WORK( I, J )
600: 170 CONTINUE
601: 180 CONTINUE
602: *
603: END IF
604: *
605: ELSE
606: *
607: * Let V = ( V1 V2 ) (V2: last K columns)
608: * where V2 is unit lower triangular.
609: *
610: IF( LSAME( SIDE, 'L' ) ) THEN
611: *
1.8 bertrand 612: * Form H * C or H**H * C where C = ( C1 )
613: * ( C2 )
1.1 bertrand 614: *
1.8 bertrand 615: * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
1.1 bertrand 616: *
1.8 bertrand 617: * W := C2**H
1.1 bertrand 618: *
619: DO 190 J = 1, K
1.14 bertrand 620: CALL ZCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
621: CALL ZLACGV( N, WORK( 1, J ), 1 )
1.1 bertrand 622: 190 CONTINUE
623: *
1.8 bertrand 624: * W := W * V2**H
1.1 bertrand 625: *
626: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
1.14 bertrand 627: $ 'Unit', N, K, ONE, V( 1, M-K+1 ), LDV, WORK,
628: $ LDWORK )
1.12 bertrand 629: IF( M.GT.K ) THEN
1.1 bertrand 630: *
1.8 bertrand 631: * W := W + C1**H * V1**H
1.1 bertrand 632: *
633: CALL ZGEMM( 'Conjugate transpose',
1.14 bertrand 634: $ 'Conjugate transpose', N, K, M-K, ONE, C,
635: $ LDC, V, LDV, ONE, WORK, LDWORK )
1.1 bertrand 636: END IF
637: *
1.8 bertrand 638: * W := W * T**H or W * T
1.1 bertrand 639: *
1.14 bertrand 640: CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
641: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 642: *
1.8 bertrand 643: * C := C - V**H * W**H
1.1 bertrand 644: *
1.12 bertrand 645: IF( M.GT.K ) THEN
1.1 bertrand 646: *
1.8 bertrand 647: * C1 := C1 - V1**H * W**H
1.1 bertrand 648: *
649: CALL ZGEMM( 'Conjugate transpose',
1.14 bertrand 650: $ 'Conjugate transpose', M-K, N, K, -ONE, V,
651: $ LDV, WORK, LDWORK, ONE, C, LDC )
1.1 bertrand 652: END IF
653: *
654: * W := W * V2
655: *
1.14 bertrand 656: CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
657: $ K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
1.1 bertrand 658: *
1.8 bertrand 659: * C2 := C2 - W**H
1.1 bertrand 660: *
661: DO 210 J = 1, K
1.14 bertrand 662: DO 200 I = 1, N
1.12 bertrand 663: C( M-K+J, I ) = C( M-K+J, I ) -
1.1 bertrand 664: $ DCONJG( WORK( I, J ) )
665: 200 CONTINUE
666: 210 CONTINUE
667: *
668: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
669: *
1.8 bertrand 670: * Form C * H or C * H**H where C = ( C1 C2 )
1.1 bertrand 671: *
1.8 bertrand 672: * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
1.1 bertrand 673: *
674: * W := C2
675: *
676: DO 220 J = 1, K
1.14 bertrand 677: CALL ZCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
1.1 bertrand 678: 220 CONTINUE
679: *
1.8 bertrand 680: * W := W * V2**H
1.1 bertrand 681: *
682: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
1.14 bertrand 683: $ 'Unit', M, K, ONE, V( 1, N-K+1 ), LDV, WORK,
684: $ LDWORK )
1.12 bertrand 685: IF( N.GT.K ) THEN
1.1 bertrand 686: *
1.8 bertrand 687: * W := W + C1 * V1**H
1.1 bertrand 688: *
1.14 bertrand 689: CALL ZGEMM( 'No transpose', 'Conjugate transpose', M,
690: $ K, N-K, ONE, C, LDC, V, LDV, ONE, WORK,
691: $ LDWORK )
1.1 bertrand 692: END IF
693: *
1.8 bertrand 694: * W := W * T or W * T**H
1.1 bertrand 695: *
1.14 bertrand 696: CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
697: $ ONE, T, LDT, WORK, LDWORK )
1.1 bertrand 698: *
699: * C := C - W * V
700: *
1.12 bertrand 701: IF( N.GT.K ) THEN
1.1 bertrand 702: *
703: * C1 := C1 - W * V1
704: *
1.14 bertrand 705: CALL ZGEMM( 'No transpose', 'No transpose', M, N-K, K,
706: $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
1.1 bertrand 707: END IF
708: *
709: * W := W * V2
710: *
1.14 bertrand 711: CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
712: $ K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
1.1 bertrand 713: *
714: * C1 := C1 - W
715: *
716: DO 240 J = 1, K
1.14 bertrand 717: DO 230 I = 1, M
1.12 bertrand 718: C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
1.1 bertrand 719: 230 CONTINUE
720: 240 CONTINUE
721: *
722: END IF
723: *
724: END IF
725: END IF
726: *
727: RETURN
728: *
729: * End of ZLARFB
730: *
731: 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