Functions
svd Namespace Reference

Functions

template<unsigned int Precision>
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)
 
template<unsigned int Precision>
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)
 

Function Documentation

◆ rmatrixsvd()

template<unsigned int Precision>
bool svd::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 at line 140 of file svd.h.

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  }
int j
Definition: facHensel.cc:105
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
static void make_assertion(bool bClause)
Definition: ap.h:49
int m
Definition: cfEzgcd.cc:121
void tau(int **points, int sizePoints, int k)
int i
Definition: cfEzgcd.cc:125
int minint(int m1, int m2)
Definition: ap.cpp:167
int maxint(int m1, int m2)
Definition: ap.cpp:162
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
const CanonicalForm & w
Definition: facAbsFact.cc:55
Definition: amp.h:81
return result
Definition: facAbsBiFact.cc:76

◆ svddecomposition()

template<unsigned int Precision>
bool svd::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 at line 408 of file svd.h.

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  }
int j
Definition: facHensel.cc:105
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
static void make_assertion(bool bClause)
Definition: ap.h:49
int m
Definition: cfEzgcd.cc:121
void tau(int **points, int sizePoints, int k)
int i
Definition: cfEzgcd.cc:125
int minint(int m1, int m2)
Definition: ap.cpp:167
int maxint(int m1, int m2)
Definition: ap.cpp:162
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
const CanonicalForm & w
Definition: facAbsFact.cc:55
Definition: amp.h:81
return result
Definition: facAbsBiFact.cc:76