Mercurial > octave-nkf

comparison libcruft/lapack/clarfb.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
2 $ T, LDT, C, LDC, WORK, LDWORK )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIRECT, SIDE, STOREV, TRANS
10 INTEGER K, LDC, LDT, LDV, LDWORK, M, N
11 * ..
12 * .. Array Arguments ..
13 COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
14 $ WORK( LDWORK, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CLARFB applies a complex block reflector H or its transpose H' to a
21 * complex M-by-N matrix C, from either the left or the right.
22 *
23 * Arguments
24 * =========
25 *
26 * SIDE (input) CHARACTER*1
27 * = 'L': apply H or H' from the Left
28 * = 'R': apply H or H' from the Right
29 *
30 * TRANS (input) CHARACTER*1
31 * = 'N': apply H (No transpose)
32 * = 'C': apply H' (Conjugate transpose)
33 *
34 * DIRECT (input) CHARACTER*1
35 * Indicates how H is formed from a product of elementary
36 * reflectors
37 * = 'F': H = H(1) H(2) . . . H(k) (Forward)
38 * = 'B': H = H(k) . . . H(2) H(1) (Backward)
39 *
40 * STOREV (input) CHARACTER*1
41 * Indicates how the vectors which define the elementary
42 * reflectors are stored:
43 * = 'C': Columnwise
44 * = 'R': Rowwise
45 *
46 * M (input) INTEGER
47 * The number of rows of the matrix C.
48 *
49 * N (input) INTEGER
50 * The number of columns of the matrix C.
51 *
52 * K (input) INTEGER
53 * The order of the matrix T (= the number of elementary
54 * reflectors whose product defines the block reflector).
55 *
56 * V (input) COMPLEX array, dimension
57 * (LDV,K) if STOREV = 'C'
58 * (LDV,M) if STOREV = 'R' and SIDE = 'L'
59 * (LDV,N) if STOREV = 'R' and SIDE = 'R'
60 * The matrix V. See further details.
61 *
62 * LDV (input) INTEGER
63 * The leading dimension of the array V.
64 * If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
65 * if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
66 * if STOREV = 'R', LDV >= K.
67 *
68 * T (input) COMPLEX array, dimension (LDT,K)
69 * The triangular K-by-K matrix T in the representation of the
70 * block reflector.
71 *
72 * LDT (input) INTEGER
73 * The leading dimension of the array T. LDT >= K.
74 *
75 * C (input/output) COMPLEX array, dimension (LDC,N)
76 * On entry, the M-by-N matrix C.
77 * On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
78 *
79 * LDC (input) INTEGER
80 * The leading dimension of the array C. LDC >= max(1,M).
81 *
82 * WORK (workspace) COMPLEX array, dimension (LDWORK,K)
83 *
84 * LDWORK (input) INTEGER
85 * The leading dimension of the array WORK.
86 * If SIDE = 'L', LDWORK >= max(1,N);
87 * if SIDE = 'R', LDWORK >= max(1,M).
88 *
89 * =====================================================================
90 *
91 * .. Parameters ..
92 COMPLEX ONE
93 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
94 * ..
95 * .. Local Scalars ..
96 CHARACTER TRANST
97 INTEGER I, J
98 * ..
99 * .. External Functions ..
100 LOGICAL LSAME
101 EXTERNAL LSAME
102 * ..
103 * .. External Subroutines ..
104 EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM
105 * ..
106 * .. Intrinsic Functions ..
107 INTRINSIC CONJG
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 = 'C'
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 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
136 *
137 * W := C1'
138 *
139 DO 10 J = 1, K
140 CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
141 CALL CLACGV( N, WORK( 1, J ), 1 )
142 10 CONTINUE
143 *
144 * W := W * V1
145 *
146 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
147 $ K, ONE, V, LDV, WORK, LDWORK )
148 IF( M.GT.K ) THEN
149 *
150 * W := W + C2'*V2
151 *
152 CALL CGEMM( 'Conjugate transpose', 'No transpose', N,
153 $ K, M-K, ONE, C( K+1, 1 ), LDC,
154 $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
155 END IF
156 *
157 * W := W * T' or W * T
158 *
159 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
160 $ ONE, T, LDT, WORK, LDWORK )
161 *
162 * C := C - V * W'
163 *
164 IF( M.GT.K ) THEN
165 *
166 * C2 := C2 - V2 * W'
167 *
168 CALL CGEMM( 'No transpose', 'Conjugate transpose',
169 $ M-K, N, K, -ONE, V( K+1, 1 ), LDV, WORK,
170 $ LDWORK, ONE, C( K+1, 1 ), LDC )
171 END IF
172 *
173 * W := W * V1'
174 *
175 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
176 $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
177 *
178 * C1 := C1 - W'
179 *
180 DO 30 J = 1, K
181 DO 20 I = 1, N
182 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
183 20 CONTINUE
184 30 CONTINUE
185 *
186 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
187 *
188 * Form C * H or C * H' where C = ( C1 C2 )
189 *
190 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
191 *
192 * W := C1
193 *
194 DO 40 J = 1, K
195 CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
196 40 CONTINUE
197 *
198 * W := W * V1
199 *
200 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
201 $ K, ONE, V, LDV, WORK, LDWORK )
202 IF( N.GT.K ) THEN
203 *
204 * W := W + C2 * V2
205 *
206 CALL CGEMM( 'No transpose', 'No transpose', M, K, N-K,
207 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
208 $ ONE, WORK, LDWORK )
209 END IF
210 *
211 * W := W * T or W * T'
212 *
213 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
214 $ ONE, T, LDT, WORK, LDWORK )
215 *
216 * C := C - W * V'
217 *
218 IF( N.GT.K ) THEN
219 *
220 * C2 := C2 - W * V2'
221 *
222 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
223 $ N-K, K, -ONE, WORK, LDWORK, V( K+1, 1 ),
224 $ LDV, ONE, C( 1, K+1 ), LDC )
225 END IF
226 *
227 * W := W * V1'
228 *
229 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
230 $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
231 *
232 * C1 := C1 - W
233 *
234 DO 60 J = 1, K
235 DO 50 I = 1, M
236 C( I, J ) = C( I, J ) - WORK( I, J )
237 50 CONTINUE
238 60 CONTINUE
239 END IF
240 *
241 ELSE
242 *
243 * Let V = ( V1 )
244 * ( V2 ) (last K rows)
245 * where V2 is unit upper triangular.
246 *
247 IF( LSAME( SIDE, 'L' ) ) THEN
248 *
249 * Form H * C or H' * C where C = ( C1 )
250 * ( C2 )
251 *
252 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
253 *
254 * W := C2'
255 *
256 DO 70 J = 1, K
257 CALL CCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
258 CALL CLACGV( N, WORK( 1, J ), 1 )
259 70 CONTINUE
260 *
261 * W := W * V2
262 *
263 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
264 $ K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
265 IF( M.GT.K ) THEN
266 *
267 * W := W + C1'*V1
268 *
269 CALL CGEMM( 'Conjugate transpose', 'No transpose', N,
270 $ K, M-K, ONE, C, LDC, V, LDV, ONE, WORK,
271 $ LDWORK )
272 END IF
273 *
274 * W := W * T' or W * T
275 *
276 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
277 $ ONE, T, LDT, WORK, LDWORK )
278 *
279 * C := C - V * W'
280 *
281 IF( M.GT.K ) THEN
282 *
283 * C1 := C1 - V1 * W'
284 *
285 CALL CGEMM( 'No transpose', 'Conjugate transpose',
286 $ M-K, N, K, -ONE, V, LDV, WORK, LDWORK,
287 $ ONE, C, LDC )
288 END IF
289 *
290 * W := W * V2'
291 *
292 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
293 $ 'Unit', N, K, ONE, V( M-K+1, 1 ), LDV, WORK,
294 $ LDWORK )
295 *
296 * C2 := C2 - W'
297 *
298 DO 90 J = 1, K
299 DO 80 I = 1, N
300 C( M-K+J, I ) = C( M-K+J, I ) -
301 $ CONJG( WORK( I, J ) )
302 80 CONTINUE
303 90 CONTINUE
304 *
305 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
306 *
307 * Form C * H or C * H' where C = ( C1 C2 )
308 *
309 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
310 *
311 * W := C2
312 *
313 DO 100 J = 1, K
314 CALL CCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
315 100 CONTINUE
316 *
317 * W := W * V2
318 *
319 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
320 $ K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
321 IF( N.GT.K ) THEN
322 *
323 * W := W + C1 * V1
324 *
325 CALL CGEMM( 'No transpose', 'No transpose', M, K, N-K,
326 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
327 END IF
328 *
329 * W := W * T or W * T'
330 *
331 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
332 $ ONE, T, LDT, WORK, LDWORK )
333 *
334 * C := C - W * V'
335 *
336 IF( N.GT.K ) THEN
337 *
338 * C1 := C1 - W * V1'
339 *
340 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
341 $ N-K, K, -ONE, WORK, LDWORK, V, LDV, ONE,
342 $ C, LDC )
343 END IF
344 *
345 * W := W * V2'
346 *
347 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
348 $ 'Unit', M, K, ONE, V( N-K+1, 1 ), LDV, WORK,
349 $ LDWORK )
350 *
351 * C2 := C2 - W
352 *
353 DO 120 J = 1, K
354 DO 110 I = 1, M
355 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
356 110 CONTINUE
357 120 CONTINUE
358 END IF
359 END IF
360 *
361 ELSE IF( LSAME( STOREV, 'R' ) ) THEN
362 *
363 IF( LSAME( DIRECT, 'F' ) ) THEN
364 *
365 * Let V = ( V1 V2 ) (V1: first K columns)
366 * where V1 is unit upper triangular.
367 *
368 IF( LSAME( SIDE, 'L' ) ) THEN
369 *
370 * Form H * C or H' * C where C = ( C1 )
371 * ( C2 )
372 *
373 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
374 *
375 * W := C1'
376 *
377 DO 130 J = 1, K
378 CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
379 CALL CLACGV( N, WORK( 1, J ), 1 )
380 130 CONTINUE
381 *
382 * W := W * V1'
383 *
384 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
385 $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
386 IF( M.GT.K ) THEN
387 *
388 * W := W + C2'*V2'
389 *
390 CALL CGEMM( 'Conjugate transpose',
391 $ 'Conjugate transpose', N, K, M-K, ONE,
392 $ C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
393 $ WORK, LDWORK )
394 END IF
395 *
396 * W := W * T' or W * T
397 *
398 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
399 $ ONE, T, LDT, WORK, LDWORK )
400 *
401 * C := C - V' * W'
402 *
403 IF( M.GT.K ) THEN
404 *
405 * C2 := C2 - V2' * W'
406 *
407 CALL CGEMM( 'Conjugate transpose',
408 $ 'Conjugate transpose', M-K, N, K, -ONE,
409 $ V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
410 $ C( K+1, 1 ), LDC )
411 END IF
412 *
413 * W := W * V1
414 *
415 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
416 $ K, ONE, V, LDV, WORK, LDWORK )
417 *
418 * C1 := C1 - W'
419 *
420 DO 150 J = 1, K
421 DO 140 I = 1, N
422 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
423 140 CONTINUE
424 150 CONTINUE
425 *
426 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
427 *
428 * Form C * H or C * H' where C = ( C1 C2 )
429 *
430 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
431 *
432 * W := C1
433 *
434 DO 160 J = 1, K
435 CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
436 160 CONTINUE
437 *
438 * W := W * V1'
439 *
440 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
441 $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
442 IF( N.GT.K ) THEN
443 *
444 * W := W + C2 * V2'
445 *
446 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
447 $ K, N-K, ONE, C( 1, K+1 ), LDC,
448 $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
449 END IF
450 *
451 * W := W * T or W * T'
452 *
453 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
454 $ ONE, T, LDT, WORK, LDWORK )
455 *
456 * C := C - W * V
457 *
458 IF( N.GT.K ) THEN
459 *
460 * C2 := C2 - W * V2
461 *
462 CALL CGEMM( 'No transpose', 'No transpose', M, N-K, K,
463 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
464 $ C( 1, K+1 ), LDC )
465 END IF
466 *
467 * W := W * V1
468 *
469 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
470 $ K, ONE, V, LDV, WORK, LDWORK )
471 *
472 * C1 := C1 - W
473 *
474 DO 180 J = 1, K
475 DO 170 I = 1, M
476 C( I, J ) = C( I, J ) - WORK( I, J )
477 170 CONTINUE
478 180 CONTINUE
479 *
480 END IF
481 *
482 ELSE
483 *
484 * Let V = ( V1 V2 ) (V2: last K columns)
485 * where V2 is unit lower triangular.
486 *
487 IF( LSAME( SIDE, 'L' ) ) THEN
488 *
489 * Form H * C or H' * C where C = ( C1 )
490 * ( C2 )
491 *
492 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
493 *
494 * W := C2'
495 *
496 DO 190 J = 1, K
497 CALL CCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
498 CALL CLACGV( N, WORK( 1, J ), 1 )
499 190 CONTINUE
500 *
501 * W := W * V2'
502 *
503 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
504 $ 'Unit', N, K, ONE, V( 1, M-K+1 ), LDV, WORK,
505 $ LDWORK )
506 IF( M.GT.K ) THEN
507 *
508 * W := W + C1'*V1'
509 *
510 CALL CGEMM( 'Conjugate transpose',
511 $ 'Conjugate transpose', N, K, M-K, ONE, C,
512 $ LDC, V, LDV, ONE, WORK, LDWORK )
513 END IF
514 *
515 * W := W * T' or W * T
516 *
517 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
518 $ ONE, T, LDT, WORK, LDWORK )
519 *
520 * C := C - V' * W'
521 *
522 IF( M.GT.K ) THEN
523 *
524 * C1 := C1 - V1' * W'
525 *
526 CALL CGEMM( 'Conjugate transpose',
527 $ 'Conjugate transpose', M-K, N, K, -ONE, V,
528 $ LDV, WORK, LDWORK, ONE, C, LDC )
529 END IF
530 *
531 * W := W * V2
532 *
533 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
534 $ K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
535 *
536 * C2 := C2 - W'
537 *
538 DO 210 J = 1, K
539 DO 200 I = 1, N
540 C( M-K+J, I ) = C( M-K+J, I ) -
541 $ CONJG( WORK( I, J ) )
542 200 CONTINUE
543 210 CONTINUE
544 *
545 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
546 *
547 * Form C * H or C * H' where C = ( C1 C2 )
548 *
549 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
550 *
551 * W := C2
552 *
553 DO 220 J = 1, K
554 CALL CCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
555 220 CONTINUE
556 *
557 * W := W * V2'
558 *
559 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
560 $ 'Unit', M, K, ONE, V( 1, N-K+1 ), LDV, WORK,
561 $ LDWORK )
562 IF( N.GT.K ) THEN
563 *
564 * W := W + C1 * V1'
565 *
566 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
567 $ K, N-K, ONE, C, LDC, V, LDV, ONE, WORK,
568 $ LDWORK )
569 END IF
570 *
571 * W := W * T or W * T'
572 *
573 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
574 $ ONE, T, LDT, WORK, LDWORK )
575 *
576 * C := C - W * V
577 *
578 IF( N.GT.K ) THEN
579 *
580 * C1 := C1 - W * V1
581 *
582 CALL CGEMM( 'No transpose', 'No transpose', M, N-K, K,
583 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
584 END IF
585 *
586 * W := W * V2
587 *
588 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
589 $ K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
590 *
591 * C1 := C1 - W
592 *
593 DO 240 J = 1, K
594 DO 230 I = 1, M
595 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
596 230 CONTINUE
597 240 CONTINUE
598 *
599 END IF
600 *
601 END IF
602 END IF
603 *
604 RETURN
605 *
606 * End of CLARFB
607 *
608 END