My Project  debian-1:4.1.1-p2+ds-4
calcSVD.cc
Go to the documentation of this file.
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 
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 *************************************************************************/
28 {
29 
30  const unsigned int Precision=300;
31 
32  int max_i=M->nrows;
33  int max_j=M->ncols;
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;
51 
52 
53  materr = 0;
54  orterr = 0;
55  othererr = 0;
56  wsorted = true;
57  wfailed = false;
58  waserrors = false;
59  maxmn = 50;
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);
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
si_min
static int si_min(const int a, const int b)
Definition: auxiliary.h:139
toString
std::string toString(const gfan::ZCone *const c)
Definition: bbcone.cc:27
ip_smatrix
Definition: matpol.h:13
j
int j
Definition: facHensel.cc:105
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition: amp.h:566
MATELEM
#define MATELEM(mat, i, j)
Definition: matpol.h:28
ap::template_1d_array
Definition: ap.h:641
ap::template_2d_array
Definition: ap.h:790
lists.h
svd_si.h
amp::ampf
Definition: amp.h:82
polys.h
ap::template_2d_array::setbounds
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:875
n_Read
static const FORCE_INLINE char * n_Read(const char *s, number *a, const coeffs r)
!!! Recommendation: This method is too cryptic to be part of the user- !!! interface....
Definition: coeffs.h:598
w
const CanonicalForm & w
Definition: facAbsFact.cc:55
pString
char * pString(poly p)
Definition: polys.h:280
MATRIX_CMD
Definition: grammar.cc:285
testsvd
lists testsvd(matrix M)
Definition: calcSVD.cc:26
currRing
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
p_svdInit
poly p_svdInit(char *s)
Definition: calcSVD.cc:17
i
int i
Definition: cfEzgcd.cc:125
matpol.h
M
#define M
Definition: sirandom.c:24
structs.h
mod2.h
pInit
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:59
sleftv::data
void * data
Definition: subexpr.h:87
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:208
slists::m
sleftv * m
Definition: lists.h:44
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:35
slists
Definition: lists.h:21
Print
#define Print
Definition: emacs.cc:79
m
int m
Definition: cfEzgcd.cc:121
sleftv::rtyp
int rtyp
Definition: subexpr.h:90
lists
slists * lists
Definition: mpr_numeric.h:145
slists::Init
INLINE_THIS void Init(int l=0)
p
int p
Definition: cfModGcd.cc:4019
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:48