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,
56 template<unsigned int Precision>
57 void applyrotationsfromtheright(bool isforward,
58 int m1,
59 int m2,
60 int n1,
61 int n2,
66 template<unsigned int Precision>
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,
109 {
110 int j;
111 int jp1;
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,
247 {
248 int j;
249 int jp1;
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>
360 {
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
Definition: cfModGcd.cc:4090
FILE * f
Definition: checklibs.c:9
Definition: amp.h:82
const CanonicalForm int s
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:237
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 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
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
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