My Project  debian-1:4.1.1-p2+ds-4build3
qr.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 _qr_h
40 #define _qr_h
41 
42 #include "ap.h"
43 #include "amp.h"
44 #include "reflections.h"
45 namespace qr
46 {
47  template<unsigned int Precision>
49  int m,
50  int n,
52  template<unsigned int Precision>
54  int m,
55  int n,
57  int qcolumns,
59  template<unsigned int Precision>
61  int m,
62  int n,
64  template<unsigned int Precision>
66  int m,
67  int n,
69  template<unsigned int Precision>
71  int m,
72  int n,
74  int qcolumns,
76  template<unsigned int Precision>
78  int m,
79  int n,
82 
83 
84  /*************************************************************************
85  QR decomposition of a rectangular matrix of size MxN
86 
87  Input parameters:
88  A - matrix A whose indexes range within [0..M-1, 0..N-1].
89  M - number of rows in matrix A.
90  N - number of columns in matrix A.
91 
92  Output parameters:
93  A - matrices Q and R in compact form (see below).
94  Tau - array of scalar factors which are used to form
95  matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].
96 
97  Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
98  MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.
99 
100  The elements of matrix R are located on and above the main diagonal of
101  matrix A. The elements which are located in Tau array and below the main
102  diagonal of matrix A are used to form matrix Q as follows:
103 
104  Matrix Q is represented as a product of elementary reflections
105 
106  Q = H(0)*H(2)*...*H(k-1),
107 
108  where k = min(m,n), and each H(i) is in the form
109 
110  H(i) = 1 - tau * v * (v^T)
111 
112  where tau is a scalar stored in Tau[I]; v - real vector,
113  so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).
114 
115  -- LAPACK routine (version 3.0) --
116  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
117  Courant Institute, Argonne National Lab, and Rice University
118  February 29, 1992.
119  Translation from FORTRAN to pseudocode (AlgoPascal)
120  by Sergey Bochkanov, ALGLIB project, 2005-2007.
121  *************************************************************************/
122  template<unsigned int Precision>
124  int m,
125  int n,
127  {
130  int i;
131  int k;
132  int minmn;
134 
135 
136  if( m<=0 || n<=0 )
137  {
138  return;
139  }
140  minmn = ap::minint(m, n);
141  work.setbounds(0, n-1);
142  t.setbounds(1, m);
143  tau.setbounds(0, minmn-1);
144 
145  //
146  // Test the input arguments
147  //
148  k = minmn;
149  for(i=0; i<=k-1; i++)
150  {
151 
152  //
153  // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
154  //
155  ap::vmove(t.getvector(1, m-i), a.getcolumn(i, i, m-1));
156  reflections::generatereflection<Precision>(t, m-i, tmp);
157  tau(i) = tmp;
158  ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
159  t(1) = 1;
160  if( i<n )
161  {
162 
163  //
164  // Apply H(i) to A(i:m-1,i+1:n-1) from the left
165  //
166  reflections::applyreflectionfromtheleft<Precision>(a, tau(i), t, i, m-1, i+1, n-1, work);
167  }
168  }
169  }
170 
171 
172  /*************************************************************************
173  Partial unpacking of matrix Q from the QR decomposition of a matrix A
174 
175  Input parameters:
176  A - matrices Q and R in compact form.
177  Output of RMatrixQR subroutine.
178  M - number of rows in given matrix A. M>=0.
179  N - number of columns in given matrix A. N>=0.
180  Tau - scalar factors which are used to form Q.
181  Output of the RMatrixQR subroutine.
182  QColumns - required number of columns of matrix Q. M>=QColumns>=0.
183 
184  Output parameters:
185  Q - first QColumns columns of matrix Q.
186  Array whose indexes range within [0..M-1, 0..QColumns-1].
187  If QColumns=0, the array remains unchanged.
188 
189  -- ALGLIB --
190  Copyright 2005 by Bochkanov Sergey
191  *************************************************************************/
192  template<unsigned int Precision>
194  int m,
195  int n,
197  int qcolumns,
199  {
200  int i;
201  int j;
202  int k;
203  int minmn;
206 
207 
208  ap::ap_error::make_assertion(qcolumns<=m);
209  if( m<=0 || n<=0 || qcolumns<=0 )
210  {
211  return;
212  }
213 
214  //
215  // init
216  //
217  minmn = ap::minint(m, n);
218  k = ap::minint(minmn, qcolumns);
219  q.setbounds(0, m-1, 0, qcolumns-1);
220  v.setbounds(1, m);
221  work.setbounds(0, qcolumns-1);
222  for(i=0; i<=m-1; i++)
223  {
224  for(j=0; j<=qcolumns-1; j++)
225  {
226  if( i==j )
227  {
228  q(i,j) = 1;
229  }
230  else
231  {
232  q(i,j) = 0;
233  }
234  }
235  }
236 
237  //
238  // unpack Q
239  //
240  for(i=k-1; i>=0; i--)
241  {
242 
243  //
244  // Apply H(i)
245  //
246  ap::vmove(v.getvector(1, m-i), a.getcolumn(i, i, m-1));
247  v(1) = 1;
248  reflections::applyreflectionfromtheleft<Precision>(q, tau(i), v, i, m-1, 0, qcolumns-1, work);
249  }
250  }
251 
252 
253  /*************************************************************************
254  Unpacking of matrix R from the QR decomposition of a matrix A
255 
256  Input parameters:
257  A - matrices Q and R in compact form.
258  Output of RMatrixQR subroutine.
259  M - number of rows in given matrix A. M>=0.
260  N - number of columns in given matrix A. N>=0.
261 
262  Output parameters:
263  R - matrix R, array[0..M-1, 0..N-1].
264 
265  -- ALGLIB --
266  Copyright 2005 by Bochkanov Sergey
267  *************************************************************************/
268  template<unsigned int Precision>
270  int m,
271  int n,
273  {
274  int i;
275  int k;
276 
277 
278  if( m<=0 || n<=0 )
279  {
280  return;
281  }
282  k = ap::minint(m, n);
283  r.setbounds(0, m-1, 0, n-1);
284  for(i=0; i<=n-1; i++)
285  {
286  r(0,i) = 0;
287  }
288  for(i=1; i<=m-1; i++)
289  {
290  ap::vmove(r.getrow(i, 0, n-1), r.getrow(0, 0, n-1));
291  }
292  for(i=0; i<=k-1; i++)
293  {
294  ap::vmove(r.getrow(i, i, n-1), a.getrow(i, i, n-1));
295  }
296  }
297 
298 
299  /*************************************************************************
300  Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
301  *************************************************************************/
302  template<unsigned int Precision>
304  int m,
305  int n,
307  {
310  int i;
311  int k;
312  int mmip1;
313  int minmn;
315 
316 
317  minmn = ap::minint(m, n);
318  work.setbounds(1, n);
319  t.setbounds(1, m);
320  tau.setbounds(1, minmn);
321 
322  //
323  // Test the input arguments
324  //
325  k = ap::minint(m, n);
326  for(i=1; i<=k; i++)
327  {
328 
329  //
330  // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
331  //
332  mmip1 = m-i+1;
333  ap::vmove(t.getvector(1, mmip1), a.getcolumn(i, i, m));
334  reflections::generatereflection<Precision>(t, mmip1, tmp);
335  tau(i) = tmp;
336  ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
337  t(1) = 1;
338  if( i<n )
339  {
340 
341  //
342  // Apply H(i) to A(i:m,i+1:n) from the left
343  //
344  reflections::applyreflectionfromtheleft<Precision>(a, tau(i), t, i, m, i+1, n, work);
345  }
346  }
347  }
348 
349 
350  /*************************************************************************
351  Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement.
352  *************************************************************************/
353  template<unsigned int Precision>
355  int m,
356  int n,
358  int qcolumns,
360  {
361  int i;
362  int j;
363  int k;
364  int minmn;
367  int vm;
368 
369 
370  ap::ap_error::make_assertion(qcolumns<=m);
371  if( m==0 || n==0 || qcolumns==0 )
372  {
373  return;
374  }
375 
376  //
377  // init
378  //
379  minmn = ap::minint(m, n);
380  k = ap::minint(minmn, qcolumns);
381  q.setbounds(1, m, 1, qcolumns);
382  v.setbounds(1, m);
383  work.setbounds(1, qcolumns);
384  for(i=1; i<=m; i++)
385  {
386  for(j=1; j<=qcolumns; j++)
387  {
388  if( i==j )
389  {
390  q(i,j) = 1;
391  }
392  else
393  {
394  q(i,j) = 0;
395  }
396  }
397  }
398 
399  //
400  // unpack Q
401  //
402  for(i=k; i>=1; i--)
403  {
404 
405  //
406  // Apply H(i)
407  //
408  vm = m-i+1;
409  ap::vmove(v.getvector(1, vm), a.getcolumn(i, i, m));
410  v(1) = 1;
411  reflections::applyreflectionfromtheleft<Precision>(q, tau(i), v, i, m, 1, qcolumns, work);
412  }
413  }
414 
415 
416  /*************************************************************************
417  Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
418  *************************************************************************/
419  template<unsigned int Precision>
421  int m,
422  int n,
425  {
426  int i;
427  int k;
431 
432 
433  k = ap::minint(m, n);
434  if( n<=0 )
435  {
436  return;
437  }
438  work.setbounds(1, m);
439  v.setbounds(1, m);
440  q.setbounds(1, m, 1, m);
441  r.setbounds(1, m, 1, n);
442 
443  //
444  // QRDecomposition
445  //
446  qrdecomposition<Precision>(a, m, n, tau);
447 
448  //
449  // R
450  //
451  for(i=1; i<=n; i++)
452  {
453  r(1,i) = 0;
454  }
455  for(i=2; i<=m; i++)
456  {
457  ap::vmove(r.getrow(i, 1, n), r.getrow(1, 1, n));
458  }
459  for(i=1; i<=k; i++)
460  {
461  ap::vmove(r.getrow(i, i, n), a.getrow(i, i, n));
462  }
463 
464  //
465  // Q
466  //
467  unpackqfromqr<Precision>(a, m, n, tau, m, q);
468  }
469 } // namespace
470 
471 #endif
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:242
j
int j
Definition: facHensel.cc:105
k
int k
Definition: cfEzgcd.cc:92
ap::template_1d_array::getvector
raw_vector< T > getvector(int iStart, int iEnd)
Definition: ap.h:783
ap::template_1d_array
Definition: ap.h:662
ap::template_2d_array
Definition: ap.h:812
qr::unpackqfromqr
void unpackqfromqr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: qr.h:390
amp::ampf
Definition: amp.h:83
qr::rmatrixqrunpackr
void rmatrixqrunpackr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &r)
Definition: qr.h:305
ap.h
i
int i
Definition: cfEzgcd.cc:125
ap::minint
int minint(int m1, int m2)
Definition: ap.cpp:167
ap::ap_error::make_assertion
static void make_assertion(bool bClause)
Definition: ap.h:61
qr::qrdecompositionunpacked
void qrdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &q, ap::template_2d_array< amp::ampf< Precision > > &r)
Definition: qr.h:456
tau
void tau(int **points, int sizePoints, int k)
Definition: cfNewtonPolygon.cc:461
m
int m
Definition: cfEzgcd.cc:121
qr::rmatrixqr
void rmatrixqr(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
Definition: qr.h:159
ap::template_1d_array::setbounds
void setbounds(int iLow, int iHigh)
Definition: ap.h:742
qr::qrdecomposition
void qrdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
Definition: qr.h:339
reflections.h
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
qr::rmatrixqrunpackq
void rmatrixqrunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: qr.h:229
qr
Definition: qr.h:46
amp.h