My Project  debian-1:4.1.1-p2+ds-4build4
svd.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 _svd_h
34 #define _svd_h
35 
36 /*MAKEHEADER*/
37 #include "ap.h"
38 
39 /*MAKEHEADER*/
40 #include "amp.h"
41 
42 /*MAKEHEADER*/
43 #include "reflections.h"
44 
45 /*MAKEHEADER*/
46 #include "bidiagonal.h"
47 
48 /*MAKEHEADER*/
49 #include "qr.h"
50 
51 /*MAKEHEADER*/
52 #include "lq.h"
53 
54 /*MAKEHEADER*/
55 #include "blas.h"
56 
57 /*MAKEHEADER*/
58 #include "rotations.h"
59 
60 /*MAKEHEADER*/
61 #include "bdsvd.h"
62 
63 namespace svd
64 {
65  template<unsigned int Precision>
67  int m,
68  int n,
69  int uneeded,
70  int vtneeded,
71  int additionalmemory,
75  template<unsigned int Precision>
77  int m,
78  int n,
79  int uneeded,
80  int vtneeded,
81  int additionalmemory,
85 
86 
87  /*************************************************************************
88  Singular value decomposition of a rectangular matrix.
89 
90  The algorithm calculates the singular value decomposition of a matrix of
91  size MxN: A = U * S * V^T
92 
93  The algorithm finds the singular values and, optionally, matrices U and V^T.
94  The algorithm can find both first min(M,N) columns of matrix U and rows of
95  matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
96  and NxN respectively).
97 
98  Take into account that the subroutine does not return matrix V but V^T.
99 
100  Input parameters:
101  A - matrix to be decomposed.
102  Array whose indexes range within [0..M-1, 0..N-1].
103  M - number of rows in matrix A.
104  N - number of columns in matrix A.
105  UNeeded - 0, 1 or 2. See the description of the parameter U.
106  VTNeeded - 0, 1 or 2. See the description of the parameter VT.
107  AdditionalMemory -
108  If the parameter:
109  * equals 0, the algorithm doesn’t use additional
110  memory (lower requirements, lower performance).
111  * equals 1, the algorithm uses additional
112  memory of size min(M,N)*min(M,N) of real numbers.
113  It often speeds up the algorithm.
114  * equals 2, the algorithm uses additional
115  memory of size M*min(M,N) of real numbers.
116  It allows to get a maximum performance.
117  The recommended value of the parameter is 2.
118 
119  Output parameters:
120  W - contains singular values in descending order.
121  U - if UNeeded=0, U isn't changed, the left singular vectors
122  are not calculated.
123  if Uneeded=1, U contains left singular vectors (first
124  min(M,N) columns of matrix U). Array whose indexes range
125  within [0..M-1, 0..Min(M,N)-1].
126  if UNeeded=2, U contains matrix U wholly. Array whose
127  indexes range within [0..M-1, 0..M-1].
128  VT - if VTNeeded=0, VT isn’t changed, the right singular vectors
129  are not calculated.
130  if VTNeeded=1, VT contains right singular vectors (first
131  min(M,N) rows of matrix V^T). Array whose indexes range
132  within [0..min(M,N)-1, 0..N-1].
133  if VTNeeded=2, VT contains matrix V^T wholly. Array whose
134  indexes range within [0..N-1, 0..N-1].
135 
136  -- ALGLIB --
137  Copyright 2005 by Bochkanov Sergey
138  *************************************************************************/
139  template<unsigned int Precision>
141  int m,
142  int n,
143  int uneeded,
144  int vtneeded,
145  int additionalmemory,
149  {
150  bool result;
157  bool isupper;
158  int minmn;
159  int ncu;
160  int nrvt;
161  int nru;
162  int ncvt;
163  int i;
164  int j;
165  int im1;
167 
168 
169  result = true;
170  if( m==0 || n==0 )
171  {
172  return result;
173  }
174  ap::ap_error::make_assertion(uneeded>=0 && uneeded<=2);
175  ap::ap_error::make_assertion(vtneeded>=0 && vtneeded<=2);
176  ap::ap_error::make_assertion(additionalmemory>=0 && additionalmemory<=2);
177 
178  //
179  // initialize
180  //
181  minmn = ap::minint(m, n);
182  w.setbounds(1, minmn);
183  ncu = 0;
184  nru = 0;
185  if( uneeded==1 )
186  {
187  nru = m;
188  ncu = minmn;
189  u.setbounds(0, nru-1, 0, ncu-1);
190  }
191  if( uneeded==2 )
192  {
193  nru = m;
194  ncu = m;
195  u.setbounds(0, nru-1, 0, ncu-1);
196  }
197  nrvt = 0;
198  ncvt = 0;
199  if( vtneeded==1 )
200  {
201  nrvt = minmn;
202  ncvt = n;
203  vt.setbounds(0, nrvt-1, 0, ncvt-1);
204  }
205  if( vtneeded==2 )
206  {
207  nrvt = n;
208  ncvt = n;
209  vt.setbounds(0, nrvt-1, 0, ncvt-1);
210  }
211 
212  //
213  // M much larger than N
214  // Use bidiagonal reduction with QR-decomposition
215  //
216  if( m>amp::ampf<Precision>("1.6")*n )
217  {
218  if( uneeded==0 )
219  {
220 
221  //
222  // No left singular vectors to be computed
223  //
224  qr::rmatrixqr<Precision>(a, m, n, tau);
225  for(i=0; i<=n-1; i++)
226  {
227  for(j=0; j<=i-1; j++)
228  {
229  a(i,j) = 0;
230  }
231  }
232  bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
233  bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
234  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper, w, e);
235  result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
236  return result;
237  }
238  else
239  {
240 
241  //
242  // Left singular vectors (may be full matrix U) to be computed
243  //
244  qr::rmatrixqr<Precision>(a, m, n, tau);
245  qr::rmatrixqrunpackq<Precision>(a, m, n, tau, ncu, u);
246  for(i=0; i<=n-1; i++)
247  {
248  for(j=0; j<=i-1; j++)
249  {
250  a(i,j) = 0;
251  }
252  }
253  bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
254  bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
255  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper, w, e);
256  if( additionalmemory<1 )
257  {
258 
259  //
260  // No additional memory can be used
261  //
262  bidiagonal::rmatrixbdmultiplybyq<Precision>(a, n, n, tauq, u, m, n, true, false);
263  result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
264  }
265  else
266  {
267 
268  //
269  // Large U. Transforming intermediate matrix T2
270  //
271  work.setbounds(1, ap::maxint(m, n));
272  bidiagonal::rmatrixbdunpackq<Precision>(a, n, n, tauq, n, t2);
273  blas::copymatrix<Precision>(u, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
274  blas::inplacetranspose<Precision>(t2, 0, n-1, 0, n-1, work);
275  result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
276  blas::matrixmatrixmultiply<Precision>(a, 0, m-1, 0, n-1, false, t2, 0, n-1, 0, n-1, true, amp::ampf<Precision>("1.0"), u, 0, m-1, 0, n-1, amp::ampf<Precision>("0.0"), work);
277  }
278  return result;
279  }
280  }
281 
282  //
283  // N much larger than M
284  // Use bidiagonal reduction with LQ-decomposition
285  //
286  if( n>amp::ampf<Precision>("1.6")*m )
287  {
288  if( vtneeded==0 )
289  {
290 
291  //
292  // No right singular vectors to be computed
293  //
294  lq::rmatrixlq<Precision>(a, m, n, tau);
295  for(i=0; i<=m-1; i++)
296  {
297  for(j=i+1; j<=m-1; j++)
298  {
299  a(i,j) = 0;
300  }
301  }
302  bidiagonal::rmatrixbd<Precision>(a, m, m, tauq, taup);
303  bidiagonal::rmatrixbdunpackq<Precision>(a, m, m, tauq, ncu, u);
304  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, m, isupper, w, e);
305  work.setbounds(1, m);
306  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
307  result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
308  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
309  return result;
310  }
311  else
312  {
313 
314  //
315  // Right singular vectors (may be full matrix VT) to be computed
316  //
317  lq::rmatrixlq<Precision>(a, m, n, tau);
318  lq::rmatrixlqunpackq<Precision>(a, m, n, tau, nrvt, vt);
319  for(i=0; i<=m-1; i++)
320  {
321  for(j=i+1; j<=m-1; j++)
322  {
323  a(i,j) = 0;
324  }
325  }
326  bidiagonal::rmatrixbd<Precision>(a, m, m, tauq, taup);
327  bidiagonal::rmatrixbdunpackq<Precision>(a, m, m, tauq, ncu, u);
328  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, m, isupper, w, e);
329  work.setbounds(1, ap::maxint(m, n));
330  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
331  if( additionalmemory<1 )
332  {
333 
334  //
335  // No additional memory available
336  //
337  bidiagonal::rmatrixbdmultiplybyp<Precision>(a, m, m, taup, vt, m, n, false, true);
338  result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, n);
339  }
340  else
341  {
342 
343  //
344  // Large VT. Transforming intermediate matrix T2
345  //
346  bidiagonal::rmatrixbdunpackpt<Precision>(a, m, m, taup, m, t2);
347  result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, t2, m);
348  blas::copymatrix<Precision>(vt, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
349  blas::matrixmatrixmultiply<Precision>(t2, 0, m-1, 0, m-1, false, a, 0, m-1, 0, n-1, false, amp::ampf<Precision>("1.0"), vt, 0, m-1, 0, n-1, amp::ampf<Precision>("0.0"), work);
350  }
351  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
352  return result;
353  }
354  }
355 
356  //
357  // M<=N
358  // We can use inplace transposition of U to get rid of columnwise operations
359  //
360  if( m<=n )
361  {
362  bidiagonal::rmatrixbd<Precision>(a, m, n, tauq, taup);
363  bidiagonal::rmatrixbdunpackq<Precision>(a, m, n, tauq, ncu, u);
364  bidiagonal::rmatrixbdunpackpt<Precision>(a, m, n, taup, nrvt, vt);
365  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, n, isupper, w, e);
366  work.setbounds(1, m);
367  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
368  result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
369  blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
370  return result;
371  }
372 
373  //
374  // Simple bidiagonal reduction
375  //
376  bidiagonal::rmatrixbd<Precision>(a, m, n, tauq, taup);
377  bidiagonal::rmatrixbdunpackq<Precision>(a, m, n, tauq, ncu, u);
378  bidiagonal::rmatrixbdunpackpt<Precision>(a, m, n, taup, nrvt, vt);
379  bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, n, isupper, w, e);
380  if( additionalmemory<2 || uneeded==0 )
381  {
382 
383  //
384  // We cant use additional memory or there is no need in such operations
385  //
386  result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
387  }
388  else
389  {
390 
391  //
392  // We can use additional memory
393  //
394  t2.setbounds(0, minmn-1, 0, m-1);
395  blas::copyandtranspose<Precision>(u, 0, m-1, 0, minmn-1, t2, 0, minmn-1, 0, m-1);
396  result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
397  blas::copyandtranspose<Precision>(t2, 0, minmn-1, 0, m-1, u, 0, m-1, 0, minmn-1);
398  }
399  return result;
400  }
401 
402 
403  /*************************************************************************
404  Obsolete 1-based subroutine.
405  See RMatrixSVD for 0-based replacement.
406  *************************************************************************/
407  template<unsigned int Precision>
409  int m,
410  int n,
411  int uneeded,
412  int vtneeded,
413  int additionalmemory,
417  {
418  bool result;
425  bool isupper;
426  int minmn;
427  int ncu;
428  int nrvt;
429  int nru;
430  int ncvt;
431  int i;
432  int j;
433  int im1;
435 
436 
437  result = true;
438  if( m==0 || n==0 )
439  {
440  return result;
441  }
442  ap::ap_error::make_assertion(uneeded>=0 && uneeded<=2);
443  ap::ap_error::make_assertion(vtneeded>=0 && vtneeded<=2);
444  ap::ap_error::make_assertion(additionalmemory>=0 && additionalmemory<=2);
445 
446  //
447  // initialize
448  //
449  minmn = ap::minint(m, n);
450  w.setbounds(1, minmn);
451  ncu = 0;
452  nru = 0;
453  if( uneeded==1 )
454  {
455  nru = m;
456  ncu = minmn;
457  u.setbounds(1, nru, 1, ncu);
458  }
459  if( uneeded==2 )
460  {
461  nru = m;
462  ncu = m;
463  u.setbounds(1, nru, 1, ncu);
464  }
465  nrvt = 0;
466  ncvt = 0;
467  if( vtneeded==1 )
468  {
469  nrvt = minmn;
470  ncvt = n;
471  vt.setbounds(1, nrvt, 1, ncvt);
472  }
473  if( vtneeded==2 )
474  {
475  nrvt = n;
476  ncvt = n;
477  vt.setbounds(1, nrvt, 1, ncvt);
478  }
479 
480  //
481  // M much larger than N
482  // Use bidiagonal reduction with QR-decomposition
483  //
484  if( m>amp::ampf<Precision>("1.6")*n )
485  {
486  if( uneeded==0 )
487  {
488 
489  //
490  // No left singular vectors to be computed
491  //
492  qr::qrdecomposition<Precision>(a, m, n, tau);
493  for(i=2; i<=n; i++)
494  {
495  for(j=1; j<=i-1; j++)
496  {
497  a(i,j) = 0;
498  }
499  }
500  bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
501  bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
502  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper, w, e);
503  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
504  return result;
505  }
506  else
507  {
508 
509  //
510  // Left singular vectors (may be full matrix U) to be computed
511  //
512  qr::qrdecomposition<Precision>(a, m, n, tau);
513  qr::unpackqfromqr<Precision>(a, m, n, tau, ncu, u);
514  for(i=2; i<=n; i++)
515  {
516  for(j=1; j<=i-1; j++)
517  {
518  a(i,j) = 0;
519  }
520  }
521  bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
522  bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
523  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper, w, e);
524  if( additionalmemory<1 )
525  {
526 
527  //
528  // No additional memory can be used
529  //
530  bidiagonal::multiplybyqfrombidiagonal<Precision>(a, n, n, tauq, u, m, n, true, false);
531  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
532  }
533  else
534  {
535 
536  //
537  // Large U. Transforming intermediate matrix T2
538  //
539  work.setbounds(1, ap::maxint(m, n));
540  bidiagonal::unpackqfrombidiagonal<Precision>(a, n, n, tauq, n, t2);
541  blas::copymatrix<Precision>(u, 1, m, 1, n, a, 1, m, 1, n);
542  blas::inplacetranspose<Precision>(t2, 1, n, 1, n, work);
543  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
544  blas::matrixmatrixmultiply<Precision>(a, 1, m, 1, n, false, t2, 1, n, 1, n, true, amp::ampf<Precision>("1.0"), u, 1, m, 1, n, amp::ampf<Precision>("0.0"), work);
545  }
546  return result;
547  }
548  }
549 
550  //
551  // N much larger than M
552  // Use bidiagonal reduction with LQ-decomposition
553  //
554  if( n>amp::ampf<Precision>("1.6")*m )
555  {
556  if( vtneeded==0 )
557  {
558 
559  //
560  // No right singular vectors to be computed
561  //
562  lq::lqdecomposition<Precision>(a, m, n, tau);
563  for(i=1; i<=m-1; i++)
564  {
565  for(j=i+1; j<=m; j++)
566  {
567  a(i,j) = 0;
568  }
569  }
570  bidiagonal::tobidiagonal<Precision>(a, m, m, tauq, taup);
571  bidiagonal::unpackqfrombidiagonal<Precision>(a, m, m, tauq, ncu, u);
572  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, m, isupper, w, e);
573  work.setbounds(1, m);
574  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
575  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
576  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
577  return result;
578  }
579  else
580  {
581 
582  //
583  // Right singular vectors (may be full matrix VT) to be computed
584  //
585  lq::lqdecomposition<Precision>(a, m, n, tau);
586  lq::unpackqfromlq<Precision>(a, m, n, tau, nrvt, vt);
587  for(i=1; i<=m-1; i++)
588  {
589  for(j=i+1; j<=m; j++)
590  {
591  a(i,j) = 0;
592  }
593  }
594  bidiagonal::tobidiagonal<Precision>(a, m, m, tauq, taup);
595  bidiagonal::unpackqfrombidiagonal<Precision>(a, m, m, tauq, ncu, u);
596  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, m, isupper, w, e);
597  work.setbounds(1, ap::maxint(m, n));
598  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
599  if( additionalmemory<1 )
600  {
601 
602  //
603  // No additional memory available
604  //
605  bidiagonal::multiplybypfrombidiagonal<Precision>(a, m, m, taup, vt, m, n, false, true);
606  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, n);
607  }
608  else
609  {
610 
611  //
612  // Large VT. Transforming intermediate matrix T2
613  //
614  bidiagonal::unpackptfrombidiagonal<Precision>(a, m, m, taup, m, t2);
615  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, t2, m);
616  blas::copymatrix<Precision>(vt, 1, m, 1, n, a, 1, m, 1, n);
617  blas::matrixmatrixmultiply<Precision>(t2, 1, m, 1, m, false, a, 1, m, 1, n, false, amp::ampf<Precision>("1.0"), vt, 1, m, 1, n, amp::ampf<Precision>("0.0"), work);
618  }
619  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
620  return result;
621  }
622  }
623 
624  //
625  // M<=N
626  // We can use inplace transposition of U to get rid of columnwise operations
627  //
628  if( m<=n )
629  {
630  bidiagonal::tobidiagonal<Precision>(a, m, n, tauq, taup);
631  bidiagonal::unpackqfrombidiagonal<Precision>(a, m, n, tauq, ncu, u);
632  bidiagonal::unpackptfrombidiagonal<Precision>(a, m, n, taup, nrvt, vt);
633  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, n, isupper, w, e);
634  work.setbounds(1, m);
635  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
636  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
637  blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
638  return result;
639  }
640 
641  //
642  // Simple bidiagonal reduction
643  //
644  bidiagonal::tobidiagonal<Precision>(a, m, n, tauq, taup);
645  bidiagonal::unpackqfrombidiagonal<Precision>(a, m, n, tauq, ncu, u);
646  bidiagonal::unpackptfrombidiagonal<Precision>(a, m, n, taup, nrvt, vt);
647  bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, n, isupper, w, e);
648  if( additionalmemory<2 || uneeded==0 )
649  {
650 
651  //
652  // We cant use additional memory or there is no need in such operations
653  //
654  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
655  }
656  else
657  {
658 
659  //
660  // We can use additional memory
661  //
662  t2.setbounds(1, minmn, 1, m);
663  blas::copyandtranspose<Precision>(u, 1, m, 1, minmn, t2, 1, minmn, 1, m);
664  result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
665  blas::copyandtranspose<Precision>(t2, 1, minmn, 1, m, u, 1, m, 1, minmn);
666  }
667  return result;
668  }
669 } // namespace
670 
671 #endif
int m
Definition: cfEzgcd.cc:121
int i
Definition: cfEzgcd.cc:125
void tau(int **points, int sizePoints, int k)
Definition: amp.h:83
static void make_assertion(bool bClause)
Definition: ap.h:49
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
return result
Definition: facAbsBiFact.cc:76
const CanonicalForm & w
Definition: facAbsFact.cc:55
int j
Definition: facHensel.cc:105
int maxint(int m1, int m2)
Definition: ap.cpp:162
int minint(int m1, int m2)
Definition: ap.cpp:167
Definition: svd.h:64
bool rmatrixsvd(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition: svd.h:140
bool svddecomposition(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition: svd.h:408