Testing SVD decomposition subroutine *************************************************************************/
| 25 | Testing SVD decomposition subroutine |
| 26 | *************************************************************************/ |
| 27 | lists testsvd(matrix M) |
| 28 | { |
| 29 | |
| 30 | const unsigned int Precision=300; |
| 31 | |
| 32 | int max_i=M->nrows; |
| 33 | int max_j=M->ncols; |
| 34 | ap::template_2d_array< amp::ampf<Precision> > a; |
| 35 | int m; |
| 36 | int n; |
| 37 | int maxmn; |
| 38 | int i; |
| 39 | int j; |
| 40 | int gpass; |
| 41 | int pass; |
| 42 | bool waserrors; |
| 43 | bool wsorted; |
| 44 | bool wfailed; |
| 45 | amp::ampf<Precision> materr; |
| 46 | amp::ampf<Precision> orterr; |
| 47 | amp::ampf<Precision> othererr; |
| 48 | amp::ampf<Precision> threshold; |
| 49 | amp::ampf<Precision> failthreshold; |
| 50 | amp::ampf<Precision> failr; |
| 51 | |
| 52 | |
| 53 | materr = 0; |
| 54 | orterr = 0; |
| 55 | othererr = 0; |
| 56 | wsorted = true; |
| 57 | wfailed = false; |
| 58 | waserrors = false; |
| 59 | maxmn = 50; |
| 60 | threshold = 5*100*amp::ampf<Precision>::getAlgoPascalEpsilon(); |
| 61 | failthreshold = amp::ampf<Precision>("5.0E-3"); |
| 62 | a.setbounds(1, max_i, 1, max_j); |
| 63 | |
| 64 | |
| 65 | // |
| 66 | // fill matrix a entries from M |
| 67 | // |
| 68 | for(i=1; i<=max_i; i++) |
| 69 | { |
| 70 | for(j=1; j<=max_j; j++) |
| 71 | { |
| 72 | char *str=pString(MATELEM(M,i,j)); |
| 73 | Print(" to svd:%d,%d=%s\n",i,j,str); |
| 74 | a(i,j) = amp::ampf<Precision>(str); |
| 75 | } |
| 76 | } |
| 77 | //testsvdproblem<Precision>(a, max_i, max_j, materr, orterr, othererr, wsorted, wfailed); |
| 78 | ap::template_2d_array< amp::ampf<Precision> > u; |
| 79 | ap::template_2d_array< amp::ampf<Precision> > vt; |
| 80 | ap::template_1d_array< amp::ampf<Precision> > w; |
| 81 | svd::svddecomposition<Precision>(a, max_i, max_j, 2, 2, 2, w, u, vt); |
| 82 | matrix Mu,Mw,Mvt; |
| 83 | Mu=mpNew(max_i,max_i); |
| 84 | Mw=mpNew(max_i,max_j); |