My Project
Loading...
Searching...
No Matches
reflections.h
Go to the documentation of this file.
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
3
4Contributors:
5 * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6 pseudocode.
7
8See subroutines comments for additional copyrights.
9
10Redistribution and use in source and binary forms, with or without
11modification, are permitted provided that the following conditions are
12met:
13
14- Redistributions of source code must retain the above copyright
15 notice, this list of conditions and the following disclaimer.
16
17- Redistributions in binary form must reproduce the above copyright
18 notice, this list of conditions and the following disclaimer listed
19 in this license in the documentation and/or other materials
20 provided with the distribution.
21
22- Neither the name of the copyright holders nor the names of its
23 contributors may be used to endorse or promote products derived from
24 this software without specific prior written permission.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37*************************************************************************/
38
39#ifndef _reflections_h
40#define _reflections_h
41
42#include "ap.h"
43#include "amp.h"
44namespace reflections
45{
46 template<unsigned int Precision>
48 int n,
50 template<unsigned int Precision>
54 int m1,
55 int m2,
56 int n1,
57 int n2,
59 template<unsigned int Precision>
63 int m1,
64 int m2,
65 int n1,
66 int n2,
68
69
70 /*************************************************************************
71 Generation of an elementary reflection transformation
72
73 The subroutine generates elementary reflection H of order N, so that, for
74 a given X, the following equality holds true:
75
76 ( X(1) ) ( Beta )
77 H * ( .. ) = ( 0 )
78 ( X(n) ) ( 0 )
79
80 where
81 ( V(1) )
82 H = 1 - Tau * ( .. ) * ( V(1), ..., V(n) )
83 ( V(n) )
84
85 where the first component of vector V equals 1.
86
87 Input parameters:
88 X - vector. Array whose index ranges within [1..N].
89 N - reflection order.
90
91 Output parameters:
92 X - components from 2 to N are replaced with vector V.
93 The first component is replaced with parameter Beta.
94 Tau - scalar value Tau. If X is a null vector, Tau equals 0,
95 otherwise 1 <= Tau <= 2.
96
97 This subroutine is the modification of the DLARFG subroutines from
98 the LAPACK library. It has a similar functionality except for the
99 fact that it doesn't handle errors when the intermediate results
100 cause an overflow.
101
102
103 MODIFICATIONS:
104 24.12.2005 sign(Alpha) was replaced with an analogous to the Fortran SIGN code.
105
106 -- LAPACK auxiliary routine (version 3.0) --
107 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
108 Courant Institute, Argonne National Lab, and Rice University
109 September 30, 1994
110 *************************************************************************/
111 template<unsigned int Precision>
113 int n,
115 {
116 int j;
122
123
124
125 //
126 // Executable Statements ..
127 //
128 if( n<=1 )
129 {
130 tau = 0;
131 return;
132 }
133
134 //
135 // XNORM = DNRM2( N-1, X, INCX )
136 //
137 alpha = x(1);
138 mx = 0;
139 for(j=2; j<=n; j++)
140 {
141 mx = amp::maximum<Precision>(amp::abs<Precision>(x(j)), mx);
142 }
143 xnorm = 0;
144 if( mx!=0 )
145 {
146 for(j=2; j<=n; j++)
147 {
148 xnorm = xnorm+amp::sqr<Precision>(x(j)/mx);
149 }
150 xnorm = amp::sqrt<Precision>(xnorm)*mx;
151 }
152 if( xnorm==0 )
153 {
154
155 //
156 // H = I
157 //
158 tau = 0;
159 return;
160 }
161
162 //
163 // general case
164 //
165 mx = amp::maximum<Precision>(amp::abs<Precision>(alpha), amp::abs<Precision>(xnorm));
166 beta = -mx*amp::sqrt<Precision>(amp::sqr<Precision>(alpha/mx)+amp::sqr<Precision>(xnorm/mx));
167 if( alpha<0 )
168 {
169 beta = -beta;
170 }
171 tau = (beta-alpha)/beta;
172 v = 1/(alpha-beta);
173 ap::vmul(x.getvector(2, n), v);
174 x(1) = beta;
175 }
176
177
178 /*************************************************************************
179 Application of an elementary reflection to a rectangular matrix of size MxN
180
181 The algorithm pre-multiplies the matrix by an elementary reflection transformation
182 which is given by column V and scalar Tau (see the description of the
183 GenerateReflection procedure). Not the whole matrix but only a part of it
184 is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements
185 of this submatrix are changed.
186
187 Input parameters:
188 C - matrix to be transformed.
189 Tau - scalar defining the transformation.
190 V - column defining the transformation.
191 Array whose index ranges within [1..M2-M1+1].
192 M1, M2 - range of rows to be transformed.
193 N1, N2 - range of columns to be transformed.
194 WORK - working array whose indexes goes from N1 to N2.
195
196 Output parameters:
197 C - the result of multiplying the input matrix C by the
198 transformation matrix which is given by Tau and V.
199 If N1>N2 or M1>M2, C is not modified.
200
201 -- LAPACK auxiliary routine (version 3.0) --
202 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
203 Courant Institute, Argonne National Lab, and Rice University
204 September 30, 1994
205 *************************************************************************/
206 template<unsigned int Precision>
210 int m1,
211 int m2,
212 int n1,
213 int n2,
215 {
217 int i;
218 int vm;
219
220
221 if( tau==0 || n1>n2 || m1>m2 )
222 {
223 return;
224 }
225
226 //
227 // w := C' * v
228 //
229 vm = m2-m1+1;
230 for(i=n1; i<=n2; i++)
231 {
232 work(i) = 0;
233 }
234 for(i=m1; i<=m2; i++)
235 {
236 t = v(i+1-m1);
237 ap::vadd(work.getvector(n1, n2), c.getrow(i, n1, n2), t);
238 }
239
240 //
241 // C := C - tau * v * w'
242 //
243 for(i=m1; i<=m2; i++)
244 {
245 t = v(i-m1+1)*tau;
246 ap::vsub(c.getrow(i, n1, n2), work.getvector(n1, n2), t);
247 }
248 }
249
250
251 /*************************************************************************
252 Application of an elementary reflection to a rectangular matrix of size MxN
253
254 The algorithm post-multiplies the matrix by an elementary reflection transformation
255 which is given by column V and scalar Tau (see the description of the
256 GenerateReflection procedure). Not the whole matrix but only a part of it
257 is transformed (rows from M1 to M2, columns from N1 to N2). Only the
258 elements of this submatrix are changed.
259
260 Input parameters:
261 C - matrix to be transformed.
262 Tau - scalar defining the transformation.
263 V - column defining the transformation.
264 Array whose index ranges within [1..N2-N1+1].
265 M1, M2 - range of rows to be transformed.
266 N1, N2 - range of columns to be transformed.
267 WORK - working array whose indexes goes from M1 to M2.
268
269 Output parameters:
270 C - the result of multiplying the input matrix C by the
271 transformation matrix which is given by Tau and V.
272 If N1>N2 or M1>M2, C is not modified.
273
274 -- LAPACK auxiliary routine (version 3.0) --
275 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
276 Courant Institute, Argonne National Lab, and Rice University
277 September 30, 1994
278 *************************************************************************/
279 template<unsigned int Precision>
283 int m1,
284 int m2,
285 int n1,
286 int n2,
288 {
290 int i;
291 int vm;
292
293
294 if( tau==0 || n1>n2 || m1>m2 )
295 {
296 return;
297 }
298
299 //
300 // w := C * v
301 //
302 vm = n2-n1+1;
303 for(i=m1; i<=m2; i++)
304 {
305 t = ap::vdotproduct(c.getrow(i, n1, n2), v.getvector(1, vm));
306 work(i) = t;
307 }
308
309 //
310 // C := C - w * v'
311 //
312 for(i=m1; i<=m2; i++)
313 {
314 t = work(i)*tau;
315 ap::vsub(c.getrow(i, n1, n2), v.getvector(1, vm), t);
316 }
317 }
318} // namespace
319
320#endif
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
void tau(int **points, int sizePoints, int k)
Definition: amp.h:82
Variable alpha
Definition: facAbsBiFact.cc:51
Variable beta
Definition: facAbsFact.cc:95
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
T vdotproduct(const_raw_vector< T > v1, const_raw_vector< T > v2)
Definition: ap.h:181
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:413
void vmul(raw_vector< T > vdst, T2 alpha)
Definition: ap.h:603
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:533
void generatereflection(ap::template_1d_array< amp::ampf< Precision > > &x, int n, amp::ampf< Precision > &tau)
Definition: reflections.h:112
void applyreflectionfromtheright(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition: reflections.h:280
void applyreflectionfromtheleft(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition: reflections.h:207