Loading...
Searching...
No Matches
calcSVD.cc
Go to the documentation of this file.
#include <stdio.h>
#include "kernel/mod2.h"
#ifdef HAVE_SVD
#include "Singular/svd_si.h"
#include "kernel/structs.h"
#include "kernel/polys.h"
#include "polys/matpol.h"
#include "Singular/lists.h"
template class std::vector< amp::mpfr_record* >;
poly p_svdInit(char *s)
{
poly p=pInit();
n_Read(s, &pGetCoeff(p), currRing->cf);
return p;
}
/*************************************************************************
Testing SVD decomposition subroutine
*************************************************************************/
lists testsvd(matrix M)
{
const unsigned int Precision=300;
int max_i=M->nrows;
int max_j=M->ncols;
ap::template_2d_array< amp::ampf<Precision> > a;
int m;
int n;
int maxmn;
int i;
int j;
int gpass;
int pass;
bool waserrors;
bool wsorted;
bool wfailed;
amp::ampf<Precision> materr;
amp::ampf<Precision> orterr;
amp::ampf<Precision> othererr;
amp::ampf<Precision> threshold;
amp::ampf<Precision> failthreshold;
amp::ampf<Precision> failr;
materr = 0;
orterr = 0;
othererr = 0;
wsorted = true;
wfailed = false;
waserrors = false;
maxmn = 50;
threshold = 5*100*amp::ampf<Precision>::getAlgoPascalEpsilon();
failthreshold = amp::ampf<Precision>("5.0E-3");
a.setbounds(1, max_i, 1, max_j);
//
// fill matrix a entries from M
//
for(i=1; i<=max_i; i++)
{
for(j=1; j<=max_j; j++)
{
char *str=pString(MATELEM(M,i,j));
Print(" to svd:%d,%d=%s\n",i,j,str);
a(i,j) = amp::ampf<Precision>(str);
}
}
//testsvdproblem<Precision>(a, max_i, max_j, materr, orterr, othererr, wsorted, wfailed);
ap::template_2d_array< amp::ampf<Precision> > u;
ap::template_2d_array< amp::ampf<Precision> > vt;
ap::template_1d_array< amp::ampf<Precision> > w;
svd::svddecomposition<Precision>(a, max_i, max_j, 2, 2, 2, w, u, vt);
matrix Mu,Mw,Mvt;
Mu=mpNew(max_i,max_i);
Mw=mpNew(max_i,max_j);
Mvt=mpNew(max_j,max_j);
for(i=1;i<=max_i;i++)
{
for(j=1;j<=max_i;j++)
{
MATELEM(Mu,i,j)=p_svdInit(u(i,j).toString());
// Print(" after svd:%d,%d=%s\n",i,j,u(i,j).toString());
}
}
for(i=1;i<=si_min(max_i,max_j);i++)
{
MATELEM(Mw,i,i)=p_svdInit(w(i).toString());
//Print(" after svd:%d,%d=%s\n",i,w(i).toString());
}
for(i=1;i<=max_j;i++)
{
for(j=1;j<=max_j;j++)
{
MATELEM(Mvt,i,j)=p_svdInit(vt(i,j).toString());
//Print(" after svd:%d,%d=%s\n",i,j,vt(i,j).toString());
}
}
lists L=(lists)omAlloc(sizeof(slists));
L->Init(3);
L->m[0].rtyp=MATRIX_CMD;
L->m[1].rtyp=MATRIX_CMD;
L->m[2].rtyp=MATRIX_CMD;
L->m[0].data=(char*)Mu;
L->m[1].data=(char*)Mw;
L->m[2].data=(char*)Mvt;
return L;
}
#endif
si_min
static int si_min(const int a, const int b)
Definition auxiliary.h:126
toString
std::string toString(const gfan::ZCone *const c)
Definition bbcone.cc:27
p_svdInit
poly p_svdInit(char *s)
Definition calcSVD.cc:17
testsvd
lists testsvd(matrix M)
Definition calcSVD.cc:27
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
p
int p
Definition cfModGcd.cc:4086
amp::ampf
Definition amp.h:82
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition amp.h:565
ap::template_1d_array
Definition ap.h:657
ap::template_2d_array
Definition ap.h:807
ap::template_2d_array::setbounds
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition ap.h:890
sleftv::rtyp
int rtyp
Definition subexpr.h:91
sleftv::data
void * data
Definition subexpr.h:88
slists
Definition lists.h:24
slists::m
sleftv * m
Definition lists.h:46
slists::Init
INLINE_THIS void Init(int l=0)
n_Read
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:601
Print
#define Print
Definition emacs.cc:80
s
const CanonicalForm int s
Definition facAbsFact.cc:51
w
const CanonicalForm & w
Definition facAbsFact.cc:51
j
int j
Definition facHensel.cc:110
MATRIX_CMD
@ MATRIX_CMD
Definition grammar.cc:287
lists.h
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition matpol.cc:37
matpol.h
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
matrix
ip_smatrix * matrix
Definition matpol.h:43
mod2.h
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:44
lists
slists * lists
Definition mpr_numeric.h:146
svd::svddecomposition
bool svddecomposition(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition svd.h:408
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
polys.h
Compatibility layer for legacy polynomial operations (over currRing).
pInit
#define pInit()
allocates a new monomial and initializes everything to 0
Definition polys.h:62
pString
char * pString(poly p)
Definition polys.h:307
M
#define M
Definition sirandom.c:25
structs.h
svd_si.h