My Project  debian-1:4.1.2-p1+ds-2
lq.h
Go to the documentation of this file.
1 /*************************************************************************
2 Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3 
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions are
6 met:
7 
8 - Redistributions of source code must retain the above copyright
9  notice, this list of conditions and the following disclaimer.
10 
11 - Redistributions in binary form must reproduce the above copyright
12  notice, this list of conditions and the following disclaimer listed
13  in this license in the documentation and/or other materials
14  provided with the distribution.
15 
16 - Neither the name of the copyright holders nor the names of its
17  contributors may be used to endorse or promote products derived from
18  this software without specific prior written permission.
19 
20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 *************************************************************************/
32 
33 #ifndef _lq_h
34 #define _lq_h
35 
36 #include "ap.h"
37 #include "amp.h"
38 #include "reflections.h"
39 namespace lq
40 {
41  template<unsigned int Precision>
43  int m,
44  int n,
46  template<unsigned int Precision>
48  int m,
49  int n,
51  int qrows,
53  template<unsigned int Precision>
55  int m,
56  int n,
58  template<unsigned int Precision>
60  int m,
61  int n,
63  template<unsigned int Precision>
65  int m,
66  int n,
68  int qrows,
70  template<unsigned int Precision>
72  int m,
73  int n,
76 
77 
78  /*************************************************************************
79  LQ decomposition of a rectangular matrix of size MxN
80 
81  Input parameters:
82  A - matrix A whose indexes range within [0..M-1, 0..N-1].
83  M - number of rows in matrix A.
84  N - number of columns in matrix A.
85 
86  Output parameters:
87  A - matrices L and Q in compact form (see below)
88  Tau - array of scalar factors which are used to form
89  matrix Q. Array whose index ranges within [0..Min(M,N)-1].
90 
91  Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size
92  MxM, L - lower triangular (or lower trapezoid) matrix of size M x N.
93 
94  The elements of matrix L are located on and below the main diagonal of
95  matrix A. The elements which are located in Tau array and above the main
96  diagonal of matrix A are used to form matrix Q as follows:
97 
98  Matrix Q is represented as a product of elementary reflections
99 
100  Q = H(k-1)*H(k-2)*...*H(1)*H(0),
101 
102  where k = min(m,n), and each H(i) is of the form
103 
104  H(i) = 1 - tau * v * (v^T)
105 
106  where tau is a scalar stored in Tau[I]; v - real vector, so that v(0:i-1)=0,
107  v(i) = 1, v(i+1:n-1) stored in A(i,i+1:n-1).
108 
109  -- ALGLIB --
110  Copyright 2005-2007 by Bochkanov Sergey
111  *************************************************************************/
112  template<unsigned int Precision>
114  int m,
115  int n,
117  {
120  int i;
121  int k;
122  int minmn;
123  int maxmn;
125 
126 
127  minmn = ap::minint(m, n);
128  maxmn = ap::maxint(m, n);
129  work.setbounds(0, m);
130  t.setbounds(0, n);
131  tau.setbounds(0, minmn-1);
132  k = ap::minint(m, n);
133  for(i=0; i<=k-1; i++)
134  {
135 
136  //
137  // Generate elementary reflector H(i) to annihilate A(i,i+1:n-1)
138  //
139  ap::vmove(t.getvector(1, n-i), a.getrow(i, i, n-1));
140  reflections::generatereflection<Precision>(t, n-i, tmp);
141  tau(i) = tmp;
142  ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
143  t(1) = 1;
144  if( i<n )
145  {
146 
147  //
148  // Apply H(i) to A(i+1:m,i:n) from the right
149  //
150  reflections::applyreflectionfromtheright<Precision>(a, tau(i), t, i+1, m-1, i, n-1, work);
151  }
152  }
153  }
154 
155 
156  /*************************************************************************
157  Partial unpacking of matrix Q from the LQ decomposition of a matrix A
158 
159  Input parameters:
160  A - matrices L and Q in compact form.
161  Output of RMatrixLQ subroutine.
162  M - number of rows in given matrix A. M>=0.
163  N - number of columns in given matrix A. N>=0.
164  Tau - scalar factors which are used to form Q.
165  Output of the RMatrixLQ subroutine.
166  QRows - required number of rows in matrix Q. N>=QRows>=0.
167 
168  Output parameters:
169  Q - first QRows rows of matrix Q. Array whose indexes range
170  within [0..QRows-1, 0..N-1]. If QRows=0, the array remains
171  unchanged.
172 
173  -- ALGLIB --
174  Copyright 2005 by Bochkanov Sergey
175  *************************************************************************/
176  template<unsigned int Precision>
178  int m,
179  int n,
181  int qrows,
183  {
184  int i;
185  int j;
186  int k;
187  int minmn;
190 
191 
193  if( m<=0 || n<=0 || qrows<=0 )
194  {
195  return;
196  }
197 
198  //
199  // init
200  //
201  minmn = ap::minint(m, n);
202  k = ap::minint(minmn, qrows);
203  q.setbounds(0, qrows-1, 0, n-1);
204  v.setbounds(0, n);
205  work.setbounds(0, qrows);
206  for(i=0; i<=qrows-1; i++)
207  {
208  for(j=0; j<=n-1; j++)
209  {
210  if( i==j )
211  {
212  q(i,j) = 1;
213  }
214  else
215  {
216  q(i,j) = 0;
217  }
218  }
219  }
220 
221  //
222  // unpack Q
223  //
224  for(i=k-1; i>=0; i--)
225  {
226 
227  //
228  // Apply H(i)
229  //
230  ap::vmove(v.getvector(1, n-i), a.getrow(i, i, n-1));
231  v(1) = 1;
232  reflections::applyreflectionfromtheright<Precision>(q, tau(i), v, 0, qrows-1, i, n-1, work);
233  }
234  }
235 
236 
237  /*************************************************************************
238  Unpacking of matrix L from the LQ decomposition of a matrix A
239 
240  Input parameters:
241  A - matrices Q and L in compact form.
242  Output of RMatrixLQ subroutine.
243  M - number of rows in given matrix A. M>=0.
244  N - number of columns in given matrix A. N>=0.
245 
246  Output parameters:
247  L - matrix L, array[0..M-1, 0..N-1].
248 
249  -- ALGLIB --
250  Copyright 2005 by Bochkanov Sergey
251  *************************************************************************/
252  template<unsigned int Precision>
254  int m,
255  int n,
257  {
258  int i;
259  int k;
260 
261 
262  if( m<=0 || n<=0 )
263  {
264  return;
265  }
266  l.setbounds(0, m-1, 0, n-1);
267  for(i=0; i<=n-1; i++)
268  {
269  l(0,i) = 0;
270  }
271  for(i=1; i<=m-1; i++)
272  {
273  ap::vmove(l.getrow(i, 0, n-1), l.getrow(0, 0, n-1));
274  }
275  for(i=0; i<=m-1; i++)
276  {
277  k = ap::minint(i, n-1);
278  ap::vmove(l.getrow(i, 0, k), a.getrow(i, 0, k));
279  }
280  }
281 
282 
283  /*************************************************************************
284  Obsolete 1-based subroutine
285  See RMatrixLQ for 0-based replacement.
286  *************************************************************************/
287  template<unsigned int Precision>
289  int m,
290  int n,
292  {
295  int i;
296  int k;
297  int nmip1;
298  int minmn;
299  int maxmn;
301 
302 
303  minmn = ap::minint(m, n);
304  maxmn = ap::maxint(m, n);
305  work.setbounds(1, m);
306  t.setbounds(1, n);
307  tau.setbounds(1, minmn);
308 
309  //
310  // Test the input arguments
311  //
312  k = ap::minint(m, n);
313  for(i=1; i<=k; i++)
314  {
315 
316  //
317  // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
318  //
319  nmip1 = n-i+1;
320  ap::vmove(t.getvector(1, nmip1), a.getrow(i, i, n));
321  reflections::generatereflection<Precision>(t, nmip1, tmp);
322  tau(i) = tmp;
323  ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
324  t(1) = 1;
325  if( i<n )
326  {
327 
328  //
329  // Apply H(i) to A(i+1:m,i:n) from the right
330  //
331  reflections::applyreflectionfromtheright<Precision>(a, tau(i), t, i+1, m, i, n, work);
332  }
333  }
334  }
335 
336 
337  /*************************************************************************
338  Obsolete 1-based subroutine
339  See RMatrixLQUnpackQ for 0-based replacement.
340  *************************************************************************/
341  template<unsigned int Precision>
343  int m,
344  int n,
346  int qrows,
348  {
349  int i;
350  int j;
351  int k;
352  int minmn;
355  int vm;
356 
357 
359  if( m==0 || n==0 || qrows==0 )
360  {
361  return;
362  }
363 
364  //
365  // init
366  //
367  minmn = ap::minint(m, n);
368  k = ap::minint(minmn, qrows);
369  q.setbounds(1, qrows, 1, n);
370  v.setbounds(1, n);
371  work.setbounds(1, qrows);
372  for(i=1; i<=qrows; i++)
373  {
374  for(j=1; j<=n; j++)
375  {
376  if( i==j )
377  {
378  q(i,j) = 1;
379  }
380  else
381  {
382  q(i,j) = 0;
383  }
384  }
385  }
386 
387  //
388  // unpack Q
389  //
390  for(i=k; i>=1; i--)
391  {
392 
393  //
394  // Apply H(i)
395  //
396  vm = n-i+1;
397  ap::vmove(v.getvector(1, vm), a.getrow(i, i, n));
398  v(1) = 1;
399  reflections::applyreflectionfromtheright<Precision>(q, tau(i), v, 1, qrows, i, n, work);
400  }
401  }
402 
403 
404  /*************************************************************************
405  Obsolete 1-based subroutine
406  *************************************************************************/
407  template<unsigned int Precision>
409  int m,
410  int n,
413  {
414  int i;
415  int j;
417 
418 
419  if( n<=0 )
420  {
421  return;
422  }
423  q.setbounds(1, n, 1, n);
424  l.setbounds(1, m, 1, n);
425 
426  //
427  // LQDecomposition
428  //
429  lqdecomposition<Precision>(a, m, n, tau);
430 
431  //
432  // L
433  //
434  for(i=1; i<=m; i++)
435  {
436  for(j=1; j<=n; j++)
437  {
438  if( j>i )
439  {
440  l(i,j) = 0;
441  }
442  else
443  {
444  l(i,j) = a(i,j);
445  }
446  }
447  }
448 
449  //
450  // Q
451  //
452  unpackqfromlq<Precision>(a, m, n, tau, n, q);
453  }
454 } // namespace
455 
456 #endif
int l
Definition: cfEzgcd.cc:93
int m
Definition: cfEzgcd.cc:121
int i
Definition: cfEzgcd.cc:125
int k
Definition: cfEzgcd.cc:92
void tau(int **points, int sizePoints, int k)
Definition: amp.h:82
static void make_assertion(bool bClause)
Definition: ap.h:49
raw_vector< T > getvector(int iStart, int iEnd)
Definition: ap.h:776
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
int j
Definition: facHensel.cc:105
int maxint(int m1, int m2)
Definition: ap.cpp:162
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:237
int minint(int m1, int m2)
Definition: ap.cpp:167
Definition: lq.h:40
void rmatrixlq(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
Definition: lq.h:113
void rmatrixlqunpackl(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l)
Definition: lq.h:253
void lqdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
Definition: lq.h:288
void rmatrixlqunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: lq.h:177
void unpackqfromlq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: lq.h:342
void lqdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: lq.h:408