My Project  debian-1:4.1.2-p1+ds-2
rotations.h
Go to the documentation of this file.
1 /*************************************************************************
2 Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
3 
4 Contributors:
5  * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6  pseudocode.
7 
8 See subroutines comments for additional copyrights.
9 
10 Redistribution and use in source and binary forms, with or without
11 modification, are permitted provided that the following conditions are
12 met:
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 
26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36 OF 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"
44 namespace 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;
112  amp::ampf<Precision> ctemp;
113  amp::ampf<Precision> stemp;
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;
250  amp::ampf<Precision> ctemp;
251  amp::ampf<Precision> stemp;
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
rotations
Definition: rotations.h:43
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:386
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:230
j
int j
Definition: facHensel.cc:105
f
FILE * f
Definition: checklibs.c:9
ap::vadd
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:403
ap::template_1d_array
Definition: ap.h:641
ap::template_2d_array
Definition: ap.h:790
amp::ampf
Definition: amp.h:81
g
g
Definition: cfModGcd.cc:4031
ap::vmul
void vmul(raw_vector< T > vdst, T2 alpha)
Definition: ap.h:589
ap.h
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:270
ap::vsub
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:521
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
amp.h
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:133