My Project  debian-1:4.1.1-p2+ds-4build3
bidiagonal.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 _bidiagonal_h
40 #define _bidiagonal_h
41 
42 #include "ap.h"
43 #include "amp.h"
44 #include "reflections.h"
45 namespace bidiagonal
46 {
47  template<unsigned int Precision>
49  int m,
50  int n,
53  template<unsigned int Precision>
55  int m,
56  int n,
58  int qcolumns,
60  template<unsigned int Precision>
62  int m,
63  int n,
66  int zrows,
67  int zcolumns,
68  bool fromtheright,
69  bool dotranspose);
70  template<unsigned int Precision>
72  int m,
73  int n,
75  int ptrows,
77  template<unsigned int Precision>
79  int m,
80  int n,
83  int zrows,
84  int zcolumns,
85  bool fromtheright,
86  bool dotranspose);
87  template<unsigned int Precision>
89  int m,
90  int n,
91  bool& isupper,
94  template<unsigned int Precision>
96  int m,
97  int n,
100  template<unsigned int Precision>
102  int m,
103  int n,
105  int qcolumns,
107  template<unsigned int Precision>
109  int m,
110  int n,
113  int zrows,
114  int zcolumns,
115  bool fromtheright,
116  bool dotranspose);
117  template<unsigned int Precision>
119  int m,
120  int n,
122  int ptrows,
124  template<unsigned int Precision>
126  int m,
127  int n,
130  int zrows,
131  int zcolumns,
132  bool fromtheright,
133  bool dotranspose);
134  template<unsigned int Precision>
136  int m,
137  int n,
138  bool& isupper,
141 
142 
143  /*************************************************************************
144  Reduction of a rectangular matrix to bidiagonal form
145 
146  The algorithm reduces the rectangular matrix A to bidiagonal form by
147  orthogonal transformations P and Q: A = Q*B*P.
148 
149  Input parameters:
150  A - source matrix. array[0..M-1, 0..N-1]
151  M - number of rows in matrix A.
152  N - number of columns in matrix A.
153 
154  Output parameters:
155  A - matrices Q, B, P in compact form (see below).
156  TauQ - scalar factors which are used to form matrix Q.
157  TauP - scalar factors which are used to form matrix P.
158 
159  The main diagonal and one of the secondary diagonals of matrix A are
160  replaced with bidiagonal matrix B. Other elements contain elementary
161  reflections which form MxM matrix Q and NxN matrix P, respectively.
162 
163  If M>=N, B is the upper bidiagonal MxN matrix and is stored in the
164  corresponding elements of matrix A. Matrix Q is represented as a
165  product of elementary reflections Q = H(0)*H(1)*...*H(n-1), where
166  H(i) = 1-tau*v*v'. Here tau is a scalar which is stored in TauQ[i], and
167  vector v has the following structure: v(0:i-1)=0, v(i)=1, v(i+1:m-1) is
168  stored in elements A(i+1:m-1,i). Matrix P is as follows: P =
169  G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. Tau is stored in TauP[i],
170  u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).
171 
172  If M<N, B is the lower bidiagonal MxN matrix and is stored in the
173  corresponding elements of matrix A. Q = H(0)*H(1)*...*H(m-2), where
174  H(i) = 1 - tau*v*v', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
175  is stored in elements A(i+2:m-1,i). P = G(0)*G(1)*...*G(m-1),
176  G(i) = 1-tau*u*u', tau is stored in TauP, u(0:i-1)=0, u(i)=1, u(i+1:n-1)
177  is stored in A(i,i+1:n-1).
178 
179  EXAMPLE:
180 
181  m=6, n=5 (m > n): m=5, n=6 (m < n):
182 
183  ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
184  ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
185  ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
186  ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
187  ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
188  ( v1 v2 v3 v4 v5 )
189 
190  Here vi and ui are vectors which form H(i) and G(i), and d and e -
191  are the diagonal and off-diagonal elements of matrix B.
192  *************************************************************************/
193  template<unsigned int Precision>
195  int m,
196  int n,
199  {
202  int minmn;
203  int maxmn;
204  int i;
205  int j;
207 
208 
209 
210  //
211  // Prepare
212  //
213  if( n<=0 || m<=0 )
214  {
215  return;
216  }
217  minmn = ap::minint(m, n);
218  maxmn = ap::maxint(m, n);
219  work.setbounds(0, maxmn);
220  t.setbounds(0, maxmn);
221  if( m>=n )
222  {
223  tauq.setbounds(0, n-1);
224  taup.setbounds(0, n-1);
225  }
226  else
227  {
228  tauq.setbounds(0, m-1);
229  taup.setbounds(0, m-1);
230  }
231  if( m>=n )
232  {
233 
234  //
235  // Reduce to upper bidiagonal form
236  //
237  for(i=0; i<=n-1; i++)
238  {
239 
240  //
241  // Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
242  //
243  ap::vmove(t.getvector(1, m-i), a.getcolumn(i, i, m-1));
244  reflections::generatereflection<Precision>(t, m-i, ltau);
245  tauq(i) = ltau;
246  ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
247  t(1) = 1;
248 
249  //
250  // Apply H(i) to A(i:m-1,i+1:n-1) from the left
251  //
252  reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i, m-1, i+1, n-1, work);
253  if( i<n-1 )
254  {
255 
256  //
257  // Generate elementary reflector G(i) to annihilate
258  // A(i,i+2:n-1)
259  //
260  ap::vmove(t.getvector(1, n-i-1), a.getrow(i, i+1, n-1));
261  reflections::generatereflection<Precision>(t, n-1-i, ltau);
262  taup(i) = ltau;
263  ap::vmove(a.getrow(i, i+1, n-1), t.getvector(1, n-1-i));
264  t(1) = 1;
265 
266  //
267  // Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
268  //
269  reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
270  }
271  else
272  {
273  taup(i) = 0;
274  }
275  }
276  }
277  else
278  {
279 
280  //
281  // Reduce to lower bidiagonal form
282  //
283  for(i=0; i<=m-1; i++)
284  {
285 
286  //
287  // Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
288  //
289  ap::vmove(t.getvector(1, n-i), a.getrow(i, i, n-1));
290  reflections::generatereflection<Precision>(t, n-i, ltau);
291  taup(i) = ltau;
292  ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
293  t(1) = 1;
294 
295  //
296  // Apply G(i) to A(i+1:m-1,i:n-1) from the right
297  //
298  reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m-1, i, n-1, work);
299  if( i<m-1 )
300  {
301 
302  //
303  // Generate elementary reflector H(i) to annihilate
304  // A(i+2:m-1,i)
305  //
306  ap::vmove(t.getvector(1, m-1-i), a.getcolumn(i, i+1, m-1));
307  reflections::generatereflection<Precision>(t, m-1-i, ltau);
308  tauq(i) = ltau;
309  ap::vmove(a.getcolumn(i, i+1, m-1), t.getvector(1, m-1-i));
310  t(1) = 1;
311 
312  //
313  // Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
314  //
315  reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1, m-1, i+1, n-1, work);
316  }
317  else
318  {
319  tauq(i) = 0;
320  }
321  }
322  }
323  }
324 
325 
326  /*************************************************************************
327  Unpacking matrix Q which reduces a matrix to bidiagonal form.
328 
329  Input parameters:
330  QP - matrices Q and P in compact form.
331  Output of ToBidiagonal subroutine.
332  M - number of rows in matrix A.
333  N - number of columns in matrix A.
334  TAUQ - scalar factors which are used to form Q.
335  Output of ToBidiagonal subroutine.
336  QColumns - required number of columns in matrix Q.
337  M>=QColumns>=0.
338 
339  Output parameters:
340  Q - first QColumns columns of matrix Q.
341  Array[0..M-1, 0..QColumns-1]
342  If QColumns=0, the array is not modified.
343 
344  -- ALGLIB --
345  Copyright 2005 by Bochkanov Sergey
346  *************************************************************************/
347  template<unsigned int Precision>
349  int m,
350  int n,
352  int qcolumns,
354  {
355  int i;
356  int j;
357 
358 
359  ap::ap_error::make_assertion(qcolumns<=m);
360  ap::ap_error::make_assertion(qcolumns>=0);
361  if( m==0 || n==0 || qcolumns==0 )
362  {
363  return;
364  }
365 
366  //
367  // prepare Q
368  //
369  q.setbounds(0, m-1, 0, qcolumns-1);
370  for(i=0; i<=m-1; i++)
371  {
372  for(j=0; j<=qcolumns-1; j++)
373  {
374  if( i==j )
375  {
376  q(i,j) = 1;
377  }
378  else
379  {
380  q(i,j) = 0;
381  }
382  }
383  }
384 
385  //
386  // Calculate
387  //
388  rmatrixbdmultiplybyq<Precision>(qp, m, n, tauq, q, m, qcolumns, false, false);
389  }
390 
391 
392  /*************************************************************************
393  Multiplication by matrix Q which reduces matrix A to bidiagonal form.
394 
395  The algorithm allows pre- or post-multiply by Q or Q'.
396 
397  Input parameters:
398  QP - matrices Q and P in compact form.
399  Output of ToBidiagonal subroutine.
400  M - number of rows in matrix A.
401  N - number of columns in matrix A.
402  TAUQ - scalar factors which are used to form Q.
403  Output of ToBidiagonal subroutine.
404  Z - multiplied matrix.
405  array[0..ZRows-1,0..ZColumns-1]
406  ZRows - number of rows in matrix Z. If FromTheRight=False,
407  ZRows=M, otherwise ZRows can be arbitrary.
408  ZColumns - number of columns in matrix Z. If FromTheRight=True,
409  ZColumns=M, otherwise ZColumns can be arbitrary.
410  FromTheRight - pre- or post-multiply.
411  DoTranspose - multiply by Q or Q'.
412 
413  Output parameters:
414  Z - product of Z and Q.
415  Array[0..ZRows-1,0..ZColumns-1]
416  If ZRows=0 or ZColumns=0, the array is not modified.
417 
418  -- ALGLIB --
419  Copyright 2005 by Bochkanov Sergey
420  *************************************************************************/
421  template<unsigned int Precision>
423  int m,
424  int n,
427  int zrows,
428  int zcolumns,
429  bool fromtheright,
430  bool dotranspose)
431  {
432  int i;
433  int i1;
434  int i2;
435  int istep;
438  int mx;
439 
440 
441  if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
442  {
443  return;
444  }
445  ap::ap_error::make_assertion(fromtheright && zcolumns==m || !fromtheright && zrows==m);
446 
447  //
448  // init
449  //
450  mx = ap::maxint(m, n);
451  mx = ap::maxint(mx, zrows);
452  mx = ap::maxint(mx, zcolumns);
453  v.setbounds(0, mx);
454  work.setbounds(0, mx);
455  if( m>=n )
456  {
457 
458  //
459  // setup
460  //
461  if( fromtheright )
462  {
463  i1 = 0;
464  i2 = n-1;
465  istep = +1;
466  }
467  else
468  {
469  i1 = n-1;
470  i2 = 0;
471  istep = -1;
472  }
473  if( dotranspose )
474  {
475  i = i1;
476  i1 = i2;
477  i2 = i;
478  istep = -istep;
479  }
480 
481  //
482  // Process
483  //
484  i = i1;
485  do
486  {
487  ap::vmove(v.getvector(1, m-i), qp.getcolumn(i, i, m-1));
488  v(1) = 1;
489  if( fromtheright )
490  {
491  reflections::applyreflectionfromtheright<Precision>(z, tauq(i), v, 0, zrows-1, i, m-1, work);
492  }
493  else
494  {
495  reflections::applyreflectionfromtheleft<Precision>(z, tauq(i), v, i, m-1, 0, zcolumns-1, work);
496  }
497  i = i+istep;
498  }
499  while( i!=i2+istep );
500  }
501  else
502  {
503 
504  //
505  // setup
506  //
507  if( fromtheright )
508  {
509  i1 = 0;
510  i2 = m-2;
511  istep = +1;
512  }
513  else
514  {
515  i1 = m-2;
516  i2 = 0;
517  istep = -1;
518  }
519  if( dotranspose )
520  {
521  i = i1;
522  i1 = i2;
523  i2 = i;
524  istep = -istep;
525  }
526 
527  //
528  // Process
529  //
530  if( m-1>0 )
531  {
532  i = i1;
533  do
534  {
535  ap::vmove(v.getvector(1, m-i-1), qp.getcolumn(i, i+1, m-1));
536  v(1) = 1;
537  if( fromtheright )
538  {
539  reflections::applyreflectionfromtheright<Precision>(z, tauq(i), v, 0, zrows-1, i+1, m-1, work);
540  }
541  else
542  {
543  reflections::applyreflectionfromtheleft<Precision>(z, tauq(i), v, i+1, m-1, 0, zcolumns-1, work);
544  }
545  i = i+istep;
546  }
547  while( i!=i2+istep );
548  }
549  }
550  }
551 
552 
553  /*************************************************************************
554  Unpacking matrix P which reduces matrix A to bidiagonal form.
555  The subroutine returns transposed matrix P.
556 
557  Input parameters:
558  QP - matrices Q and P in compact form.
559  Output of ToBidiagonal subroutine.
560  M - number of rows in matrix A.
561  N - number of columns in matrix A.
562  TAUP - scalar factors which are used to form P.
563  Output of ToBidiagonal subroutine.
564  PTRows - required number of rows of matrix P^T. N >= PTRows >= 0.
565 
566  Output parameters:
567  PT - first PTRows columns of matrix P^T
568  Array[0..PTRows-1, 0..N-1]
569  If PTRows=0, the array is not modified.
570 
571  -- ALGLIB --
572  Copyright 2005-2007 by Bochkanov Sergey
573  *************************************************************************/
574  template<unsigned int Precision>
576  int m,
577  int n,
579  int ptrows,
581  {
582  int i;
583  int j;
584 
585 
586  ap::ap_error::make_assertion(ptrows<=n);
587  ap::ap_error::make_assertion(ptrows>=0);
588  if( m==0 || n==0 || ptrows==0 )
589  {
590  return;
591  }
592 
593  //
594  // prepare PT
595  //
596  pt.setbounds(0, ptrows-1, 0, n-1);
597  for(i=0; i<=ptrows-1; i++)
598  {
599  for(j=0; j<=n-1; j++)
600  {
601  if( i==j )
602  {
603  pt(i,j) = 1;
604  }
605  else
606  {
607  pt(i,j) = 0;
608  }
609  }
610  }
611 
612  //
613  // Calculate
614  //
615  rmatrixbdmultiplybyp<Precision>(qp, m, n, taup, pt, ptrows, n, true, true);
616  }
617 
618 
619  /*************************************************************************
620  Multiplication by matrix P which reduces matrix A to bidiagonal form.
621 
622  The algorithm allows pre- or post-multiply by P or P'.
623 
624  Input parameters:
625  QP - matrices Q and P in compact form.
626  Output of RMatrixBD subroutine.
627  M - number of rows in matrix A.
628  N - number of columns in matrix A.
629  TAUP - scalar factors which are used to form P.
630  Output of RMatrixBD subroutine.
631  Z - multiplied matrix.
632  Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
633  ZRows - number of rows in matrix Z. If FromTheRight=False,
634  ZRows=N, otherwise ZRows can be arbitrary.
635  ZColumns - number of columns in matrix Z. If FromTheRight=True,
636  ZColumns=N, otherwise ZColumns can be arbitrary.
637  FromTheRight - pre- or post-multiply.
638  DoTranspose - multiply by P or P'.
639 
640  Output parameters:
641  Z - product of Z and P.
642  Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
643  If ZRows=0 or ZColumns=0, the array is not modified.
644 
645  -- ALGLIB --
646  Copyright 2005-2007 by Bochkanov Sergey
647  *************************************************************************/
648  template<unsigned int Precision>
650  int m,
651  int n,
654  int zrows,
655  int zcolumns,
656  bool fromtheright,
657  bool dotranspose)
658  {
659  int i;
662  int mx;
663  int i1;
664  int i2;
665  int istep;
666 
667 
668  if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
669  {
670  return;
671  }
672  ap::ap_error::make_assertion(fromtheright && zcolumns==n || !fromtheright && zrows==n);
673 
674  //
675  // init
676  //
677  mx = ap::maxint(m, n);
678  mx = ap::maxint(mx, zrows);
679  mx = ap::maxint(mx, zcolumns);
680  v.setbounds(0, mx);
681  work.setbounds(0, mx);
682  v.setbounds(0, mx);
683  work.setbounds(0, mx);
684  if( m>=n )
685  {
686 
687  //
688  // setup
689  //
690  if( fromtheright )
691  {
692  i1 = n-2;
693  i2 = 0;
694  istep = -1;
695  }
696  else
697  {
698  i1 = 0;
699  i2 = n-2;
700  istep = +1;
701  }
702  if( !dotranspose )
703  {
704  i = i1;
705  i1 = i2;
706  i2 = i;
707  istep = -istep;
708  }
709 
710  //
711  // Process
712  //
713  if( n-1>0 )
714  {
715  i = i1;
716  do
717  {
718  ap::vmove(v.getvector(1, n-1-i), qp.getrow(i, i+1, n-1));
719  v(1) = 1;
720  if( fromtheright )
721  {
722  reflections::applyreflectionfromtheright<Precision>(z, taup(i), v, 0, zrows-1, i+1, n-1, work);
723  }
724  else
725  {
726  reflections::applyreflectionfromtheleft<Precision>(z, taup(i), v, i+1, n-1, 0, zcolumns-1, work);
727  }
728  i = i+istep;
729  }
730  while( i!=i2+istep );
731  }
732  }
733  else
734  {
735 
736  //
737  // setup
738  //
739  if( fromtheright )
740  {
741  i1 = m-1;
742  i2 = 0;
743  istep = -1;
744  }
745  else
746  {
747  i1 = 0;
748  i2 = m-1;
749  istep = +1;
750  }
751  if( !dotranspose )
752  {
753  i = i1;
754  i1 = i2;
755  i2 = i;
756  istep = -istep;
757  }
758 
759  //
760  // Process
761  //
762  i = i1;
763  do
764  {
765  ap::vmove(v.getvector(1, n-i), qp.getrow(i, i, n-1));
766  v(1) = 1;
767  if( fromtheright )
768  {
769  reflections::applyreflectionfromtheright<Precision>(z, taup(i), v, 0, zrows-1, i, n-1, work);
770  }
771  else
772  {
773  reflections::applyreflectionfromtheleft<Precision>(z, taup(i), v, i, n-1, 0, zcolumns-1, work);
774  }
775  i = i+istep;
776  }
777  while( i!=i2+istep );
778  }
779  }
780 
781 
782  /*************************************************************************
783  Unpacking of the main and secondary diagonals of bidiagonal decomposition
784  of matrix A.
785 
786  Input parameters:
787  B - output of RMatrixBD subroutine.
788  M - number of rows in matrix B.
789  N - number of columns in matrix B.
790 
791  Output parameters:
792  IsUpper - True, if the matrix is upper bidiagonal.
793  otherwise IsUpper is False.
794  D - the main diagonal.
795  Array whose index ranges within [0..Min(M,N)-1].
796  E - the secondary diagonal (upper or lower, depending on
797  the value of IsUpper).
798  Array index ranges within [0..Min(M,N)-1], the last
799  element is not used.
800 
801  -- ALGLIB --
802  Copyright 2005-2007 by Bochkanov Sergey
803  *************************************************************************/
804  template<unsigned int Precision>
806  int m,
807  int n,
808  bool& isupper,
811  {
812  int i;
813 
814 
815  isupper = m>=n;
816  if( m<=0 || n<=0 )
817  {
818  return;
819  }
820  if( isupper )
821  {
822  d.setbounds(0, n-1);
823  e.setbounds(0, n-1);
824  for(i=0; i<=n-2; i++)
825  {
826  d(i) = b(i,i);
827  e(i) = b(i,i+1);
828  }
829  d(n-1) = b(n-1,n-1);
830  }
831  else
832  {
833  d.setbounds(0, m-1);
834  e.setbounds(0, m-1);
835  for(i=0; i<=m-2; i++)
836  {
837  d(i) = b(i,i);
838  e(i) = b(i+1,i);
839  }
840  d(m-1) = b(m-1,m-1);
841  }
842  }
843 
844 
845  /*************************************************************************
846  Obsolete 1-based subroutine.
847  See RMatrixBD for 0-based replacement.
848  *************************************************************************/
849  template<unsigned int Precision>
851  int m,
852  int n,
855  {
858  int minmn;
859  int maxmn;
860  int i;
862  int mmip1;
863  int nmi;
864  int ip1;
865  int nmip1;
866  int mmi;
867 
868 
869  minmn = ap::minint(m, n);
870  maxmn = ap::maxint(m, n);
871  work.setbounds(1, maxmn);
872  t.setbounds(1, maxmn);
873  taup.setbounds(1, minmn);
874  tauq.setbounds(1, minmn);
875  if( m>=n )
876  {
877 
878  //
879  // Reduce to upper bidiagonal form
880  //
881  for(i=1; i<=n; i++)
882  {
883 
884  //
885  // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
886  //
887  mmip1 = m-i+1;
888  ap::vmove(t.getvector(1, mmip1), a.getcolumn(i, i, m));
889  reflections::generatereflection<Precision>(t, mmip1, ltau);
890  tauq(i) = ltau;
891  ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
892  t(1) = 1;
893 
894  //
895  // Apply H(i) to A(i:m,i+1:n) from the left
896  //
897  reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i, m, i+1, n, work);
898  if( i<n )
899  {
900 
901  //
902  // Generate elementary reflector G(i) to annihilate
903  // A(i,i+2:n)
904  //
905  nmi = n-i;
906  ip1 = i+1;
907  ap::vmove(t.getvector(1, nmi), a.getrow(i, ip1, n));
908  reflections::generatereflection<Precision>(t, nmi, ltau);
909  taup(i) = ltau;
910  ap::vmove(a.getrow(i, ip1, n), t.getvector(1, nmi));
911  t(1) = 1;
912 
913  //
914  // Apply G(i) to A(i+1:m,i+1:n) from the right
915  //
916  reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m, i+1, n, work);
917  }
918  else
919  {
920  taup(i) = 0;
921  }
922  }
923  }
924  else
925  {
926 
927  //
928  // Reduce to lower bidiagonal form
929  //
930  for(i=1; i<=m; i++)
931  {
932 
933  //
934  // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
935  //
936  nmip1 = n-i+1;
937  ap::vmove(t.getvector(1, nmip1), a.getrow(i, i, n));
938  reflections::generatereflection<Precision>(t, nmip1, ltau);
939  taup(i) = ltau;
940  ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
941  t(1) = 1;
942 
943  //
944  // Apply G(i) to A(i+1:m,i:n) from the right
945  //
946  reflections::applyreflectionfromtheright<Precision>(a, ltau, t, i+1, m, i, n, work);
947  if( i<m )
948  {
949 
950  //
951  // Generate elementary reflector H(i) to annihilate
952  // A(i+2:m,i)
953  //
954  mmi = m-i;
955  ip1 = i+1;
956  ap::vmove(t.getvector(1, mmi), a.getcolumn(i, ip1, m));
957  reflections::generatereflection<Precision>(t, mmi, ltau);
958  tauq(i) = ltau;
959  ap::vmove(a.getcolumn(i, ip1, m), t.getvector(1, mmi));
960  t(1) = 1;
961 
962  //
963  // Apply H(i) to A(i+1:m,i+1:n) from the left
964  //
965  reflections::applyreflectionfromtheleft<Precision>(a, ltau, t, i+1, m, i+1, n, work);
966  }
967  else
968  {
969  tauq(i) = 0;
970  }
971  }
972  }
973  }
974 
975 
976  /*************************************************************************
977  Obsolete 1-based subroutine.
978  See RMatrixBDUnpackQ for 0-based replacement.
979  *************************************************************************/
980  template<unsigned int Precision>
982  int m,
983  int n,
985  int qcolumns,
987  {
988  int i;
989  int j;
990  int ip1;
993  int vm;
994 
995 
996  ap::ap_error::make_assertion(qcolumns<=m);
997  if( m==0 || n==0 || qcolumns==0 )
998  {
999  return;
1000  }
1001 
1002  //
1003  // init
1004  //
1005  q.setbounds(1, m, 1, qcolumns);
1006  v.setbounds(1, m);
1007  work.setbounds(1, qcolumns);
1008 
1009  //
1010  // prepare Q
1011  //
1012  for(i=1; i<=m; i++)
1013  {
1014  for(j=1; j<=qcolumns; j++)
1015  {
1016  if( i==j )
1017  {
1018  q(i,j) = 1;
1019  }
1020  else
1021  {
1022  q(i,j) = 0;
1023  }
1024  }
1025  }
1026  if( m>=n )
1027  {
1028  for(i=ap::minint(n, qcolumns); i>=1; i--)
1029  {
1030  vm = m-i+1;
1031  ap::vmove(v.getvector(1, vm), qp.getcolumn(i, i, m));
1032  v(1) = 1;
1033  reflections::applyreflectionfromtheleft<Precision>(q, tauq(i), v, i, m, 1, qcolumns, work);
1034  }
1035  }
1036  else
1037  {
1038  for(i=ap::minint(m-1, qcolumns-1); i>=1; i--)
1039  {
1040  vm = m-i;
1041  ip1 = i+1;
1042  ap::vmove(v.getvector(1, vm), qp.getcolumn(i, ip1, m));
1043  v(1) = 1;
1044  reflections::applyreflectionfromtheleft<Precision>(q, tauq(i), v, i+1, m, 1, qcolumns, work);
1045  }
1046  }
1047  }
1048 
1049 
1050  /*************************************************************************
1051  Obsolete 1-based subroutine.
1052  See RMatrixBDMultiplyByQ for 0-based replacement.
1053  *************************************************************************/
1054  template<unsigned int Precision>
1056  int m,
1057  int n,
1060  int zrows,
1061  int zcolumns,
1062  bool fromtheright,
1063  bool dotranspose)
1064  {
1065  int i;
1066  int ip1;
1067  int i1;
1068  int i2;
1069  int istep;
1072  int vm;
1073  int mx;
1074 
1075 
1076  if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
1077  {
1078  return;
1079  }
1080  ap::ap_error::make_assertion(fromtheright && zcolumns==m || !fromtheright && zrows==m);
1081 
1082  //
1083  // init
1084  //
1085  mx = ap::maxint(m, n);
1086  mx = ap::maxint(mx, zrows);
1087  mx = ap::maxint(mx, zcolumns);
1088  v.setbounds(1, mx);
1089  work.setbounds(1, mx);
1090  if( m>=n )
1091  {
1092 
1093  //
1094  // setup
1095  //
1096  if( fromtheright )
1097  {
1098  i1 = 1;
1099  i2 = n;
1100  istep = +1;
1101  }
1102  else
1103  {
1104  i1 = n;
1105  i2 = 1;
1106  istep = -1;
1107  }
1108  if( dotranspose )
1109  {
1110  i = i1;
1111  i1 = i2;
1112  i2 = i;
1113  istep = -istep;
1114  }
1115 
1116  //
1117  // Process
1118  //
1119  i = i1;
1120  do
1121  {
1122  vm = m-i+1;
1123  ap::vmove(v.getvector(1, vm), qp.getcolumn(i, i, m));
1124  v(1) = 1;
1125  if( fromtheright )
1126  {
1127  reflections::applyreflectionfromtheright<Precision>(z, tauq(i), v, 1, zrows, i, m, work);
1128  }
1129  else
1130  {
1131  reflections::applyreflectionfromtheleft<Precision>(z, tauq(i), v, i, m, 1, zcolumns, work);
1132  }
1133  i = i+istep;
1134  }
1135  while( i!=i2+istep );
1136  }
1137  else
1138  {
1139 
1140  //
1141  // setup
1142  //
1143  if( fromtheright )
1144  {
1145  i1 = 1;
1146  i2 = m-1;
1147  istep = +1;
1148  }
1149  else
1150  {
1151  i1 = m-1;
1152  i2 = 1;
1153  istep = -1;
1154  }
1155  if( dotranspose )
1156  {
1157  i = i1;
1158  i1 = i2;
1159  i2 = i;
1160  istep = -istep;
1161  }
1162 
1163  //
1164  // Process
1165  //
1166  if( m-1>0 )
1167  {
1168  i = i1;
1169  do
1170  {
1171  vm = m-i;
1172  ip1 = i+1;
1173  ap::vmove(v.getvector(1, vm), qp.getcolumn(i, ip1, m));
1174  v(1) = 1;
1175  if( fromtheright )
1176  {
1177  reflections::applyreflectionfromtheright<Precision>(z, tauq(i), v, 1, zrows, i+1, m, work);
1178  }
1179  else
1180  {
1181  reflections::applyreflectionfromtheleft<Precision>(z, tauq(i), v, i+1, m, 1, zcolumns, work);
1182  }
1183  i = i+istep;
1184  }
1185  while( i!=i2+istep );
1186  }
1187  }
1188  }
1189 
1190 
1191  /*************************************************************************
1192  Obsolete 1-based subroutine.
1193  See RMatrixBDUnpackPT for 0-based replacement.
1194  *************************************************************************/
1195  template<unsigned int Precision>
1197  int m,
1198  int n,
1200  int ptrows,
1202  {
1203  int i;
1204  int j;
1205  int ip1;
1208  int vm;
1209 
1210 
1211  ap::ap_error::make_assertion(ptrows<=n);
1212  if( m==0 || n==0 || ptrows==0 )
1213  {
1214  return;
1215  }
1216 
1217  //
1218  // init
1219  //
1220  pt.setbounds(1, ptrows, 1, n);
1221  v.setbounds(1, n);
1222  work.setbounds(1, ptrows);
1223 
1224  //
1225  // prepare PT
1226  //
1227  for(i=1; i<=ptrows; i++)
1228  {
1229  for(j=1; j<=n; j++)
1230  {
1231  if( i==j )
1232  {
1233  pt(i,j) = 1;
1234  }
1235  else
1236  {
1237  pt(i,j) = 0;
1238  }
1239  }
1240  }
1241  if( m>=n )
1242  {
1243  for(i=ap::minint(n-1, ptrows-1); i>=1; i--)
1244  {
1245  vm = n-i;
1246  ip1 = i+1;
1247  ap::vmove(v.getvector(1, vm), qp.getrow(i, ip1, n));
1248  v(1) = 1;
1249  reflections::applyreflectionfromtheright<Precision>(pt, taup(i), v, 1, ptrows, i+1, n, work);
1250  }
1251  }
1252  else
1253  {
1254  for(i=ap::minint(m, ptrows); i>=1; i--)
1255  {
1256  vm = n-i+1;
1257  ap::vmove(v.getvector(1, vm), qp.getrow(i, i, n));
1258  v(1) = 1;
1259  reflections::applyreflectionfromtheright<Precision>(pt, taup(i), v, 1, ptrows, i, n, work);
1260  }
1261  }
1262  }
1263 
1264 
1265  /*************************************************************************
1266  Obsolete 1-based subroutine.
1267  See RMatrixBDMultiplyByP for 0-based replacement.
1268  *************************************************************************/
1269  template<unsigned int Precision>
1271  int m,
1272  int n,
1275  int zrows,
1276  int zcolumns,
1277  bool fromtheright,
1278  bool dotranspose)
1279  {
1280  int i;
1281  int ip1;
1284  int vm;
1285  int mx;
1286  int i1;
1287  int i2;
1288  int istep;
1289 
1290 
1291  if( m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
1292  {
1293  return;
1294  }
1295  ap::ap_error::make_assertion(fromtheright && zcolumns==n || !fromtheright && zrows==n);
1296 
1297  //
1298  // init
1299  //
1300  mx = ap::maxint(m, n);
1301  mx = ap::maxint(mx, zrows);
1302  mx = ap::maxint(mx, zcolumns);
1303  v.setbounds(1, mx);
1304  work.setbounds(1, mx);
1305  v.setbounds(1, mx);
1306  work.setbounds(1, mx);
1307  if( m>=n )
1308  {
1309 
1310  //
1311  // setup
1312  //
1313  if( fromtheright )
1314  {
1315  i1 = n-1;
1316  i2 = 1;
1317  istep = -1;
1318  }
1319  else
1320  {
1321  i1 = 1;
1322  i2 = n-1;
1323  istep = +1;
1324  }
1325  if( !dotranspose )
1326  {
1327  i = i1;
1328  i1 = i2;
1329  i2 = i;
1330  istep = -istep;
1331  }
1332 
1333  //
1334  // Process
1335  //
1336  if( n-1>0 )
1337  {
1338  i = i1;
1339  do
1340  {
1341  vm = n-i;
1342  ip1 = i+1;
1343  ap::vmove(v.getvector(1, vm), qp.getrow(i, ip1, n));
1344  v(1) = 1;
1345  if( fromtheright )
1346  {
1347  reflections::applyreflectionfromtheright<Precision>(z, taup(i), v, 1, zrows, i+1, n, work);
1348  }
1349  else
1350  {
1351  reflections::applyreflectionfromtheleft<Precision>(z, taup(i), v, i+1, n, 1, zcolumns, work);
1352  }
1353  i = i+istep;
1354  }
1355  while( i!=i2+istep );
1356  }
1357  }
1358  else
1359  {
1360 
1361  //
1362  // setup
1363  //
1364  if( fromtheright )
1365  {
1366  i1 = m;
1367  i2 = 1;
1368  istep = -1;
1369  }
1370  else
1371  {
1372  i1 = 1;
1373  i2 = m;
1374  istep = +1;
1375  }
1376  if( !dotranspose )
1377  {
1378  i = i1;
1379  i1 = i2;
1380  i2 = i;
1381  istep = -istep;
1382  }
1383 
1384  //
1385  // Process
1386  //
1387  i = i1;
1388  do
1389  {
1390  vm = n-i+1;
1391  ap::vmove(v.getvector(1, vm), qp.getrow(i, i, n));
1392  v(1) = 1;
1393  if( fromtheright )
1394  {
1395  reflections::applyreflectionfromtheright<Precision>(z, taup(i), v, 1, zrows, i, n, work);
1396  }
1397  else
1398  {
1399  reflections::applyreflectionfromtheleft<Precision>(z, taup(i), v, i, n, 1, zcolumns, work);
1400  }
1401  i = i+istep;
1402  }
1403  while( i!=i2+istep );
1404  }
1405  }
1406 
1407 
1408  /*************************************************************************
1409  Obsolete 1-based subroutine.
1410  See RMatrixBDUnpackDiagonals for 0-based replacement.
1411  *************************************************************************/
1412  template<unsigned int Precision>
1414  int m,
1415  int n,
1416  bool& isupper,
1419  {
1420  int i;
1421 
1422 
1423  isupper = m>=n;
1424  if( m==0 || n==0 )
1425  {
1426  return;
1427  }
1428  if( isupper )
1429  {
1430  d.setbounds(1, n);
1431  e.setbounds(1, n);
1432  for(i=1; i<=n-1; i++)
1433  {
1434  d(i) = b(i,i);
1435  e(i) = b(i,i+1);
1436  }
1437  d(n) = b(n,n);
1438  }
1439  else
1440  {
1441  d.setbounds(1, m);
1442  e.setbounds(1, m);
1443  for(i=1; i<=m-1; i++)
1444  {
1445  d(i) = b(i,i);
1446  e(i) = b(i+1,i);
1447  }
1448  d(m) = b(m,m);
1449  }
1450  }
1451 } // namespace
1452 
1453 #endif
bidiagonal::unpackptfrombidiagonal
void unpackptfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
Definition: bidiagonal.h:1232
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:242
j
int j
Definition: facHensel.cc:105
bidiagonal::rmatrixbdunpackq
void rmatrixbdunpackq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: bidiagonal.h:384
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
amp::ampf
Definition: amp.h:83
bidiagonal::tobidiagonal
void tobidiagonal(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
Definition: bidiagonal.h:886
b
CanonicalForm b
Definition: cfModGcd.cc:4044
ap.h
i
int i
Definition: cfEzgcd.cc:125
bidiagonal::multiplybyqfrombidiagonal
void multiplybyqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:1091
bidiagonal::multiplybypfrombidiagonal
void multiplybypfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:1306
bidiagonal::rmatrixbdunpackdiagonals
void rmatrixbdunpackdiagonals(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
Definition: bidiagonal.h:841
bidiagonal::rmatrixbdmultiplybyp
void rmatrixbdmultiplybyp(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:685
bidiagonal::unpackqfrombidiagonal
void unpackqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: bidiagonal.h:1017
bidiagonal::rmatrixbdunpackpt
void rmatrixbdunpackpt(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
Definition: bidiagonal.h:611
bidiagonal::rmatrixbd
void rmatrixbd(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
Definition: bidiagonal.h:230
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
bidiagonal::unpackdiagonalsfrombidiagonal
void unpackdiagonalsfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
Definition: bidiagonal.h:1449
bidiagonal
Definition: bidiagonal.h:46
ap::maxint
int maxint(int m1, int m2)
Definition: ap.cpp:162
m
int m
Definition: cfEzgcd.cc:121
ap::template_1d_array::setbounds
void setbounds(int iLow, int iHigh)
Definition: ap.h:742
reflections.h
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
bidiagonal::rmatrixbdmultiplybyq
void rmatrixbdmultiplybyq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:458
amp.h