Mercurial > hg > octave-nkf

comparison libcruft/lapack/slarfb.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 SLARFB( 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 REAL C( LDC, * ), T( LDT, * ), V( LDV, * ),
14 $ WORK( LDWORK, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * SLARFB applies a real block reflector H or its transpose H' to a
21 * real 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 * = 'T': apply H' (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) REAL 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) REAL 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) REAL 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. LDA >= max(1,M).
81 *
82 * WORK (workspace) REAL 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 REAL ONE
93 PARAMETER ( ONE = 1.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 SCOPY, SGEMM, STRMM
105 * ..
106 * .. Executable Statements ..
107 *
108 * Quick return if possible
109 *
110 IF( M.LE.0 .OR. N.LE.0 )
111 $ RETURN
112 *
113 IF( LSAME( TRANS, 'N' ) ) THEN
114 TRANST = 'T'
115 ELSE
116 TRANST = 'N'
117 END IF
118 *
119 IF( LSAME( STOREV, 'C' ) ) THEN
120 *
121 IF( LSAME( DIRECT, 'F' ) ) THEN
122 *
123 * Let V = ( V1 ) (first K rows)
124 * ( V2 )
125 * where V1 is unit lower triangular.
126 *
127 IF( LSAME( SIDE, 'L' ) ) THEN
128 *
129 * Form H * C or H' * C where C = ( C1 )
130 * ( C2 )
131 *
132 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
133 *
134 * W := C1'
135 *
136 DO 10 J = 1, K
137 CALL SCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
138 10 CONTINUE
139 *
140 * W := W * V1
141 *
142 CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
143 $ K, ONE, V, LDV, WORK, LDWORK )
144 IF( M.GT.K ) THEN
145 *
146 * W := W + C2'*V2
147 *
148 CALL SGEMM( 'Transpose', 'No transpose', N, K, M-K,
149 $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
150 $ ONE, WORK, LDWORK )
151 END IF
152 *
153 * W := W * T' or W * T
154 *
155 CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
156 $ ONE, T, LDT, WORK, LDWORK )
157 *
158 * C := C - V * W'
159 *
160 IF( M.GT.K ) THEN
161 *
162 * C2 := C2 - V2 * W'
163 *
164 CALL SGEMM( 'No transpose', 'Transpose', M-K, N, K,
165 $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
166 $ C( K+1, 1 ), LDC )
167 END IF
168 *
169 * W := W * V1'
170 *
171 CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
172 $ ONE, V, LDV, WORK, LDWORK )
173 *
174 * C1 := C1 - W'
175 *
176 DO 30 J = 1, K
177 DO 20 I = 1, N
178 C( J, I ) = C( J, I ) - WORK( I, J )
179 20 CONTINUE
180 30 CONTINUE
181 *
182 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
183 *
184 * Form C * H or C * H' where C = ( C1 C2 )
185 *
186 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
187 *
188 * W := C1
189 *
190 DO 40 J = 1, K
191 CALL SCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
192 40 CONTINUE
193 *
194 * W := W * V1
195 *
196 CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
197 $ K, ONE, V, LDV, WORK, LDWORK )
198 IF( N.GT.K ) THEN
199 *
200 * W := W + C2 * V2
201 *
202 CALL SGEMM( 'No transpose', 'No transpose', M, K, N-K,
203 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
204 $ ONE, WORK, LDWORK )
205 END IF
206 *
207 * W := W * T or W * T'
208 *
209 CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
210 $ ONE, T, LDT, WORK, LDWORK )
211 *
212 * C := C - W * V'
213 *
214 IF( N.GT.K ) THEN
215 *
216 * C2 := C2 - W * V2'
217 *
218 CALL SGEMM( 'No transpose', 'Transpose', M, N-K, K,
219 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
220 $ C( 1, K+1 ), LDC )
221 END IF
222 *
223 * W := W * V1'
224 *
225 CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
226 $ ONE, V, LDV, WORK, LDWORK )
227 *
228 * C1 := C1 - W
229 *
230 DO 60 J = 1, K
231 DO 50 I = 1, M
232 C( I, J ) = C( I, J ) - WORK( I, J )
233 50 CONTINUE
234 60 CONTINUE
235 END IF
236 *
237 ELSE
238 *
239 * Let V = ( V1 )
240 * ( V2 ) (last K rows)
241 * where V2 is unit upper triangular.
242 *
243 IF( LSAME( SIDE, 'L' ) ) THEN
244 *
245 * Form H * C or H' * C where C = ( C1 )
246 * ( C2 )
247 *
248 * W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK)
249 *
250 * W := C2'
251 *
252 DO 70 J = 1, K
253 CALL SCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
254 70 CONTINUE
255 *
256 * W := W * V2
257 *
258 CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
259 $ K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
260 IF( M.GT.K ) THEN
261 *
262 * W := W + C1'*V1
263 *
264 CALL SGEMM( 'Transpose', 'No transpose', N, K, M-K,
265 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
266 END IF
267 *
268 * W := W * T' or W * T
269 *
270 CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
271 $ ONE, T, LDT, WORK, LDWORK )
272 *
273 * C := C - V * W'
274 *
275 IF( M.GT.K ) THEN
276 *
277 * C1 := C1 - V1 * W'
278 *
279 CALL SGEMM( 'No transpose', 'Transpose', M-K, N, K,
280 $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
281 END IF
282 *
283 * W := W * V2'
284 *
285 CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
286 $ ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
287 *
288 * C2 := C2 - W'
289 *
290 DO 90 J = 1, K
291 DO 80 I = 1, N
292 C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
293 80 CONTINUE
294 90 CONTINUE
295 *
296 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
297 *
298 * Form C * H or C * H' where C = ( C1 C2 )
299 *
300 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
301 *
302 * W := C2
303 *
304 DO 100 J = 1, K
305 CALL SCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
306 100 CONTINUE
307 *
308 * W := W * V2
309 *
310 CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
311 $ K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
312 IF( N.GT.K ) THEN
313 *
314 * W := W + C1 * V1
315 *
316 CALL SGEMM( 'No transpose', 'No transpose', M, K, N-K,
317 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
318 END IF
319 *
320 * W := W * T or W * T'
321 *
322 CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
323 $ ONE, T, LDT, WORK, LDWORK )
324 *
325 * C := C - W * V'
326 *
327 IF( N.GT.K ) THEN
328 *
329 * C1 := C1 - W * V1'
330 *
331 CALL SGEMM( 'No transpose', 'Transpose', M, N-K, K,
332 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
333 END IF
334 *
335 * W := W * V2'
336 *
337 CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
338 $ ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
339 *
340 * C2 := C2 - W
341 *
342 DO 120 J = 1, K
343 DO 110 I = 1, M
344 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
345 110 CONTINUE
346 120 CONTINUE
347 END IF
348 END IF
349 *
350 ELSE IF( LSAME( STOREV, 'R' ) ) THEN
351 *
352 IF( LSAME( DIRECT, 'F' ) ) THEN
353 *
354 * Let V = ( V1 V2 ) (V1: first K columns)
355 * where V1 is unit upper triangular.
356 *
357 IF( LSAME( SIDE, 'L' ) ) THEN
358 *
359 * Form H * C or H' * C where C = ( C1 )
360 * ( C2 )
361 *
362 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
363 *
364 * W := C1'
365 *
366 DO 130 J = 1, K
367 CALL SCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
368 130 CONTINUE
369 *
370 * W := W * V1'
371 *
372 CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
373 $ ONE, V, LDV, WORK, LDWORK )
374 IF( M.GT.K ) THEN
375 *
376 * W := W + C2'*V2'
377 *
378 CALL SGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
379 $ C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
380 $ WORK, LDWORK )
381 END IF
382 *
383 * W := W * T' or W * T
384 *
385 CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
386 $ ONE, T, LDT, WORK, LDWORK )
387 *
388 * C := C - V' * W'
389 *
390 IF( M.GT.K ) THEN
391 *
392 * C2 := C2 - V2' * W'
393 *
394 CALL SGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
395 $ V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
396 $ C( K+1, 1 ), LDC )
397 END IF
398 *
399 * W := W * V1
400 *
401 CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
402 $ K, ONE, V, LDV, WORK, LDWORK )
403 *
404 * C1 := C1 - W'
405 *
406 DO 150 J = 1, K
407 DO 140 I = 1, N
408 C( J, I ) = C( J, I ) - WORK( I, J )
409 140 CONTINUE
410 150 CONTINUE
411 *
412 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
413 *
414 * Form C * H or C * H' where C = ( C1 C2 )
415 *
416 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
417 *
418 * W := C1
419 *
420 DO 160 J = 1, K
421 CALL SCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
422 160 CONTINUE
423 *
424 * W := W * V1'
425 *
426 CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
427 $ ONE, V, LDV, WORK, LDWORK )
428 IF( N.GT.K ) THEN
429 *
430 * W := W + C2 * V2'
431 *
432 CALL SGEMM( 'No transpose', 'Transpose', M, K, N-K,
433 $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
434 $ ONE, WORK, LDWORK )
435 END IF
436 *
437 * W := W * T or W * T'
438 *
439 CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
440 $ ONE, T, LDT, WORK, LDWORK )
441 *
442 * C := C - W * V
443 *
444 IF( N.GT.K ) THEN
445 *
446 * C2 := C2 - W * V2
447 *
448 CALL SGEMM( 'No transpose', 'No transpose', M, N-K, K,
449 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
450 $ C( 1, K+1 ), LDC )
451 END IF
452 *
453 * W := W * V1
454 *
455 CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
456 $ K, ONE, V, LDV, WORK, LDWORK )
457 *
458 * C1 := C1 - W
459 *
460 DO 180 J = 1, K
461 DO 170 I = 1, M
462 C( I, J ) = C( I, J ) - WORK( I, J )
463 170 CONTINUE
464 180 CONTINUE
465 *
466 END IF
467 *
468 ELSE
469 *
470 * Let V = ( V1 V2 ) (V2: last K columns)
471 * where V2 is unit lower triangular.
472 *
473 IF( LSAME( SIDE, 'L' ) ) THEN
474 *
475 * Form H * C or H' * C where C = ( C1 )
476 * ( C2 )
477 *
478 * W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK)
479 *
480 * W := C2'
481 *
482 DO 190 J = 1, K
483 CALL SCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
484 190 CONTINUE
485 *
486 * W := W * V2'
487 *
488 CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
489 $ ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
490 IF( M.GT.K ) THEN
491 *
492 * W := W + C1'*V1'
493 *
494 CALL SGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
495 $ C, LDC, V, LDV, ONE, WORK, LDWORK )
496 END IF
497 *
498 * W := W * T' or W * T
499 *
500 CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
501 $ ONE, T, LDT, WORK, LDWORK )
502 *
503 * C := C - V' * W'
504 *
505 IF( M.GT.K ) THEN
506 *
507 * C1 := C1 - V1' * W'
508 *
509 CALL SGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
510 $ V, LDV, WORK, LDWORK, ONE, C, LDC )
511 END IF
512 *
513 * W := W * V2
514 *
515 CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
516 $ K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
517 *
518 * C2 := C2 - W'
519 *
520 DO 210 J = 1, K
521 DO 200 I = 1, N
522 C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
523 200 CONTINUE
524 210 CONTINUE
525 *
526 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
527 *
528 * Form C * H or C * H' where C = ( C1 C2 )
529 *
530 * W := C * V' = (C1*V1' + C2*V2') (stored in WORK)
531 *
532 * W := C2
533 *
534 DO 220 J = 1, K
535 CALL SCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
536 220 CONTINUE
537 *
538 * W := W * V2'
539 *
540 CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
541 $ ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
542 IF( N.GT.K ) THEN
543 *
544 * W := W + C1 * V1'
545 *
546 CALL SGEMM( 'No transpose', 'Transpose', M, K, N-K,
547 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
548 END IF
549 *
550 * W := W * T or W * T'
551 *
552 CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
553 $ ONE, T, LDT, WORK, LDWORK )
554 *
555 * C := C - W * V
556 *
557 IF( N.GT.K ) THEN
558 *
559 * C1 := C1 - W * V1
560 *
561 CALL SGEMM( 'No transpose', 'No transpose', M, N-K, K,
562 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
563 END IF
564 *
565 * W := W * V2
566 *
567 CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
568 $ K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
569 *
570 * C1 := C1 - W
571 *
572 DO 240 J = 1, K
573 DO 230 I = 1, M
574 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
575 230 CONTINUE
576 240 CONTINUE
577 *
578 END IF
579 *
580 END IF
581 END IF
582 *
583 RETURN
584 *
585 * End of SLARFB
586 *
587 END