My Project
Loading...
Searching...
No Matches
rotations.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 _rotations_h
40#define _rotations_h
41
42#include "ap.h"
43#include "amp.h"
44namespace rotations
45{
46 template<unsigned int Precision>
47 void applyrotationsfromtheleft(bool isforward,
48 int m1,
49 int m2,
50 int n1,
51 int n2,
52 const ap::template_1d_array< amp::ampf<Precision> >& c,
53 const ap::template_1d_array< amp::ampf<Precision> >& s,
54 ap::template_2d_array< amp::ampf<Precision> >& a,
55 ap::template_1d_array< amp::ampf<Precision> >& work);
56 template<unsigned int Precision>
57 void applyrotationsfromtheright(bool isforward,
58 int m1,
59 int m2,
60 int n1,
61 int n2,
62 const ap::template_1d_array< amp::ampf<Precision> >& c,
63 const ap::template_1d_array< amp::ampf<Precision> >& s,
64 ap::template_2d_array< amp::ampf<Precision> >& a,
65 ap::template_1d_array< amp::ampf<Precision> >& work);
66 template<unsigned int Precision>
67 void generaterotation(amp::ampf<Precision> f,
68 amp::ampf<Precision> g,
69 amp::ampf<Precision>& cs,
70 amp::ampf<Precision>& sn,
71 amp::ampf<Precision>& r);
72
73
74 /*************************************************************************
75 Application of a sequence of elementary rotations to a matrix
76
77 The algorithm pre-multiplies the matrix by a sequence of rotation
78 transformations which is given by arrays C and S. Depending on the value
79 of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true)
80 rows are rotated, or the rows N and N-1, N-2 and N-3 and so on, are rotated.
81
82 Not the whole matrix but only a part of it is transformed (rows from M1 to
83 M2, columns from N1 to N2). Only the elements of this submatrix are changed.
84
85 Input parameters:
86 IsForward - the sequence of the rotation application.
87 M1,M2 - the range of rows to be transformed.
88 N1, N2 - the range of columns to be transformed.
89 C,S - transformation coefficients.
90 Array whose index ranges within [1..M2-M1].
91 A - processed matrix.
92 WORK - working array whose index ranges within [N1..N2].
93
94 Output parameters:
95 A - transformed matrix.
96
97 Utility subroutine.
98 *************************************************************************/
99 template<unsigned int Precision>
100 void applyrotationsfromtheleft(bool isforward,
101 int m1,
102 int m2,
103 int n1,
104 int n2,
105 const ap::template_1d_array< amp::ampf<Precision> >& c,
106 const ap::template_1d_array< amp::ampf<Precision> >& s,
107 ap::template_2d_array< amp::ampf<Precision> >& a,
108 ap::template_1d_array< amp::ampf<Precision> >& work)
109 {
110 int j;
111 int jp1;
112 amp::ampf<Precision> ctemp;
113 amp::ampf<Precision> stemp;
114 amp::ampf<Precision> temp;
115
116
117 if( m1>m2 || n1>n2 )
118 {
119 return;
120 }
121
122 //
123 // Form P * A
124 //
125 if( isforward )
126 {
127 if( n1!=n2 )
128 {
129
130 //
131 // Common case: N1<>N2
132 //
133 for(j=m1; j<=m2-1; j++)
134 {
135 ctemp = c(j-m1+1);
136 stemp = s(j-m1+1);
137 if( ctemp!=1 || stemp!=0 )
138 {
139 jp1 = j+1;
140 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
141 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
142 ap::vmul(a.getrow(j, n1, n2), ctemp);
143 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
144 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
145 }
146 }
147 }
148 else
149 {
150
151 //
152 // Special case: N1=N2
153 //
154 for(j=m1; j<=m2-1; j++)
155 {
156 ctemp = c(j-m1+1);
157 stemp = s(j-m1+1);
158 if( ctemp!=1 || stemp!=0 )
159 {
160 temp = a(j+1,n1);
161 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
162 a(j,n1) = stemp*temp+ctemp*a(j,n1);
163 }
164 }
165 }
166 }
167 else
168 {
169 if( n1!=n2 )
170 {
171
172 //
173 // Common case: N1<>N2
174 //
175 for(j=m2-1; j>=m1; j--)
176 {
177 ctemp = c(j-m1+1);
178 stemp = s(j-m1+1);
179 if( ctemp!=1 || stemp!=0 )
180 {
181 jp1 = j+1;
182 ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
183 ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
184 ap::vmul(a.getrow(j, n1, n2), ctemp);
185 ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
186 ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
187 }
188 }
189 }
190 else
191 {
192
193 //
194 // Special case: N1=N2
195 //
196 for(j=m2-1; j>=m1; j--)
197 {
198 ctemp = c(j-m1+1);
199 stemp = s(j-m1+1);
200 if( ctemp!=1 || stemp!=0 )
201 {
202 temp = a(j+1,n1);
203 a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
204 a(j,n1) = stemp*temp+ctemp*a(j,n1);
205 }
206 }
207 }
208 }
209 }
210
211
212 /*************************************************************************
213 Application of a sequence of elementary rotations to a matrix
214
215 The algorithm post-multiplies the matrix by a sequence of rotation
216 transformations which is given by arrays C and S. Depending on the value
217 of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true)
218 rows are rotated, or the rows N and N-1, N-2 and N-3 and so on are rotated.
219
220 Not the whole matrix but only a part of it is transformed (rows from M1
221 to M2, columns from N1 to N2). Only the elements of this submatrix are changed.
222
223 Input parameters:
224 IsForward - the sequence of the rotation application.
225 M1,M2 - the range of rows to be transformed.
226 N1, N2 - the range of columns to be transformed.
227 C,S - transformation coefficients.
228 Array whose index ranges within [1..N2-N1].
229 A - processed matrix.
230 WORK - working array whose index ranges within [M1..M2].
231
232 Output parameters:
233 A - transformed matrix.
234
235 Utility subroutine.
236 *************************************************************************/
237 template<unsigned int Precision>
238 void applyrotationsfromtheright(bool isforward,
239 int m1,
240 int m2,
241 int n1,
242 int n2,
243 const ap::template_1d_array< amp::ampf<Precision> >& c,
244 const ap::template_1d_array< amp::ampf<Precision> >& s,
245 ap::template_2d_array< amp::ampf<Precision> >& a,
246 ap::template_1d_array< amp::ampf<Precision> >& work)
247 {
248 int j;
249 int jp1;
250 amp::ampf<Precision> ctemp;
251 amp::ampf<Precision> stemp;
252 amp::ampf<Precision> temp;
253
254
255
256 //
257 // Form A * P'
258 //
259 if( isforward )
260 {
261 if( m1!=m2 )
262 {
263
264 //
265 // Common case: M1<>M2
266 //
267 for(j=n1; j<=n2-1; j++)
268 {
269 ctemp = c(j-n1+1);
270 stemp = s(j-n1+1);
271 if( ctemp!=1 || stemp!=0 )
272 {
273 jp1 = j+1;
274 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
275 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
276 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
277 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
278 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
279 }
280 }
281 }
282 else
283 {
284
285 //
286 // Special case: M1=M2
287 //
288 for(j=n1; j<=n2-1; j++)
289 {
290 ctemp = c(j-n1+1);
291 stemp = s(j-n1+1);
292 if( ctemp!=1 || stemp!=0 )
293 {
294 temp = a(m1,j+1);
295 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
296 a(m1,j) = stemp*temp+ctemp*a(m1,j);
297 }
298 }
299 }
300 }
301 else
302 {
303 if( m1!=m2 )
304 {
305
306 //
307 // Common case: M1<>M2
308 //
309 for(j=n2-1; j>=n1; j--)
310 {
311 ctemp = c(j-n1+1);
312 stemp = s(j-n1+1);
313 if( ctemp!=1 || stemp!=0 )
314 {
315 jp1 = j+1;
316 ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
317 ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
318 ap::vmul(a.getcolumn(j, m1, m2), ctemp);
319 ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
320 ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
321 }
322 }
323 }
324 else
325 {
326
327 //
328 // Special case: M1=M2
329 //
330 for(j=n2-1; j>=n1; j--)
331 {
332 ctemp = c(j-n1+1);
333 stemp = s(j-n1+1);
334 if( ctemp!=1 || stemp!=0 )
335 {
336 temp = a(m1,j+1);
337 a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
338 a(m1,j) = stemp*temp+ctemp*a(m1,j);
339 }
340 }
341 }
342 }
343 }
344
345
346 /*************************************************************************
347 The subroutine generates the elementary rotation, so that:
348
349 [ CS SN ] . [ F ] = [ R ]
350 [ -SN CS ] [ G ] [ 0 ]
351
352 CS**2 + SN**2 = 1
353 *************************************************************************/
354 template<unsigned int Precision>
355 void generaterotation(amp::ampf<Precision> f,
356 amp::ampf<Precision> g,
357 amp::ampf<Precision>& cs,
358 amp::ampf<Precision>& sn,
359 amp::ampf<Precision>& r)
360 {
361 amp::ampf<Precision> f1;
362 amp::ampf<Precision> g1;
363
364
365 if( g==0 )
366 {
367 cs = 1;
368 sn = 0;
369 r = f;
370 }
371 else
372 {
373 if( f==0 )
374 {
375 cs = 0;
376 sn = 1;
377 r = g;
378 }
379 else
380 {
381 f1 = f;
382 g1 = g;
383 r = amp::sqrt<Precision>(amp::sqr<Precision>(f1)+amp::sqr<Precision>(g1));
384 cs = f1/r;
385 sn = g1/r;
386 if( amp::abs<Precision>(f)>amp::abs<Precision>(g) && cs<0 )
387 {
388 cs = -cs;
389 sn = -sn;
390 r = -r;
391 }
392 }
393 }
394 }
395} // namespace
396
397#endif
g
g
Definition cfModGcd.cc:4098
f
FILE * f
Definition checklibs.c:9
ap::template_1d_array
Definition ap.h:657
ap::template_2d_array
Definition ap.h:807
s
const CanonicalForm int s
Definition facAbsFact.cc:51
j
int j
Definition facHensel.cc:110
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition amp.h:713
amp::sqrt
const ampf< Precision > sqrt(const ampf< Precision > &x)
Definition amp.h:740
amp::sqr
const ampf< Precision > sqr(const ampf< Precision > &x)
Definition amp.h:693
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:237
ap::vadd
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:413
ap::vmul
void vmul(raw_vector< T > vdst, T2 alpha)
Definition ap.h:603
ap::vsub
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:533
rotations::applyrotationsfromtheright
void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition rotations.h:238
rotations::generaterotation
void generaterotation(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > &cs, amp::ampf< Precision > &sn, amp::ampf< Precision > &r)
Definition rotations.h:355
rotations::applyrotationsfromtheleft
void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition rotations.h:100
Generated on for My Project by   doxygen 1.17.0 for Singular