1 #include <stdio.h>
2
3
4
5 #include "kernel/mod2.h"
6
7 #ifdef HAVE_SVD
8
9 #include "Singular/svd_si.h"
10 #include "kernel/structs.h"
11 #include "kernel/polys.h"
12 #include "polys/matpol.h"
13 #include "Singular/lists.h"
14
15 template class std::vector< amp::mpfr_record* >;
16
p_svdInit(char * s)17 poly p_svdInit(char *s)
18 {
19 poly p=pInit();
20 n_Read(s, &pGetCoeff(p), currRing->cf);
21 return p;
22 }
23
24 /*************************************************************************
25 Testing SVD decomposition subroutine
26 *************************************************************************/
testsvd(matrix M)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);
85 Mvt=mpNew(max_j,max_j);
86 for(i=1;i<=max_i;i++)
87 {
88 for(j=1;j<=max_i;j++)
89 {
90 MATELEM(Mu,i,j)=p_svdInit(u(i,j).toString());
91 // Print(" after svd:%d,%d=%s\n",i,j,u(i,j).toString());
92 }
93 }
94 for(i=1;i<=si_min(max_i,max_j);i++)
95 {
96 MATELEM(Mw,i,i)=p_svdInit(w(i).toString());
97 //Print(" after svd:%d,%d=%s\n",i,w(i).toString());
98 }
99 for(i=1;i<=max_j;i++)
100 {
101 for(j=1;j<=max_j;j++)
102 {
103 MATELEM(Mvt,i,j)=p_svdInit(vt(i,j).toString());
104 //Print(" after svd:%d,%d=%s\n",i,j,vt(i,j).toString());
105 }
106 }
107 lists L=(lists)omAlloc(sizeof(slists));
108 L->Init(3);
109 L->m[0].rtyp=MATRIX_CMD;
110 L->m[1].rtyp=MATRIX_CMD;
111 L->m[2].rtyp=MATRIX_CMD;
112 L->m[0].data=(char*)Mu;
113 L->m[1].data=(char*)Mw;
114 L->m[2].data=(char*)Mvt;
115 return L;
116 }
117
118 #endif
119