comparison libcruft/lapack/dlarz.f @ 6936:e92bc778c634

[project @ 2007-09-30 21:39:05 by dbateman]
author dbateman
date Sun, 30 Sep 2007 21:41:04 +0000
parents
children
comparison
equal deleted inserted replaced
6935:5cd272497aae 6936:e92bc778c634
1 SUBROUTINE DLARZ( SIDE, M, N, L, V, INCV, TAU, C, LDC, WORK )
2 *
3 * -- LAPACK routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER SIDE
9 INTEGER INCV, L, LDC, M, N
10 DOUBLE PRECISION TAU
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLARZ applies a real elementary reflector H to a real M-by-N
20 * matrix C, from either the left or the right. H is represented in the
21 * form
22 *
23 * H = I - tau * v * v'
24 *
25 * where tau is a real scalar and v is a real vector.
26 *
27 * If tau = 0, then H is taken to be the unit matrix.
28 *
29 *
30 * H is a product of k elementary reflectors as returned by DTZRZF.
31 *
32 * Arguments
33 * =========
34 *
35 * SIDE (input) CHARACTER*1
36 * = 'L': form H * C
37 * = 'R': form C * H
38 *
39 * M (input) INTEGER
40 * The number of rows of the matrix C.
41 *
42 * N (input) INTEGER
43 * The number of columns of the matrix C.
44 *
45 * L (input) INTEGER
46 * The number of entries of the vector V containing
47 * the meaningful part of the Householder vectors.
48 * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
49 *
50 * V (input) DOUBLE PRECISION array, dimension (1+(L-1)*abs(INCV))
51 * The vector v in the representation of H as returned by
52 * DTZRZF. V is not used if TAU = 0.
53 *
54 * INCV (input) INTEGER
55 * The increment between elements of v. INCV <> 0.
56 *
57 * TAU (input) DOUBLE PRECISION
58 * The value tau in the representation of H.
59 *
60 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
61 * On entry, the M-by-N matrix C.
62 * On exit, C is overwritten by the matrix H * C if SIDE = 'L',
63 * or C * H if SIDE = 'R'.
64 *
65 * LDC (input) INTEGER
66 * The leading dimension of the array C. LDC >= max(1,M).
67 *
68 * WORK (workspace) DOUBLE PRECISION array, dimension
69 * (N) if SIDE = 'L'
70 * or (M) if SIDE = 'R'
71 *
72 * Further Details
73 * ===============
74 *
75 * Based on contributions by
76 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
77 *
78 * =====================================================================
79 *
80 * .. Parameters ..
81 DOUBLE PRECISION ONE, ZERO
82 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
83 * ..
84 * .. External Subroutines ..
85 EXTERNAL DAXPY, DCOPY, DGEMV, DGER
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 EXTERNAL LSAME
90 * ..
91 * .. Executable Statements ..
92 *
93 IF( LSAME( SIDE, 'L' ) ) THEN
94 *
95 * Form H * C
96 *
97 IF( TAU.NE.ZERO ) THEN
98 *
99 * w( 1:n ) = C( 1, 1:n )
100 *
101 CALL DCOPY( N, C, LDC, WORK, 1 )
102 *
103 * w( 1:n ) = w( 1:n ) + C( m-l+1:m, 1:n )' * v( 1:l )
104 *
105 CALL DGEMV( 'Transpose', L, N, ONE, C( M-L+1, 1 ), LDC, V,
106 $ INCV, ONE, WORK, 1 )
107 *
108 * C( 1, 1:n ) = C( 1, 1:n ) - tau * w( 1:n )
109 *
110 CALL DAXPY( N, -TAU, WORK, 1, C, LDC )
111 *
112 * C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
113 * tau * v( 1:l ) * w( 1:n )'
114 *
115 CALL DGER( L, N, -TAU, V, INCV, WORK, 1, C( M-L+1, 1 ),
116 $ LDC )
117 END IF
118 *
119 ELSE
120 *
121 * Form C * H
122 *
123 IF( TAU.NE.ZERO ) THEN
124 *
125 * w( 1:m ) = C( 1:m, 1 )
126 *
127 CALL DCOPY( M, C, 1, WORK, 1 )
128 *
129 * w( 1:m ) = w( 1:m ) + C( 1:m, n-l+1:n, 1:n ) * v( 1:l )
130 *
131 CALL DGEMV( 'No transpose', M, L, ONE, C( 1, N-L+1 ), LDC,
132 $ V, INCV, ONE, WORK, 1 )
133 *
134 * C( 1:m, 1 ) = C( 1:m, 1 ) - tau * w( 1:m )
135 *
136 CALL DAXPY( M, -TAU, WORK, 1, C, 1 )
137 *
138 * C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
139 * tau * w( 1:m ) * v( 1:l )'
140 *
141 CALL DGER( M, L, -TAU, WORK, 1, V, INCV, C( 1, N-L+1 ),
142 $ LDC )
143 *
144 END IF
145 *
146 END IF
147 *
148 RETURN
149 *
150 * End of DLARZ
151 *
152 END