source: git/Singular/calcSVD.cc @ f5d2647

spielwiese
Last change on this file since f5d2647 was a9c298, checked in by Hans Schoenemann <hannes@…>, 10 years ago
format stuff
  • Property mode set to 100644
File size: 3.0 KB
Line 
1#include <stdio.h>
2#ifdef HAVE_CONFIG_H
3#include "singularconfig.h"
4#endif /* HAVE_CONFIG_H */
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#include <kernel/febase.h>
15
16template class std::vector< amp::mpfr_record* >;
17
18poly p_svdInit(char *s)
19{
20  poly p=pInit();
21  currRing->cf->nRead(s,&pGetCoeff(p));
22  return p;
23}
24
25/*************************************************************************
26Testing SVD decomposition subroutine
27*************************************************************************/
28lists testsvd(matrix M)
29{
30
31    const unsigned int Precision=300;
32
33    int max_i=M->nrows;
34    int max_j=M->ncols;
35    ap::template_2d_array< amp::ampf<Precision> > a;
36    int m;
37    int n;
38    int maxmn;
39    int i;
40    int j;
41    int gpass;
42    int pass;
43    bool waserrors;
44    bool wsorted;
45    bool wfailed;
46    amp::ampf<Precision> materr;
47    amp::ampf<Precision> orterr;
48    amp::ampf<Precision> othererr;
49    amp::ampf<Precision> threshold;
50    amp::ampf<Precision> failthreshold;
51    amp::ampf<Precision> failr;
52
53
54    materr = 0;
55    orterr = 0;
56    othererr = 0;
57    wsorted = true;
58    wfailed = false;
59    waserrors = false;
60    maxmn = 50;
61    threshold = 5*100*amp::ampf<Precision>::getAlgoPascalEpsilon();
62    failthreshold = amp::ampf<Precision>("5.0E-3");
63    a.setbounds(1, max_i, 1, max_j);
64
65
66        //
67        // fill matrix a entries from M
68        //
69        for(i=1; i<=max_i; i++)
70        {
71            for(j=1; j<=max_j; j++)
72            {
73                char *str=pString(MATELEM(M,i,j));
74                Print(" to svd:%d,%d=%s\n",i,j,str);
75                a(i,j) = amp::ampf<Precision>(str);
76            }
77        }
78        //testsvdproblem<Precision>(a, max_i, max_j, materr, orterr, othererr, wsorted, wfailed);
79    ap::template_2d_array< amp::ampf<Precision> > u;
80    ap::template_2d_array< amp::ampf<Precision> > vt;
81    ap::template_1d_array< amp::ampf<Precision> > w;
82    svd::svddecomposition<Precision>(a, max_i, max_j, 2, 2, 2, w, u, vt);       
83    matrix Mu,Mw,Mvt;
84    Mu=mpNew(max_i,max_i);
85    Mw=mpNew(max_i,max_j);
86    Mvt=mpNew(max_j,max_j);
87    for(i=1;i<=max_i;i++)
88    {
89        for(j=1;j<=max_i;j++)
90        {
91                MATELEM(Mu,i,j)=p_svdInit(u(i,j).toString());
92//                      Print(" after svd:%d,%d=%s\n",i,j,u(i,j).toString());
93        }
94    }
95    for(i=1;i<=si_min(max_i,max_j);i++)
96    {
97                MATELEM(Mw,i,i)=p_svdInit(w(i).toString());
98//Print(" after svd:%d,%d=%s\n",i,w(i).toString());                     
99    }
100    for(i=1;i<=max_j;i++)
101    {
102        for(j=1;j<=max_j;j++)
103        {
104                MATELEM(Mvt,i,j)=p_svdInit(vt(i,j).toString());
105//Print(" after svd:%d,%d=%s\n",i,j,vt(i,j).toString());                       
106        }
107    }
108    lists L=(lists)omAlloc(sizeof(slists));
109    L->Init(3);
110    L->m[0].rtyp=MATRIX_CMD;
111    L->m[1].rtyp=MATRIX_CMD;
112    L->m[2].rtyp=MATRIX_CMD;
113    L->m[0].data=(char*)Mu;
114    L->m[1].data=(char*)Mw;
115    L->m[2].data=(char*)Mvt;
116    return L;
117}
118
119#endif
Note: See TracBrowser for help on using the repository browser.