My Project
Loading...
Searching...
No Matches
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
15template class std::vector< amp::mpfr_record* >;
16
17poly p_svdInit(char *s)
18{
19 poly p=pInit();
20 n_Read(s, &pGetCoeff(p), currRing->cf);
21 return p;
22}
23
24/*************************************************************************
25Testing 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;
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
static int si_min(const int a, const int b)
Definition: auxiliary.h:125
std::string toString(const gfan::ZCone *const c)
Definition: bbcone.cc:27
poly p_svdInit(char *s)
Definition: calcSVD.cc:17
lists testsvd(matrix M)
Definition: calcSVD.cc:27
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int p
Definition: cfModGcd.cc:4078
Definition: amp.h:82
static const ampf getAlgoPascalEpsilon()
Definition: amp.h:565
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
int rtyp
Definition: subexpr.h:91
void * data
Definition: subexpr.h:88
Definition: lists.h:24
sleftv * m
Definition: lists.h:46
INLINE_THIS void Init(int l=0)
static FORCE_INLINE const 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:595
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm & w
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
@ MATRIX_CMD
Definition: grammar.cc:286
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
slists * lists
Definition: mpr_numeric.h:146
#define omAlloc(size)
Definition: omAllocDecl.h:210
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
char * pString(poly p)
Definition: polys.h:306
#define M
Definition: sirandom.c:25