My Project
Loading...
Searching...
No Matches
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;
151 ap::template_1d_array< amp::ampf<Precision> > tauq;
152 ap::template_1d_array< amp::ampf<Precision> > taup;
153 ap::template_1d_array< amp::ampf<Precision> > tau;
154 ap::template_1d_array< amp::ampf<Precision> > e;
155 ap::template_1d_array< amp::ampf<Precision> > work;
156 ap::template_2d_array< amp::ampf<Precision> > t2;
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;
166 amp::ampf<Precision> sm;
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 }
m
int m
Definition: cfEzgcd.cc:128
i
int i
Definition: cfEzgcd.cc:132
tau
void tau(int **points, int sizePoints, int k)
Definition: cfNewtonPolygon.cc:461
amp::ampf
Definition: amp.h:82
ap::ap_error::make_assertion
static void make_assertion(bool bClause)
Definition: ap.h:49
ap::template_1d_array
Definition: ap.h:657
ap::template_1d_array::setbounds
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
ap::template_2d_array
Definition: ap.h:807
ap::template_2d_array::setbounds
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
result
return result
Definition: facAbsBiFact.cc:75
w
const CanonicalForm & w
Definition: facAbsFact.cc:51
j
int j
Definition: facHensel.cc:110
ap::maxint
int maxint(int m1, int m2)
Definition: ap.cpp:162
ap::minint
int minint(int m1, int m2)
Definition: ap.cpp:167

◆ 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;
419 ap::template_1d_array< amp::ampf<Precision> > tauq;
420 ap::template_1d_array< amp::ampf<Precision> > taup;
421 ap::template_1d_array< amp::ampf<Precision> > tau;
422 ap::template_1d_array< amp::ampf<Precision> > e;
423 ap::template_1d_array< amp::ampf<Precision> > work;
424 ap::template_2d_array< amp::ampf<Precision> > t2;
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;
434 amp::ampf<Precision> sm;
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 }
Generated on Mon Feb 27 2023 10:53:54 for My Project by   doxygen 1.9.5 for Singular