Dynare tensor library tests: further modernizations

time-shift
Sébastien Villemot 2019-02-13 16:32:23 +01:00
parent 5cbc34e9de
commit 579be3c5e2
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
3 changed files with 48 additions and 51 deletions

View File

@ -3,52 +3,49 @@
#include "factory.hh" #include "factory.hh"
#include <cstdlib>
#include <cmath>
void void
Factory::init(const Symmetry &s, const IntSequence &nvs) Factory::init(const Symmetry &s, const IntSequence &nvs)
{ {
IntSequence sym(s); IntSequence sym(s);
long int seed = sym[0]; decltype(mtgen)::result_type seed = sym[0];
seed = 256*seed + nvs[0]; seed = 256*seed + nvs[0];
if (sym.size() > 1) if (sym.size() > 1)
seed = 256*seed + sym[1]; seed = 256*seed + sym[1];
if (nvs.size() > 1) if (nvs.size() > 1)
seed = 256*seed + nvs[0]; seed = 256*seed + nvs[0];
srand48(seed); mtgen.seed(seed);
} }
void void
Factory::init(int dim, int nv) Factory::init(int dim, int nv)
{ {
long int seed = dim; decltype(mtgen)::result_type seed = dim;
seed = 256*seed + nv; seed = 256*seed + nv;
srand48(seed); mtgen.seed(seed);
} }
double double
Factory::get() const Factory::get()
{ {
return 1.0*(drand48()-0.5); return dis(mtgen)-0.5;
} }
void void
Factory::fillMatrix(TwoDMatrix &m) const Factory::fillMatrix(TwoDMatrix &m)
{ {
Vector &d = m.getData(); Vector &d = m.getData();
for (int i = 0; i < d.length(); i++) for (int i = 0; i < d.length(); i++)
d[i] = get(); d[i] = get();
} }
Vector * Vector
Factory::makeVector(int n) Factory::makeVector(int n)
{ {
init(n, n*n); init(n, n*n);
auto *v = new Vector(n); Vector v(n);
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
(*v)[i] = get(); v[i] = get();
return v; return v;
} }

View File

@ -4,6 +4,8 @@
#ifndef FACTORY_H #ifndef FACTORY_H
#define FACTORY_H #define FACTORY_H
#include <random>
#include "symmetry.hh" #include "symmetry.hh"
#include "int_sequence.hh" #include "int_sequence.hh"
#include "twod_matrix.hh" #include "twod_matrix.hh"
@ -13,11 +15,14 @@
class Factory class Factory
{ {
std::mt19937 mtgen;
std::uniform_real_distribution<> dis;
void init(const Symmetry &s, const IntSequence &nvs); void init(const Symmetry &s, const IntSequence &nvs);
void init(int dim, int nv); void init(int dim, int nv);
void fillMatrix(TwoDMatrix &m) const; void fillMatrix(TwoDMatrix &m);
public: public:
double get() const; double get();
// this can be used with UGSTensor, FGSTensor // this can be used with UGSTensor, FGSTensor
template <class _Ttype> template <class _Ttype>
_Ttype * _Ttype *
@ -47,25 +52,23 @@ public:
int symnum = nvs.size(); int symnum = nvs.size();
auto *res = new _Ctype(symnum); auto *res = new _Ctype(symnum);
for (int dim = 1; dim <= maxdim; dim++) for (int dim = 1; dim <= maxdim; dim++)
{ if (symnum == 1)
if (symnum == 1) {
{ // full symmetry
// full symmetry Symmetry sym{dim};
Symmetry sym{dim}; auto *t = make<_Ttype>(r, sym, nvs);
auto *t = make<_Ttype>(r, sym, nvs); res->insert(t);
res->insert(t); }
} else
else {
{ // general symmetry
// general symmetry for (int i = 0; i <= dim; i++)
for (int i = 0; i <= dim; i++) {
{ Symmetry sym{i, dim-i};
Symmetry sym{i, dim-i}; auto *t = make<_Ttype>(r, sym, nvs);
auto *t = make<_Ttype>(r, sym, nvs); res->insert(t);
res->insert(t); }
} }
}
}
return res; return res;
} }
@ -82,7 +85,7 @@ public:
return p; return p;
} }
Vector *makeVector(int n); Vector makeVector(int n);
}; };
#endif #endif

View File

@ -431,25 +431,25 @@ bool
TestRunnable::folded_contraction(int r, int nv, int dim) TestRunnable::folded_contraction(int r, int nv, int dim)
{ {
Factory fact; Factory fact;
Vector *x = fact.makeVector(nv); Vector x{fact.makeVector(nv)};
auto *forig = fact.make<FFSTensor>(r, nv, dim); auto *forig = fact.make<FFSTensor>(r, nv, dim);
auto *f = new FFSTensor(*forig); auto *f = new FFSTensor(*forig);
clock_t ctime = clock(); clock_t ctime = clock();
for (int d = dim-1; d > 0; d--) for (int d = dim-1; d > 0; d--)
{ {
FFSTensor *fnew = new FFSTensor(*f, ConstVector(*x)); FFSTensor *fnew = new FFSTensor(*f, ConstVector(x));
delete f; delete f;
f = fnew; f = fnew;
} }
ctime = clock() - ctime; ctime = clock() - ctime;
Vector res(forig->nrows()); Vector res(forig->nrows());
res.zeros(); res.zeros();
f->multaVec(res, *x); f->multaVec(res, x);
UFSTensor u(*forig); UFSTensor u(*forig);
clock_t utime = clock(); clock_t utime = clock();
URSingleTensor ux(*x, dim); URSingleTensor ux(x, dim);
Vector v(u.nrows()); Vector v(u.nrows());
v.zeros(); v.zeros();
u.multaVec(v, ux.getData()); u.multaVec(v, ux.getData());
@ -464,7 +464,6 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
<< "\terror norm1: " << v.getNorm1() << '\n'; << "\terror norm1: " << v.getNorm1() << '\n';
delete f; delete f;
delete x;
return (v.getMax() < 1.e-10); return (v.getMax() < 1.e-10);
} }
@ -473,7 +472,7 @@ bool
TestRunnable::unfolded_contraction(int r, int nv, int dim) TestRunnable::unfolded_contraction(int r, int nv, int dim)
{ {
Factory fact; Factory fact;
Vector *x = fact.makeVector(nv); Vector x{fact.makeVector(nv)};
auto *forig = fact.make<FFSTensor>(r, nv, dim); auto *forig = fact.make<FFSTensor>(r, nv, dim);
UFSTensor uorig(*forig); UFSTensor uorig(*forig);
@ -482,17 +481,17 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
clock_t ctime = clock(); clock_t ctime = clock();
for (int d = dim-1; d > 0; d--) for (int d = dim-1; d > 0; d--)
{ {
UFSTensor *unew = new UFSTensor(*u, ConstVector(*x)); UFSTensor *unew = new UFSTensor(*u, ConstVector(x));
delete u; delete u;
u = unew; u = unew;
} }
ctime = clock() - ctime; ctime = clock() - ctime;
Vector res(uorig.nrows()); Vector res(uorig.nrows());
res.zeros(); res.zeros();
u->multaVec(res, *x); u->multaVec(res, x);
clock_t utime = clock(); clock_t utime = clock();
URSingleTensor ux(*x, dim); URSingleTensor ux(x, dim);
Vector v(uorig.nrows()); Vector v(uorig.nrows());
v.zeros(); v.zeros();
uorig.multaVec(v, ux.getData()); uorig.multaVec(v, ux.getData());
@ -507,7 +506,6 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
<< "\terror norm1: " << v.getNorm1() << '\n'; << "\terror norm1: " << v.getNorm1() << '\n';
delete u; delete u;
delete x;
return (v.getMax() < 1.e-10); return (v.getMax() < 1.e-10);
} }
@ -516,7 +514,7 @@ bool
TestRunnable::poly_eval(int r, int nv, int maxdim) TestRunnable::poly_eval(int r, int nv, int maxdim)
{ {
Factory fact; Factory fact;
Vector *x = fact.makeVector(nv); Vector x{fact.makeVector(nv)};
Vector out_ft(r); Vector out_ft(r);
out_ft.zeros(); out_ft.zeros();
@ -532,13 +530,13 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
FTensorPolynomial *fp = fact.makePoly<FFSTensor, FTensorPolynomial>(r, nv, maxdim); FTensorPolynomial *fp = fact.makePoly<FFSTensor, FTensorPolynomial>(r, nv, maxdim);
clock_t ft_cl = clock(); clock_t ft_cl = clock();
fp->evalTrad(out_ft, *x); fp->evalTrad(out_ft, x);
ft_cl = clock() - ft_cl; ft_cl = clock() - ft_cl;
std::cout << "\ttime for folded power eval: " std::cout << "\ttime for folded power eval: "
<< static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n'; << static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n';
clock_t fh_cl = clock(); clock_t fh_cl = clock();
fp->evalHorner(out_fh, *x); fp->evalHorner(out_fh, x);
fh_cl = clock() - fh_cl; fh_cl = clock() - fh_cl;
std::cout << "\ttime for folded horner eval: " std::cout << "\ttime for folded horner eval: "
<< static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n'; << static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n';
@ -548,13 +546,13 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
} }
clock_t ut_cl = clock(); clock_t ut_cl = clock();
up->evalTrad(out_ut, *x); up->evalTrad(out_ut, x);
ut_cl = clock() - ut_cl; ut_cl = clock() - ut_cl;
std::cout << "\ttime for unfolded power eval: " std::cout << "\ttime for unfolded power eval: "
<< static_cast<double>(ut_cl)/CLOCKS_PER_SEC << '\n'; << static_cast<double>(ut_cl)/CLOCKS_PER_SEC << '\n';
clock_t uh_cl = clock(); clock_t uh_cl = clock();
up->evalHorner(out_uh, *x); up->evalHorner(out_uh, x);
uh_cl = clock() - uh_cl; uh_cl = clock() - uh_cl;
std::cout << "\ttime for unfolded horner eval: " std::cout << "\ttime for unfolded horner eval: "
<< static_cast<double>(uh_cl)/CLOCKS_PER_SEC << '\n'; << static_cast<double>(uh_cl)/CLOCKS_PER_SEC << '\n';
@ -571,7 +569,6 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
<< "\tunfolded horner error norm max: " << max_uh << '\n'; << "\tunfolded horner error norm max: " << max_uh << '\n';
delete up; delete up;
delete x;
return (max_ft+max_fh+max_uh < 1.0e-10); return (max_ft+max_fh+max_uh < 1.0e-10);
} }