Loading...
Searching...
No Matches
testsvdunit.h
Go to the documentation of this file.
#ifndef _testsvdunit_h
#define _testsvdunit_h
#include <stdio.h>
#include "ap.h"
#include "amp.h"
#include "reflections.h"
#include "bidiagonal.h"
#include "qr.h"
#include "lq.h"
#include "blas.h"
#include "rotations.h"
#include "bdsvd.h"
#include "svd.h"
namespace testsvdunit
{
template<unsigned int Precision>
bool testsvd(bool silent);
template<unsigned int Precision>
void fillsparsea(ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
amp::ampf<Precision> sparcity);
template<unsigned int Precision>
void getsvderror(const ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
const ap::template_2d_array< amp::ampf<Precision> >& u,
const ap::template_1d_array< amp::ampf<Precision> >& w,
const ap::template_2d_array< amp::ampf<Precision> >& vt,
amp::ampf<Precision>& materr,
amp::ampf<Precision>& orterr,
bool& wsorted);
template<unsigned int Precision>
void testsvdproblem(const ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
amp::ampf<Precision>& materr,
amp::ampf<Precision>& orterr,
amp::ampf<Precision>& othererr,
bool& wsorted,
bool& wfailed);
template<unsigned int Precision>
bool testsvdunit_test_silent();
template<unsigned int Precision>
bool testsvdunit_test();
STATIC_VAR int failcount;
STATIC_VAR int succcount;
/*************************************************************************
Testing SVD decomposition subroutine
*************************************************************************/
template<unsigned int Precision>
bool testsvd(bool silent)
{
bool result;
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;
failcount = 0;
succcount = 0;
materr = 0;
orterr = 0;
othererr = 0;
wsorted = true;
wfailed = false;
waserrors = false;
maxmn = 30;
threshold = 5*100*amp::ampf<Precision>::getAlgoPascalEpsilon();
failthreshold = amp::ampf<Precision>("5.0E-3");
a.setbounds(0, maxmn-1, 0, maxmn-1);
//
// TODO: div by zero fail, convergence fail
//
for(gpass=1; gpass<=1; gpass++)
{
//
// zero matrix, several cases
//
for(i=0; i<=maxmn-1; i++)
{
for(j=0; j<=maxmn-1; j++)
{
a(i,j) = 0;
}
}
for(i=1; i<=ap::minint(5, maxmn); i++)
{
for(j=1; j<=ap::minint(5, maxmn); j++)
{
testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
}
}
//
// Long dense matrix
//
for(i=0; i<=maxmn-1; i++)
{
for(j=0; j<=ap::minint(5, maxmn)-1; j++)
{
a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
}
}
for(i=1; i<=maxmn; i++)
{
for(j=1; j<=ap::minint(5, maxmn); j++)
{
testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
}
}
for(i=0; i<=ap::minint(5, maxmn)-1; i++)
{
for(j=0; j<=maxmn-1; j++)
{
a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
}
}
for(i=1; i<=ap::minint(5, maxmn); i++)
{
for(j=1; j<=maxmn; j++)
{
testsvdproblem<Precision>(a, i, j, materr, orterr, othererr, wsorted, wfailed);
}
}
//
// Dense matrices
//
for(m=1; m<=ap::minint(10, maxmn); m++)
{
for(n=1; n<=ap::minint(10, maxmn); n++)
{
for(i=0; i<=m-1; i++)
{
for(j=0; j<=n-1; j++)
{
a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
}
}
testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
}
}
//
// Sparse matrices, very sparse matrices, incredible sparse matrices
//
for(m=1; m<=10; m++)
{
for(n=1; n<=10; n++)
{
for(pass=1; pass<=2; pass++)
{
fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.8"));
testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.9"));
testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
fillsparsea<Precision>(a, m, n, amp::ampf<Precision>("0.95"));
testsvdproblem<Precision>(a, m, n, materr, orterr, othererr, wsorted, wfailed);
}
}
}
}
//
// report
//
failr = amp::ampf<Precision>(failcount)/(amp::ampf<Precision>(succcount+failcount));
waserrors = materr>threshold || orterr>threshold || othererr>threshold || !wsorted || failr>failthreshold;
if( !silent )
{
printf("TESTING SVD DECOMPOSITION\n");
printf("SVD decomposition error: %5.3le\n",
double(amp::ampf<Precision>(materr).toDouble()));
printf("SVD orthogonality error: %5.3le\n",
double(amp::ampf<Precision>(orterr).toDouble()));
printf("SVD with different parameters error: %5.3le\n",
double(amp::ampf<Precision>(othererr).toDouble()));
printf("Singular values order: ");
if( wsorted )
{
printf("OK\n");
}
else
{
printf("FAILED\n");
}
printf("Always converged: ");
if( !wfailed )
{
printf("YES\n");
}
else
{
printf("NO\n");
printf("Fail ratio: %5.3lf\n",
double(amp::ampf<Precision>(failr).toDouble()));
}
printf("Threshold: %5.3le\n",
double(amp::ampf<Precision>(threshold).toDouble()));
if( waserrors )
{
printf("TEST FAILED\n");
}
else
{
printf("TEST PASSED\n");
}
printf("\n\n");
}
result = !waserrors;
return result;
}
template<unsigned int Precision>
void fillsparsea(ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
amp::ampf<Precision> sparcity)
{
int i;
int j;
for(i=0; i<=m-1; i++)
{
for(j=0; j<=n-1; j++)
{
if( amp::ampf<Precision>::getRandom()>=sparcity )
{
a(i,j) = 2*amp::ampf<Precision>::getRandom()-1;
}
else
{
a(i,j) = 0;
}
}
}
}
template<unsigned int Precision>
void getsvderror(const ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
const ap::template_2d_array< amp::ampf<Precision> >& u,
const ap::template_1d_array< amp::ampf<Precision> >& w,
const ap::template_2d_array< amp::ampf<Precision> >& vt,
amp::ampf<Precision>& materr,
amp::ampf<Precision>& orterr,
bool& wsorted)
{
int i;
int j;
int k;
int minmn;
amp::ampf<Precision> locerr;
amp::ampf<Precision> sm;
minmn = ap::minint(m, n);
//
// decomposition error
//
locerr = 0;
for(i=0; i<=m-1; i++)
{
for(j=0; j<=n-1; j++)
{
sm = 0;
for(k=0; k<=minmn-1; k++)
{
sm = sm+w(k)*u(i,k)*vt(k,j);
}
locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(a(i,j)-sm));
}
}
materr = amp::maximum<Precision>(materr, locerr);
//
// orthogonality error
//
locerr = 0;
for(i=0; i<=minmn-1; i++)
{
for(j=i; j<=minmn-1; j++)
{
sm = ap::vdotproduct(u.getcolumn(i, 0, m-1), u.getcolumn(j, 0, m-1));
if( i!=j )
{
locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm));
}
else
{
locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm-1));
}
sm = ap::vdotproduct(vt.getrow(i, 0, n-1), vt.getrow(j, 0, n-1));
if( i!=j )
{
locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm));
}
else
{
locerr = amp::maximum<Precision>(locerr, amp::abs<Precision>(sm-1));
}
}
}
orterr = amp::maximum<Precision>(orterr, locerr);
//
// values order error
//
for(i=1; i<=minmn-1; i++)
{
if( w(i)>w(i-1) )
{
wsorted = false;
}
}
}
template<unsigned int Precision>
void testsvdproblem(const ap::template_2d_array< amp::ampf<Precision> >& a,
int m,
int n,
amp::ampf<Precision>& materr,
amp::ampf<Precision>& orterr,
amp::ampf<Precision>& othererr,
bool& wsorted,
bool& wfailed)
{
ap::template_2d_array< amp::ampf<Precision> > u;
ap::template_2d_array< amp::ampf<Precision> > vt;
ap::template_2d_array< amp::ampf<Precision> > u2;
ap::template_2d_array< amp::ampf<Precision> > vt2;
ap::template_1d_array< amp::ampf<Precision> > w;
ap::template_1d_array< amp::ampf<Precision> > w2;
int i;
int j;
int k;
int ujob;
int vtjob;
int memjob;
int ucheck;
int vtcheck;
amp::ampf<Precision> v;
amp::ampf<Precision> mx;
//
// Main SVD test
//
if( !svd::rmatrixsvd<Precision>(a, m, n, 2, 2, 2, w, u, vt) )
{
failcount = failcount+1;
wfailed = true;
return;
}
getsvderror<Precision>(a, m, n, u, w, vt, materr, orterr, wsorted);
//
// Additional SVD tests
//
for(ujob=0; ujob<=2; ujob++)
{
for(vtjob=0; vtjob<=2; vtjob++)
{
for(memjob=0; memjob<=2; memjob++)
{
if( !svd::rmatrixsvd<Precision>(a, m, n, ujob, vtjob, memjob, w2, u2, vt2) )
{
failcount = failcount+1;
wfailed = true;
return;
}
ucheck = 0;
if( ujob==1 )
{
ucheck = ap::minint(m, n);
}
if( ujob==2 )
{
ucheck = m;
}
vtcheck = 0;
if( vtjob==1 )
{
vtcheck = ap::minint(m, n);
}
if( vtjob==2 )
{
vtcheck = n;
}
for(i=0; i<=m-1; i++)
{
for(j=0; j<=ucheck-1; j++)
{
othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(u(i,j)-u2(i,j)));
}
}
for(i=0; i<=vtcheck-1; i++)
{
for(j=0; j<=n-1; j++)
{
othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(vt(i,j)-vt2(i,j)));
}
}
for(i=0; i<=ap::minint(m, n)-1; i++)
{
othererr = amp::maximum<Precision>(othererr, amp::abs<Precision>(w(i)-w2(i)));
}
}
}
}
//
// update counter
//
succcount = succcount+1;
}
/*************************************************************************
Silent unit test
*************************************************************************/
template<unsigned int Precision>
bool testsvdunit_test_silent()
{
bool result;
result = testsvd<Precision>(true);
return result;
}
/*************************************************************************
Unit test
*************************************************************************/
template<unsigned int Precision>
bool testsvdunit_test()
{
bool result;
result = testsvd<Precision>(false);
return result;
}
} // namespace
#endif
amp.h
ap.h
bdsvd.h
bidiagonal.h
blas.h
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
amp::ampf
Definition amp.h:82
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition amp.h:565
amp::ampf::getRandom
static const ampf getRandom()
Definition amp.h:593
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
result
return result
Definition facAbsBiFact.cc:76
w
const CanonicalForm & w
Definition facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
STATIC_VAR
#define STATIC_VAR
Definition globaldefs.h:7
lq.h
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition amp.h:713
amp::maximum
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:722
ap::vdotproduct
T vdotproduct(const_raw_vector< T > v1, const_raw_vector< T > v2)
Definition ap.h:181
ap::minint
int minint(int m1, int m2)
Definition ap.cpp:167
svd::rmatrixsvd
bool rmatrixsvd(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:140
testsvdunit
Definition testsvdunit.h:17
testsvdunit::fillsparsea
void fillsparsea(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, amp::ampf< Precision > sparcity)
Definition testsvdunit.h:239
testsvdunit::testsvdunit_test
bool testsvdunit_test()
Definition testsvdunit.h:467
testsvdunit::testsvd
bool testsvd(bool silent)
Definition testsvdunit.h:58
testsvdunit::getsvderror
void getsvderror(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_2d_array< amp::ampf< Precision > > &u, const ap::template_1d_array< amp::ampf< Precision > > &w, const ap::template_2d_array< amp::ampf< Precision > > &vt, amp::ampf< Precision > &materr, amp::ampf< Precision > &orterr, bool &wsorted)
Definition testsvdunit.h:266
testsvdunit::testsvdunit_test_silent
bool testsvdunit_test_silent()
Definition testsvdunit.h:453
testsvdunit::failcount
STATIC_VAR int failcount
Definition testsvdunit.h:50
testsvdunit::testsvdproblem
void testsvdproblem(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, amp::ampf< Precision > &materr, amp::ampf< Precision > &orterr, amp::ampf< Precision > &othererr, bool &wsorted, bool &wfailed)
Definition testsvdunit.h:348
testsvdunit::succcount
STATIC_VAR int succcount
Definition testsvdunit.h:51
qr.h
reflections.h
rotations.h
svd.h