Mercurial > octave
comparison libcruft/lapack/clahr2.f @ 7789:82be108cc558
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 CLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) | |
2 * | |
3 * -- LAPACK auxiliary routine (version 3.1) -- | |
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. | |
5 * November 2006 | |
6 * | |
7 * .. Scalar Arguments .. | |
8 INTEGER K, LDA, LDT, LDY, N, NB | |
9 * .. | |
10 * .. Array Arguments .. | |
11 COMPLEX A( LDA, * ), T( LDT, NB ), TAU( NB ), | |
12 $ Y( LDY, NB ) | |
13 * .. | |
14 * | |
15 * Purpose | |
16 * ======= | |
17 * | |
18 * CLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) | |
19 * matrix A so that elements below the k-th subdiagonal are zero. The | |
20 * reduction is performed by an unitary similarity transformation | |
21 * Q' * A * Q. The routine returns the matrices V and T which determine | |
22 * Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. | |
23 * | |
24 * This is an auxiliary routine called by CGEHRD. | |
25 * | |
26 * Arguments | |
27 * ========= | |
28 * | |
29 * N (input) INTEGER | |
30 * The order of the matrix A. | |
31 * | |
32 * K (input) INTEGER | |
33 * The offset for the reduction. Elements below the k-th | |
34 * subdiagonal in the first NB columns are reduced to zero. | |
35 * K < N. | |
36 * | |
37 * NB (input) INTEGER | |
38 * The number of columns to be reduced. | |
39 * | |
40 * A (input/output) COMPLEX array, dimension (LDA,N-K+1) | |
41 * On entry, the n-by-(n-k+1) general matrix A. | |
42 * On exit, the elements on and above the k-th subdiagonal in | |
43 * the first NB columns are overwritten with the corresponding | |
44 * elements of the reduced matrix; the elements below the k-th | |
45 * subdiagonal, with the array TAU, represent the matrix Q as a | |
46 * product of elementary reflectors. The other columns of A are | |
47 * unchanged. See Further Details. | |
48 * | |
49 * LDA (input) INTEGER | |
50 * The leading dimension of the array A. LDA >= max(1,N). | |
51 * | |
52 * TAU (output) COMPLEX array, dimension (NB) | |
53 * The scalar factors of the elementary reflectors. See Further | |
54 * Details. | |
55 * | |
56 * T (output) COMPLEX array, dimension (LDT,NB) | |
57 * The upper triangular matrix T. | |
58 * | |
59 * LDT (input) INTEGER | |
60 * The leading dimension of the array T. LDT >= NB. | |
61 * | |
62 * Y (output) COMPLEX array, dimension (LDY,NB) | |
63 * The n-by-nb matrix Y. | |
64 * | |
65 * LDY (input) INTEGER | |
66 * The leading dimension of the array Y. LDY >= N. | |
67 * | |
68 * Further Details | |
69 * =============== | |
70 * | |
71 * The matrix Q is represented as a product of nb elementary reflectors | |
72 * | |
73 * Q = H(1) H(2) . . . H(nb). | |
74 * | |
75 * Each H(i) has the form | |
76 * | |
77 * H(i) = I - tau * v * v' | |
78 * | |
79 * where tau is a complex scalar, and v is a complex vector with | |
80 * v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in | |
81 * A(i+k+1:n,i), and tau in TAU(i). | |
82 * | |
83 * The elements of the vectors v together form the (n-k+1)-by-nb matrix | |
84 * V which is needed, with T and Y, to apply the transformation to the | |
85 * unreduced part of the matrix, using an update of the form: | |
86 * A := (I - V*T*V') * (A - Y*V'). | |
87 * | |
88 * The contents of A on exit are illustrated by the following example | |
89 * with n = 7, k = 3 and nb = 2: | |
90 * | |
91 * ( a a a a a ) | |
92 * ( a a a a a ) | |
93 * ( a a a a a ) | |
94 * ( h h a a a ) | |
95 * ( v1 h a a a ) | |
96 * ( v1 v2 a a a ) | |
97 * ( v1 v2 a a a ) | |
98 * | |
99 * where a denotes an element of the original matrix A, h denotes a | |
100 * modified element of the upper Hessenberg matrix H, and vi denotes an | |
101 * element of the vector defining H(i). | |
102 * | |
103 * This file is a slight modification of LAPACK-3.0's CLAHRD | |
104 * incorporating improvements proposed by Quintana-Orti and Van de | |
105 * Gejin. Note that the entries of A(1:K,2:NB) differ from those | |
106 * returned by the original LAPACK routine. This function is | |
107 * not backward compatible with LAPACK3.0. | |
108 * | |
109 * ===================================================================== | |
110 * | |
111 * .. Parameters .. | |
112 COMPLEX ZERO, ONE | |
113 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), | |
114 $ ONE = ( 1.0E+0, 0.0E+0 ) ) | |
115 * .. | |
116 * .. Local Scalars .. | |
117 INTEGER I | |
118 COMPLEX EI | |
119 * .. | |
120 * .. External Subroutines .. | |
121 EXTERNAL CAXPY, CCOPY, CGEMM, CGEMV, CLACPY, | |
122 $ CLARFG, CSCAL, CTRMM, CTRMV, CLACGV | |
123 * .. | |
124 * .. Intrinsic Functions .. | |
125 INTRINSIC MIN | |
126 * .. | |
127 * .. Executable Statements .. | |
128 * | |
129 * Quick return if possible | |
130 * | |
131 IF( N.LE.1 ) | |
132 $ RETURN | |
133 * | |
134 DO 10 I = 1, NB | |
135 IF( I.GT.1 ) THEN | |
136 * | |
137 * Update A(K+1:N,I) | |
138 * | |
139 * Update I-th column of A - Y * V' | |
140 * | |
141 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) | |
142 CALL CGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, | |
143 $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) | |
144 CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) | |
145 * | |
146 * Apply I - V * T' * V' to this column (call it b) from the | |
147 * left, using the last column of T as workspace | |
148 * | |
149 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) | |
150 * ( V2 ) ( b2 ) | |
151 * | |
152 * where V1 is unit lower triangular | |
153 * | |
154 * w := V1' * b1 | |
155 * | |
156 CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) | |
157 CALL CTRMV( 'Lower', 'Conjugate transpose', 'UNIT', | |
158 $ I-1, A( K+1, 1 ), | |
159 $ LDA, T( 1, NB ), 1 ) | |
160 * | |
161 * w := w + V2'*b2 | |
162 * | |
163 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, | |
164 $ ONE, A( K+I, 1 ), | |
165 $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) | |
166 * | |
167 * w := T'*w | |
168 * | |
169 CALL CTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT', | |
170 $ I-1, T, LDT, | |
171 $ T( 1, NB ), 1 ) | |
172 * | |
173 * b2 := b2 - V2*w | |
174 * | |
175 CALL CGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, | |
176 $ A( K+I, 1 ), | |
177 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) | |
178 * | |
179 * b1 := b1 - V1*w | |
180 * | |
181 CALL CTRMV( 'Lower', 'NO TRANSPOSE', | |
182 $ 'UNIT', I-1, | |
183 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) | |
184 CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) | |
185 * | |
186 A( K+I-1, I-1 ) = EI | |
187 END IF | |
188 * | |
189 * Generate the elementary reflector H(I) to annihilate | |
190 * A(K+I+1:N,I) | |
191 * | |
192 CALL CLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, | |
193 $ TAU( I ) ) | |
194 EI = A( K+I, I ) | |
195 A( K+I, I ) = ONE | |
196 * | |
197 * Compute Y(K+1:N,I) | |
198 * | |
199 CALL CGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, | |
200 $ ONE, A( K+1, I+1 ), | |
201 $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) | |
202 CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, | |
203 $ ONE, A( K+I, 1 ), LDA, | |
204 $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) | |
205 CALL CGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, | |
206 $ Y( K+1, 1 ), LDY, | |
207 $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) | |
208 CALL CSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) | |
209 * | |
210 * Compute T(1:I,I) | |
211 * | |
212 CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) | |
213 CALL CTRMV( 'Upper', 'No Transpose', 'NON-UNIT', | |
214 $ I-1, T, LDT, | |
215 $ T( 1, I ), 1 ) | |
216 T( I, I ) = TAU( I ) | |
217 * | |
218 10 CONTINUE | |
219 A( K+NB, NB ) = EI | |
220 * | |
221 * Compute Y(1:K,1:NB) | |
222 * | |
223 CALL CLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) | |
224 CALL CTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', | |
225 $ 'UNIT', K, NB, | |
226 $ ONE, A( K+1, 1 ), LDA, Y, LDY ) | |
227 IF( N.GT.K+NB ) | |
228 $ CALL CGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, | |
229 $ NB, N-K-NB, ONE, | |
230 $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, | |
231 $ LDY ) | |
232 CALL CTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', | |
233 $ 'NON-UNIT', K, NB, | |
234 $ ONE, T, LDT, Y, LDY ) | |
235 * | |
236 RETURN | |
237 * | |
238 * End of CLAHR2 | |
239 * | |
240 END |