Dynare++ tensor library: first batch of simplifications for tests
Also make test failures fatal.time-shift
Sébastien Villemot
parent
44d47ee560
commit
8e52940181
|
|
@ -5,9 +5,7 @@
|
|||
#include "tl_exception.hh"
|
||||
#include "fs_tensor.hh"
|
||||
|
||||
#include <cstdlib>
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <iostream>
|
||||
|
||||
IntGenerator intgen;
|
||||
|
||||
|
|
@ -17,23 +15,23 @@ IntGenerator::init(int nf, int ny, int nv, int nw, int nu,
|
|||
{
|
||||
maxim = mx;
|
||||
probab = prob;
|
||||
long int seed = nf;
|
||||
decltype(mtgen)::result_type seed = nf;
|
||||
seed = 256*seed + ny;
|
||||
seed = 256*seed + nv;
|
||||
seed = 256*seed + nw;
|
||||
seed = 256*seed + nu;
|
||||
srand48(seed);
|
||||
mtgen.seed(seed);
|
||||
}
|
||||
|
||||
int
|
||||
IntGenerator::get() const
|
||||
IntGenerator::get()
|
||||
{
|
||||
double d = drand48();
|
||||
auto num_inter = (int) (((double) 2*maxim)/(1.0-probab));
|
||||
double d = dis(mtgen);
|
||||
auto num_inter = static_cast<int>(static_cast<double>(2*maxim)/(1.0-probab));
|
||||
int num_zero_inter = num_inter - 2*maxim;
|
||||
if (d < ((double) num_zero_inter)/num_inter)
|
||||
if (d < static_cast<double>(num_zero_inter)/num_inter)
|
||||
return 0;
|
||||
return (int) (d*num_inter)-num_zero_inter-maxim;
|
||||
return static_cast<int>((d*num_inter)-num_zero_inter-maxim);
|
||||
}
|
||||
|
||||
Monom::Monom(int len)
|
||||
|
|
@ -81,28 +79,17 @@ Monom::multiplyWith(int ex, const Monom &m)
|
|||
void
|
||||
Monom::print() const
|
||||
{
|
||||
printf("[");
|
||||
std::cout << '[';
|
||||
for (int i = 0; i < size(); i++)
|
||||
printf("%3d", operator[](i));
|
||||
printf("]");
|
||||
std::cout << operator[](i);
|
||||
std::cout << ']';
|
||||
}
|
||||
|
||||
Monom1Vector::Monom1Vector(int nxx, int l)
|
||||
: nx(nxx), len(l), x(new Monom *[len])
|
||||
: nx(nxx), len(l)
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
x[i] = new Monom(nx);
|
||||
}
|
||||
}
|
||||
|
||||
Monom1Vector::~Monom1Vector()
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
delete x[i];
|
||||
}
|
||||
delete [] x;
|
||||
x.emplace_back(nx);
|
||||
}
|
||||
|
||||
void
|
||||
|
|
@ -112,9 +99,7 @@ Monom1Vector::deriv(const IntSequence &c, Vector &out) const
|
|||
"Wrong length of output vector in Monom1Vector::deriv");
|
||||
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
out[i] = x[i]->deriv(c);
|
||||
}
|
||||
out[i] = x[i].deriv(c);
|
||||
}
|
||||
|
||||
FGSTensor *
|
||||
|
|
@ -133,60 +118,46 @@ Monom1Vector::deriv(int dim) const
|
|||
void
|
||||
Monom1Vector::print() const
|
||||
{
|
||||
printf("Variables: x(%d)\n", nx);
|
||||
printf("Rows: %d\n", len);
|
||||
std::cout << "Variables: x(" << nx << ")\n"
|
||||
<< "Rows: " << len << '\n';
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
printf("%2d: ", i);
|
||||
x[i]->print();
|
||||
printf("\n");
|
||||
std::cout << i << ": ";
|
||||
x[i].print();
|
||||
std::cout << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
Monom2Vector::Monom2Vector(int nyy, int nuu, int l)
|
||||
: ny(nyy), nu(nuu), len(l), y(new Monom *[len]), u(new Monom *[len])
|
||||
: ny(nyy), nu(nuu), len(l)
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
y[i] = new Monom(ny);
|
||||
u[i] = new Monom(nu);
|
||||
y.emplace_back(ny);
|
||||
u.emplace_back(nu);
|
||||
}
|
||||
}
|
||||
|
||||
Monom2Vector::Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon)
|
||||
: ny(xmon.ny), nu(xmon.nu), len(g.len),
|
||||
y(new Monom *[len]), u(new Monom *[len])
|
||||
: ny(xmon.ny), nu(xmon.nu), len(g.len)
|
||||
{
|
||||
TL_RAISE_IF(xmon.len != g.nx,
|
||||
"Wrong number of x's in Monom2Vector constructor");
|
||||
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
y[i] = new Monom(ny, 0);
|
||||
u[i] = new Monom(nu, 0);
|
||||
y.emplace_back(ny, 0);
|
||||
u.emplace_back(nu, 0);
|
||||
}
|
||||
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
// multiply from xmon
|
||||
for (int j = 0; j < g.nx; j++)
|
||||
{
|
||||
int ex = g.x[i]->operator[](j);
|
||||
y[i]->multiplyWith(ex, *(xmon.y[j]));
|
||||
u[i]->multiplyWith(ex, *(xmon.u[j]));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Monom2Vector::~Monom2Vector()
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
delete y[i];
|
||||
delete u[i];
|
||||
}
|
||||
delete [] y;
|
||||
delete [] u;
|
||||
// multiply from xmon
|
||||
for (int j = 0; j < g.nx; j++)
|
||||
{
|
||||
int ex = g.x[i].operator[](j);
|
||||
y[i].multiplyWith(ex, xmon.y[j]);
|
||||
u[i].multiplyWith(ex, xmon.u[j]);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
|
|
@ -202,15 +173,13 @@ Monom2Vector::deriv(const Symmetry &s, const IntSequence &c,
|
|||
IntSequence cy(c, 0, s[0]);
|
||||
IntSequence cu(c, s[0], s.dimen());
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
out[i] = y[i]->deriv(cy) * u[i]->deriv(cu);
|
||||
}
|
||||
out[i] = y[i].deriv(cy) * u[i].deriv(cu);
|
||||
}
|
||||
|
||||
FGSTensor *
|
||||
Monom2Vector::deriv(const Symmetry &s) const
|
||||
{
|
||||
IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
|
||||
IntSequence nvs{ny, nu};
|
||||
FGSTensor *t = new FGSTensor(len, TensorDimens(s, nvs));
|
||||
for (Tensor::index it = t->begin(); it != t->end(); ++it)
|
||||
{
|
||||
|
|
@ -225,96 +194,63 @@ Monom2Vector::deriv(int maxdim) const
|
|||
{
|
||||
auto *res = new FGSContainer(2);
|
||||
for (int dim = 1; dim <= maxdim; dim++)
|
||||
{
|
||||
for (int ydim = 0; ydim <= dim; ydim++)
|
||||
{
|
||||
int udim = dim - ydim;
|
||||
Symmetry s{ydim, udim};
|
||||
res->insert(deriv(s));
|
||||
}
|
||||
}
|
||||
for (int ydim = 0; ydim <= dim; ydim++)
|
||||
{
|
||||
int udim = dim - ydim;
|
||||
Symmetry s{ydim, udim};
|
||||
res->insert(deriv(s));
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
void
|
||||
Monom2Vector::print() const
|
||||
{
|
||||
printf("Variables: y(%d) u(%d)\n", ny, nu);
|
||||
printf("Rows: %d\n", len);
|
||||
std::cout << "Variables: y(" << ny << ") u(" << nu << ")\n"
|
||||
<< "Rows: " << len << '\n';
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
printf("%2d: ", i);
|
||||
y[i]->print();
|
||||
printf(" ");
|
||||
u[i]->print();
|
||||
printf("\n");
|
||||
std::cout << i;
|
||||
y[i].print();
|
||||
std::cout << " ";
|
||||
u[i].print();
|
||||
std::cout << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
Monom4Vector::~Monom4Vector()
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
delete x1[i];
|
||||
delete x2[i];
|
||||
delete x3[i];
|
||||
delete x4[i];
|
||||
}
|
||||
delete [] x1;
|
||||
delete [] x2;
|
||||
delete [] x3;
|
||||
delete [] x4;
|
||||
}
|
||||
|
||||
void
|
||||
Monom4Vector::init_random()
|
||||
{
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
x1[i] = new Monom(nx1);
|
||||
x2[i] = new Monom(nx2);
|
||||
x3[i] = new Monom(nx3);
|
||||
x4[i] = new Monom(nx4);
|
||||
x1.emplace_back(nx1);
|
||||
x2.emplace_back(nx2);
|
||||
x3.emplace_back(nx3);
|
||||
x4.emplace_back(nx4);
|
||||
}
|
||||
}
|
||||
|
||||
Monom4Vector::Monom4Vector(int l, int ny, int nu)
|
||||
: len(l), nx1(ny), nx2(nu), nx3(0), nx4(1),
|
||||
x1(new Monom *[len]),
|
||||
x2(new Monom *[len]),
|
||||
x3(new Monom *[len]),
|
||||
x4(new Monom *[len])
|
||||
: len(l), nx1(ny), nx2(nu), nx3(0), nx4(1)
|
||||
{
|
||||
init_random();
|
||||
}
|
||||
|
||||
Monom4Vector::Monom4Vector(int l, int ny, int nu, int nup)
|
||||
: len(l), nx1(ny), nx2(nu), nx3(nup), nx4(1),
|
||||
x1(new Monom *[len]),
|
||||
x2(new Monom *[len]),
|
||||
x3(new Monom *[len]),
|
||||
x4(new Monom *[len])
|
||||
: len(l), nx1(ny), nx2(nu), nx3(nup), nx4(1)
|
||||
{
|
||||
init_random();
|
||||
}
|
||||
|
||||
Monom4Vector::Monom4Vector(int l, int nbigg, int ng, int ny, int nu)
|
||||
: len(l), nx1(nbigg), nx2(ng), nx3(ny), nx4(nu),
|
||||
x1(new Monom *[len]),
|
||||
x2(new Monom *[len]),
|
||||
x3(new Monom *[len]),
|
||||
x4(new Monom *[len])
|
||||
: len(l), nx1(nbigg), nx2(ng), nx3(ny), nx4(nu)
|
||||
{
|
||||
init_random();
|
||||
}
|
||||
|
||||
Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
|
||||
const Monom4Vector &g)
|
||||
: len(f.len), nx1(bigg.nx1), nx2(bigg.nx2), nx3(bigg.nx3), nx4(1),
|
||||
x1(new Monom *[len]),
|
||||
x2(new Monom *[len]),
|
||||
x3(new Monom *[len]),
|
||||
x4(new Monom *[len])
|
||||
: len(f.len), nx1(bigg.nx1), nx2(bigg.nx2), nx3(bigg.nx3), nx4(1)
|
||||
{
|
||||
TL_RAISE_IF(!(bigg.nx1 == g.nx1 && bigg.nx2 == g.nx2 && g.nx3 == 0
|
||||
&& bigg.nx4 == 1 && g.nx4 == 1),
|
||||
|
|
@ -325,10 +261,10 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
|
|||
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
x1[i] = new Monom(nx1, 0);
|
||||
x2[i] = new Monom(nx2, 0);
|
||||
x3[i] = new Monom(nx3, 0);
|
||||
x4[i] = new Monom(nx4, 0);
|
||||
x1.emplace_back(nx1, 0);
|
||||
x2.emplace_back(nx2, 0);
|
||||
x3.emplace_back(nx3, 0);
|
||||
x4.emplace_back(nx4, 0);
|
||||
}
|
||||
|
||||
for (int i = 0; i < len; i++)
|
||||
|
|
@ -336,24 +272,24 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
|
|||
// multiply from G (first argument)
|
||||
for (int j = 0; j < f.nx1; j++)
|
||||
{
|
||||
int ex = f.x1[i]->operator[](j);
|
||||
x1[i]->multiplyWith(ex, *(bigg.x1[j]));
|
||||
x2[i]->multiplyWith(ex, *(bigg.x2[j]));
|
||||
x3[i]->multiplyWith(ex, *(bigg.x3[j]));
|
||||
x4[i]->multiplyWith(ex, *(bigg.x4[j]));
|
||||
int ex = f.x1[i].operator[](j);
|
||||
x1[i].multiplyWith(ex, bigg.x1[j]);
|
||||
x2[i].multiplyWith(ex, bigg.x2[j]);
|
||||
x3[i].multiplyWith(ex, bigg.x3[j]);
|
||||
x4[i].multiplyWith(ex, bigg.x4[j]);
|
||||
}
|
||||
// multiply from g (second argument)
|
||||
for (int j = 0; j < f.nx2; j++)
|
||||
{
|
||||
int ex = f.x2[i]->operator[](j);
|
||||
x1[i]->multiplyWith(ex, *(g.x1[j]));
|
||||
x2[i]->multiplyWith(ex, *(g.x2[j]));
|
||||
x4[i]->multiplyWith(ex, *(g.x4[j]));
|
||||
int ex = f.x2[i].operator[](j);
|
||||
x1[i].multiplyWith(ex, g.x1[j]);
|
||||
x2[i].multiplyWith(ex, g.x2[j]);
|
||||
x4[i].multiplyWith(ex, g.x4[j]);
|
||||
}
|
||||
// add y as third argument of f
|
||||
x1[i]->add(1, *(f.x3[i]));
|
||||
x1[i].add(1, f.x3[i]);
|
||||
// add u as fourth argument of f
|
||||
x2[i]->add(1, *(f.x4[i]));
|
||||
x2[i].add(1, f.x4[i]);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -372,23 +308,20 @@ Monom4Vector::deriv(const Symmetry &s, const IntSequence &coor,
|
|||
{
|
||||
out[i] = 1;
|
||||
int off = 0;
|
||||
out[i] *= x1[i]->deriv(IntSequence(coor, off, off+s[0]));
|
||||
out[i] *= x1[i].deriv(IntSequence(coor, off, off+s[0]));
|
||||
off += s[0];
|
||||
out[i] *= x2[i]->deriv(IntSequence(coor, off, off+s[1]));
|
||||
out[i] *= x2[i].deriv(IntSequence(coor, off, off+s[1]));
|
||||
off += s[1];
|
||||
out[i] *= x3[i]->deriv(IntSequence(coor, off, off+s[2]));
|
||||
out[i] *= x3[i].deriv(IntSequence(coor, off, off+s[2]));
|
||||
off += s[2];
|
||||
out[i] *= x4[i]->deriv(IntSequence(coor, off, off+s[3]));
|
||||
out[i] *= x4[i].deriv(IntSequence(coor, off, off+s[3]));
|
||||
}
|
||||
}
|
||||
|
||||
FGSTensor *
|
||||
Monom4Vector::deriv(const Symmetry &s) const
|
||||
{
|
||||
IntSequence nvs(4);
|
||||
nvs[0] = nx1; nvs[1] = nx2;
|
||||
nvs[2] = nx3; nvs[3] = nx4;
|
||||
|
||||
IntSequence nvs{nx1, nx2, nx3, nx4};
|
||||
FGSTensor *res = new FGSTensor(len, TensorDimens(s, nvs));
|
||||
for (Tensor::index run = res->begin(); run != res->end(); ++run)
|
||||
{
|
||||
|
|
@ -401,8 +334,7 @@ Monom4Vector::deriv(const Symmetry &s) const
|
|||
FSSparseTensor *
|
||||
Monom4Vector::deriv(int dim) const
|
||||
{
|
||||
IntSequence cum(4);
|
||||
cum[0] = 0; cum[1] = nx1; cum[2] = nx1+nx2; cum[3] = nx1+nx2+nx3;
|
||||
IntSequence cum{0, nx1, nx1+nx2, nx1+nx2+nx3};
|
||||
|
||||
auto *res = new FSSparseTensor(dim, nx1+nx2+nx3+nx4, len);
|
||||
|
||||
|
|
@ -423,12 +355,8 @@ Monom4Vector::deriv(int dim) const
|
|||
Vector col(len);
|
||||
deriv(ind_sym, ind, col);
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
if (col[i] != 0.0)
|
||||
{
|
||||
res->insert(run.getCoor(), i, col[i]);
|
||||
}
|
||||
}
|
||||
if (col[i] != 0.0)
|
||||
res->insert(run.getCoor(), i, col[i]);
|
||||
}
|
||||
|
||||
return res;
|
||||
|
|
@ -437,27 +365,26 @@ Monom4Vector::deriv(int dim) const
|
|||
void
|
||||
Monom4Vector::print() const
|
||||
{
|
||||
printf("Variables: x1(%d) x2(%d) x3(%d) x4(%d)\n",
|
||||
nx1, nx2, nx3, nx4);
|
||||
printf("Rows: %d\n", len);
|
||||
std::cout << "Variables: x1(" << nx1 << ") x2(" << nx2
|
||||
<< ") x3(" << nx3 << ") x4(" << nx4 << ")\n"
|
||||
<< "Rows: " << len << '\n';
|
||||
for (int i = 0; i < len; i++)
|
||||
{
|
||||
printf("%2d: ", i);
|
||||
x1[i]->print();
|
||||
printf(" ");
|
||||
x2[i]->print();
|
||||
printf(" ");
|
||||
x3[i]->print();
|
||||
printf(" ");
|
||||
x4[i]->print();
|
||||
printf("\n");
|
||||
std::cout << i << ": ";
|
||||
x1[i].print();
|
||||
std::cout << " ";
|
||||
x2[i].print();
|
||||
std::cout << " ";
|
||||
x3[i].print();
|
||||
std::cout << " ";
|
||||
x4[i].print();
|
||||
std::cout << '\n';
|
||||
}
|
||||
}
|
||||
|
||||
SparseDerivGenerator::SparseDerivGenerator(
|
||||
int nf, int ny, int nu, int nup, int nbigg, int ng,
|
||||
SparseDerivGenerator::SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
|
||||
int mx, double prob, int maxdim)
|
||||
: maxdimen(maxdim), ts(new FSSparseTensor *[maxdimen])
|
||||
: maxdimen(maxdim), bigg(4), g(4), rcont(4), ts(new FSSparseTensor *[maxdimen])
|
||||
{
|
||||
intgen.init(nf, ny, nu, nup, nbigg, mx, prob);
|
||||
|
||||
|
|
@ -465,18 +392,15 @@ SparseDerivGenerator::SparseDerivGenerator(
|
|||
Monom4Vector g_m(ng, ny, nu);
|
||||
Monom4Vector f(nf, nbigg, ng, ny, nu);
|
||||
Monom4Vector r(f, bigg_m, g_m);
|
||||
bigg = new FGSContainer(4);
|
||||
g = new FGSContainer(4);
|
||||
rcont = new FGSContainer(4);
|
||||
|
||||
for (int dim = 1; dim <= maxdimen; dim++)
|
||||
{
|
||||
for (auto &si : SymmetrySet(dim, 4))
|
||||
{
|
||||
bigg->insert(bigg_m.deriv(si));
|
||||
rcont->insert(r.deriv(si));
|
||||
bigg.insert(bigg_m.deriv(si));
|
||||
rcont.insert(r.deriv(si));
|
||||
if (si[2] == 0)
|
||||
g->insert(g_m.deriv(si));
|
||||
g.insert(g_m.deriv(si));
|
||||
}
|
||||
|
||||
ts[dim-1] = f.deriv(dim);
|
||||
|
|
@ -485,9 +409,6 @@ SparseDerivGenerator::SparseDerivGenerator(
|
|||
|
||||
SparseDerivGenerator::~SparseDerivGenerator()
|
||||
{
|
||||
delete bigg;
|
||||
delete g;
|
||||
delete rcont;
|
||||
for (int i = 0; i < maxdimen; i++)
|
||||
delete ts[i];
|
||||
delete [] ts;
|
||||
|
|
@ -517,9 +438,7 @@ DenseDerivGenerator::unfold()
|
|||
{
|
||||
uxcont = new UGSContainer(*xcont);
|
||||
for (int i = 0; i < maxdimen; i++)
|
||||
{
|
||||
uts[i] = new UGSTensor(*(ts[i]));
|
||||
}
|
||||
uts[i] = new UGSTensor(*(ts[i]));
|
||||
}
|
||||
|
||||
DenseDerivGenerator::~DenseDerivGenerator()
|
||||
|
|
|
|||
|
|
@ -4,6 +4,9 @@
|
|||
#ifndef MONOMS_H
|
||||
#define MONOMS_H
|
||||
|
||||
#include <vector>
|
||||
#include <random>
|
||||
|
||||
#include "int_sequence.hh"
|
||||
#include "gs_tensor.hh"
|
||||
#include "t_container.hh"
|
||||
|
|
@ -14,12 +17,12 @@ class IntGenerator
|
|||
{
|
||||
int maxim{5};
|
||||
double probab{0.3};
|
||||
std::mt19937 mtgen;
|
||||
std::uniform_real_distribution<> dis;
|
||||
public:
|
||||
IntGenerator()
|
||||
|
||||
= default;
|
||||
IntGenerator() = default;
|
||||
void init(int nf, int ny, int nv, int nw, int nu, int mx, double prob);
|
||||
int get() const;
|
||||
int get();
|
||||
};
|
||||
|
||||
extern IntGenerator intgen;
|
||||
|
|
@ -41,10 +44,10 @@ class Monom1Vector
|
|||
friend class Monom2Vector;
|
||||
int nx;
|
||||
int len;
|
||||
Monom **const x;
|
||||
std::vector<Monom> x;
|
||||
public:
|
||||
Monom1Vector(int nxx, int l);
|
||||
~Monom1Vector();
|
||||
~Monom1Vector() = default;
|
||||
void deriv(const IntSequence &c, Vector &out) const;
|
||||
FGSTensor *deriv(int dim) const;
|
||||
void print() const;
|
||||
|
|
@ -53,17 +56,15 @@ public:
|
|||
//class Monom3Vector;
|
||||
class Monom2Vector
|
||||
{
|
||||
int ny;
|
||||
int nu;
|
||||
int ny, nu;
|
||||
int len;
|
||||
Monom **const y;
|
||||
Monom **const u;
|
||||
std::vector<Monom> y, u;
|
||||
public:
|
||||
// generate random vector of monom two vector
|
||||
Monom2Vector(int nyy, int nuu, int l);
|
||||
// calculate g(x(y,u))
|
||||
Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon);
|
||||
~Monom2Vector();
|
||||
~Monom2Vector() = default;
|
||||
void deriv(const Symmetry &s, const IntSequence &c, Vector &out) const;
|
||||
FGSTensor *deriv(const Symmetry &s) const;
|
||||
FGSContainer *deriv(int maxdim) const;
|
||||
|
|
@ -73,14 +74,8 @@ public:
|
|||
class Monom4Vector
|
||||
{
|
||||
int len;
|
||||
int nx1;
|
||||
int nx2;
|
||||
int nx3;
|
||||
int nx4;
|
||||
Monom **const x1;
|
||||
Monom **const x2;
|
||||
Monom **const x3;
|
||||
Monom **const x4;
|
||||
int nx1, nx2, nx3, nx4;
|
||||
std::vector<Monom> x1, x2, x3, x4;
|
||||
public:
|
||||
/* random for g(y,u,sigma) */
|
||||
Monom4Vector(int l, int ny, int nu);
|
||||
|
|
@ -91,7 +86,7 @@ public:
|
|||
/* substitution f(G(y,u,u',sigma),g(y,u,sigma),y,u) */
|
||||
Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
|
||||
const Monom4Vector &g);
|
||||
~Monom4Vector();
|
||||
~Monom4Vector() = default;
|
||||
FSSparseTensor *deriv(int dim) const;
|
||||
FGSTensor *deriv(const Symmetry &s) const;
|
||||
void deriv(const Symmetry &s, const IntSequence &coor, Vector &out) const;
|
||||
|
|
@ -103,9 +98,9 @@ protected:
|
|||
struct SparseDerivGenerator
|
||||
{
|
||||
int maxdimen;
|
||||
FGSContainer *bigg;
|
||||
FGSContainer *g;
|
||||
FGSContainer *rcont;
|
||||
FGSContainer bigg;
|
||||
FGSContainer g;
|
||||
FGSContainer rcont;
|
||||
FSSparseTensor **const ts;
|
||||
SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
|
||||
int mx, double prob, int maxdim);
|
||||
|
|
|
|||
|
|
@ -13,28 +13,28 @@
|
|||
#include "ps_tensor.hh"
|
||||
#include "tl_static.hh"
|
||||
|
||||
#include <cstdio>
|
||||
#include <cstring>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <memory>
|
||||
#include <vector>
|
||||
#include <ctime>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
|
||||
class TestRunnable
|
||||
{
|
||||
char name[100];
|
||||
public:
|
||||
const std::string name;
|
||||
int dim; // dimension of the solved problem
|
||||
int nvar; // number of variable of the solved problem
|
||||
TestRunnable(const char *n, int d, int nv)
|
||||
: dim(d), nvar(nv)
|
||||
TestRunnable(std::string name_arg, int d, int nv)
|
||||
: name{std::move(name_arg)}, dim(d), nvar(nv)
|
||||
{
|
||||
strncpy(name, n, 100);
|
||||
}
|
||||
virtual ~TestRunnable() = default;
|
||||
bool test() const;
|
||||
virtual bool run() const = 0;
|
||||
const char *
|
||||
getName() const
|
||||
{
|
||||
return name;
|
||||
}
|
||||
protected:
|
||||
template<class _Ttype>
|
||||
static bool index_forward(const Symmetry &s, const IntSequence &nvs);
|
||||
|
|
@ -92,22 +92,17 @@ protected:
|
|||
bool
|
||||
TestRunnable::test() const
|
||||
{
|
||||
printf("Running test <%s>\n", name);
|
||||
std::cout << "Running test <" << name << ">" << std::endl;
|
||||
clock_t start = clock();
|
||||
bool passed = run();
|
||||
clock_t end = clock();
|
||||
printf("CPU time %8.4g (CPU seconds)..................",
|
||||
((double) (end-start))/CLOCKS_PER_SEC);
|
||||
std::cout << "CPU time " << static_cast<double>(end-start)/CLOCKS_PER_SEC
|
||||
<< " (CPU seconds)..................";
|
||||
if (passed)
|
||||
{
|
||||
printf("passed\n\n");
|
||||
return passed;
|
||||
}
|
||||
std::cout << "passed\n\n";
|
||||
else
|
||||
{
|
||||
printf("FAILED\n\n");
|
||||
return passed;
|
||||
}
|
||||
std::cout << "FAILED\n\n";
|
||||
return passed;
|
||||
}
|
||||
|
||||
/****************************************************/
|
||||
|
|
@ -137,10 +132,10 @@ TestRunnable::index_forward(const Symmetry &s, const IntSequence &nvs)
|
|||
}
|
||||
while (run != dummy.begin());
|
||||
|
||||
printf("\tnumber of columns = %d\n", dummy.ncols());
|
||||
printf("\tnumber of increments = %d\n", nincr);
|
||||
printf("\tnumber of decrements = %d\n", ndecr);
|
||||
printf("\tnumber of failures = %d\n", fails);
|
||||
std::cout << "\tnumber of columns = " << dummy.ncols() << '\n'
|
||||
<< "\tnumber of increments = " << nincr << '\n'
|
||||
<< "\tnumber of decrements = " << ndecr << '\n'
|
||||
<< "\tnumber of failures = " << fails << '\n';
|
||||
|
||||
return fails == 0;
|
||||
}
|
||||
|
|
@ -168,10 +163,10 @@ TestRunnable::index_backward(const Symmetry &s, const IntSequence &nvs)
|
|||
nincr++;
|
||||
}
|
||||
|
||||
printf("\tnumber of columns = %d\n", dummy.ncols());
|
||||
printf("\tnumber of increments = %d\n", nincr);
|
||||
printf("\tnumber of decrements = %d\n", ndecr);
|
||||
printf("\tnumber of failures = %d\n", fails);
|
||||
std::cout << "\tnumber of columns = " << dummy.ncols() << '\n'
|
||||
<< "\tnumber of increments = " << nincr << '\n'
|
||||
<< "\tnumber of decrements = " << ndecr << '\n'
|
||||
<< "\tnumber of failures = " << fails << '\n';
|
||||
|
||||
return fails == 0;
|
||||
}
|
||||
|
|
@ -191,9 +186,9 @@ TestRunnable::index_offset(const Symmetry &s, const IntSequence &nvs)
|
|||
fails++;
|
||||
}
|
||||
|
||||
printf("\tnumber of columns = %d\n", dummy.ncols());
|
||||
printf("\tnumber of increments = %d\n", nincr);
|
||||
printf("\tnumber of failures = %d\n", fails);
|
||||
std::cout << "\tnumber of columns = " << dummy.ncols() << '\n'
|
||||
<< "\tnumber of increments = " << nincr << '\n'
|
||||
<< "\tnumber of failures = " << fails << '\n';
|
||||
|
||||
return fails == 0;
|
||||
}
|
||||
|
|
@ -206,10 +201,10 @@ TestRunnable::fold_unfold(const FTensor *folded)
|
|||
folded2->add(-1.0, *folded);
|
||||
double normInf = folded2->getNormInf();
|
||||
double norm1 = folded2->getNorm1();
|
||||
printf("\tfolded size: (%d, %d)\n", folded->nrows(), folded->ncols());
|
||||
printf("\tunfolded size: (%d, %d)\n", unfolded->nrows(), unfolded->ncols());
|
||||
printf("\tdifference normInf: %8.4g\n", normInf);
|
||||
printf("\tdifference norm1: %8.4g\n", norm1);
|
||||
std::cout << "\tfolded size: (" << folded->nrows() << ", " << folded->ncols() << ")\n"
|
||||
<< "\tunfolded size: (" << unfolded->nrows() << ", " << unfolded->ncols() << ")\n"
|
||||
<< "\tdifference normInf: " << normInf << '\n'
|
||||
<< "\tdifference norm1: " << norm1 << '\n';
|
||||
|
||||
delete folded;
|
||||
delete unfolded;
|
||||
|
|
@ -247,15 +242,12 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
|
|||
double norm1 = btmp.getNorm1();
|
||||
double normInf = btmp.getNormInf();
|
||||
|
||||
printf("\ttime for folded product: %8.4g\n",
|
||||
((double) (s2-s1))/CLOCKS_PER_SEC);
|
||||
printf("\ttime for unfolded product: %8.4g\n",
|
||||
((double) (s5-s4))/CLOCKS_PER_SEC);
|
||||
printf("\ttime for container convert: %8.4g\n",
|
||||
((double) (s3-s2))/CLOCKS_PER_SEC);
|
||||
printf("\tunfolded difference normMax: %10.6g\n", norm);
|
||||
printf("\tunfolded difference norm1: %10.6g\n", norm1);
|
||||
printf("\tunfolded difference normInf: %10.6g\n", normInf);
|
||||
std::cout << "\ttime for folded product: " << static_cast<double>(s2-s1)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\ttime for unfolded product: " << static_cast<double>(s5-s4)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\ttime for container convert: " << static_cast<double>(s3-s2)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\tunfolded difference normMax: " << norm << '\n'
|
||||
<< "\tunfolded difference norm1: " << norm1 << '\n'
|
||||
<< "\tunfolded difference normInf: " << normInf << '\n';
|
||||
|
||||
delete cont;
|
||||
delete fh;
|
||||
|
|
@ -269,28 +261,27 @@ TestRunnable::folded_monomial(int ng, int nx, int ny, int nu, int dim)
|
|||
clock_t gen_time = clock();
|
||||
DenseDerivGenerator gen(ng, nx, ny, nu, 5, 0.3, dim);
|
||||
gen_time = clock()-gen_time;
|
||||
printf("\ttime for monom generation: %8.4g\n",
|
||||
((double) gen_time)/CLOCKS_PER_SEC);
|
||||
IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
|
||||
std::cout << "\ttime for monom generation: "
|
||||
<< static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
|
||||
IntSequence nvs{ny, nu};
|
||||
double maxnorm = 0;
|
||||
for (int ydim = 0; ydim <= dim; ydim++)
|
||||
{
|
||||
Symmetry s{ydim, dim-ydim};
|
||||
printf("\tSymmetry: "); s.print();
|
||||
std::cout << "\tSymmetry: ";
|
||||
s.print();
|
||||
FGSTensor res(ng, TensorDimens(s, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int d = 1; d <= dim; d++)
|
||||
{
|
||||
gen.xcont->multAndAdd(*(gen.ts[d-1]), res);
|
||||
}
|
||||
gen.xcont->multAndAdd(*(gen.ts[d-1]), res);
|
||||
stime = clock() - stime;
|
||||
printf("\t\ttime for symmetry: %8.4g\n",
|
||||
((double) stime)/CLOCKS_PER_SEC);
|
||||
std::cout << "\t\ttime for symmetry: "
|
||||
<< static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
|
||||
const FGSTensor *mres = gen.rcont->get(s);
|
||||
res.add(-1.0, *mres);
|
||||
double normtmp = res.getData().getMax();
|
||||
printf("\t\terror normMax: %10.6g\n", normtmp);
|
||||
std::cout << "\t\terror normMax: " << normtmp << '\n';
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
|
|
@ -303,34 +294,33 @@ TestRunnable::unfolded_monomial(int ng, int nx, int ny, int nu, int dim)
|
|||
clock_t gen_time = clock();
|
||||
DenseDerivGenerator gen(ng, nx, ny, nu, 5, 0.3, dim);
|
||||
gen_time = clock()-gen_time;
|
||||
printf("\ttime for monom generation: %8.4g\n",
|
||||
((double) gen_time)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for monom generation: "
|
||||
<< static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
|
||||
clock_t u_time = clock();
|
||||
gen.unfold();
|
||||
u_time = clock() - u_time;
|
||||
printf("\ttime for monom unfolding: %8.4g\n",
|
||||
((double) u_time)/CLOCKS_PER_SEC);
|
||||
IntSequence nvs(2); nvs[0] = ny; nvs[1] = nu;
|
||||
std::cout << "\ttime for monom unfolding: "
|
||||
<< static_cast<double>(u_time)/CLOCKS_PER_SEC << '\n';
|
||||
IntSequence nvs{ny, nu};
|
||||
double maxnorm = 0;
|
||||
for (int ydim = 0; ydim <= dim; ydim++)
|
||||
{
|
||||
Symmetry s{ydim, dim-ydim};
|
||||
printf("\tSymmetry: "); s.print();
|
||||
std::cout << "\tSymmetry: ";
|
||||
s.print();
|
||||
UGSTensor res(ng, TensorDimens(s, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int d = 1; d <= dim; d++)
|
||||
{
|
||||
gen.uxcont->multAndAdd(*(gen.uts[d-1]), res);
|
||||
}
|
||||
gen.uxcont->multAndAdd(*(gen.uts[d-1]), res);
|
||||
stime = clock() - stime;
|
||||
printf("\t\ttime for symmetry: %8.4g\n",
|
||||
((double) stime)/CLOCKS_PER_SEC);
|
||||
std::cout << "\t\ttime for symmetry: "
|
||||
<< static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
|
||||
const FGSTensor *mres = gen.rcont->get(s);
|
||||
FGSTensor foldres(res);
|
||||
foldres.add(-1.0, *mres);
|
||||
double normtmp = foldres.getData().getMax();
|
||||
printf("\t\terror normMax: %10.6g\n", normtmp);
|
||||
std::cout << "\t\terror normMax: " << normtmp << '\n';
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
|
|
@ -346,42 +336,40 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg,
|
|||
5, 0.55, dim);
|
||||
gen_time = clock()-gen_time;
|
||||
for (int d = 1; d <= dim; d++)
|
||||
{
|
||||
printf("\tfill of dim=%d tensor: %3.2f %%\n",
|
||||
d, 100*dg.ts[d-1]->getFillFactor());
|
||||
}
|
||||
printf("\ttime for monom generation: %8.4g\n",
|
||||
((double) gen_time)/CLOCKS_PER_SEC);
|
||||
std::cout << "\tfill of dim=" << d << " tensor: "
|
||||
<< std::setprecision(2) << std::fixed << std::setw(6)
|
||||
<< 100*dg.ts[d-1]->getFillFactor()
|
||||
<< std::setprecision(6) << std::defaultfloat << " %\n";
|
||||
std::cout << "\ttime for monom generation: "
|
||||
<< static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
IntSequence nvs(4);
|
||||
nvs[0] = ny; nvs[1] = nu; nvs[2] = nup; nvs[3] = 1;
|
||||
IntSequence nvs{ny, nu, nup, 1};
|
||||
double maxnorm = 0.0;
|
||||
|
||||
// form ZContainer
|
||||
FoldedZContainer zc(dg.bigg, nbigg, dg.g, ng, ny, nu);
|
||||
FoldedZContainer zc(&dg.bigg, nbigg, &dg.g, ng, ny, nu);
|
||||
|
||||
for (int d = 2; d <= dim; d++)
|
||||
{
|
||||
for (auto &si : SymmetrySet(d, 4))
|
||||
{
|
||||
printf("\tSymmetry: ");
|
||||
si.print();
|
||||
FGSTensor res(nf, TensorDimens(si, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int l = 1; l <= si.dimen(); l++)
|
||||
zc.multAndAdd(*(dg.ts[l-1]), res);
|
||||
stime = clock() - stime;
|
||||
printf("\t\ttime for symmetry: %8.4g\n",
|
||||
((double) stime)/CLOCKS_PER_SEC);
|
||||
const FGSTensor *mres = dg.rcont->get(si);
|
||||
res.add(-1.0, *mres);
|
||||
double normtmp = res.getData().getMax();
|
||||
printf("\t\terror normMax: %10.6g\n", normtmp);
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
}
|
||||
for (auto &si : SymmetrySet(d, 4))
|
||||
{
|
||||
std::cout << "\tSymmetry: ";
|
||||
si.print();
|
||||
FGSTensor res(nf, TensorDimens(si, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int l = 1; l <= si.dimen(); l++)
|
||||
zc.multAndAdd(*(dg.ts[l-1]), res);
|
||||
stime = clock() - stime;
|
||||
std::cout << "\t\ttime for symmetry: "
|
||||
<< static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
|
||||
const FGSTensor *mres = dg.rcont.get(si);
|
||||
res.add(-1.0, *mres);
|
||||
double normtmp = res.getData().getMax();
|
||||
std::cout << "\t\terror normMax: " << normtmp << '\n';
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
|
||||
return maxnorm < 1.0e-10;
|
||||
}
|
||||
|
||||
|
|
@ -394,50 +382,48 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
|
|||
5, 0.55, dim);
|
||||
gen_time = clock()-gen_time;
|
||||
for (int d = 1; d <= dim; d++)
|
||||
{
|
||||
printf("\tfill of dim=%d tensor: %3.2f %%\n",
|
||||
d, 100*dg.ts[d-1]->getFillFactor());
|
||||
}
|
||||
printf("\ttime for monom generation: %8.4g\n",
|
||||
((double) gen_time)/CLOCKS_PER_SEC);
|
||||
std::cout << "\tfill of dim=" << d << " tensor: "
|
||||
<< std::setprecision(2) << std::fixed << std::setw(6)
|
||||
<< 100*dg.ts[d-1]->getFillFactor()
|
||||
<< std::setprecision(6) << std::defaultfloat << " %\n";
|
||||
std::cout << "\ttime for monom generation: "
|
||||
<< static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
clock_t con_time = clock();
|
||||
UGSContainer uG_cont(*(dg.bigg));
|
||||
UGSContainer ug_cont(*(dg.g));
|
||||
UGSContainer uG_cont(dg.bigg);
|
||||
UGSContainer ug_cont(dg.g);
|
||||
con_time = clock()-con_time;
|
||||
printf("\ttime for container unfold: %8.4g\n",
|
||||
((double) con_time)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for container unfold: "
|
||||
<< static_cast<double>(con_time)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
IntSequence nvs(4);
|
||||
nvs[0] = ny; nvs[1] = nu; nvs[2] = nup; nvs[3] = 1;
|
||||
IntSequence nvs{ny, nu, nup, 1};
|
||||
double maxnorm = 0.0;
|
||||
|
||||
// form ZContainer
|
||||
UnfoldedZContainer zc(&uG_cont, nbigg, &ug_cont, ng, ny, nu);
|
||||
|
||||
for (int d = 2; d <= dim; d++)
|
||||
{
|
||||
for (auto &si : SymmetrySet(d, 4))
|
||||
{
|
||||
printf("\tSymmetry: ");
|
||||
si.print();
|
||||
UGSTensor res(nf, TensorDimens(si, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int l = 1; l <= si.dimen(); l++)
|
||||
zc.multAndAdd(*(dg.ts[l-1]), res);
|
||||
stime = clock() - stime;
|
||||
printf("\t\ttime for symmetry: %8.4g\n",
|
||||
((double) stime)/CLOCKS_PER_SEC);
|
||||
FGSTensor fold_res(res);
|
||||
const FGSTensor *mres = dg.rcont->get(si);
|
||||
fold_res.add(-1.0, *mres);
|
||||
double normtmp = fold_res.getData().getMax();
|
||||
printf("\t\terror normMax: %10.6g\n", normtmp);
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
}
|
||||
for (auto &si : SymmetrySet(d, 4))
|
||||
{
|
||||
std::cout << "\tSymmetry: ";
|
||||
si.print();
|
||||
UGSTensor res(nf, TensorDimens(si, nvs));
|
||||
res.getData().zeros();
|
||||
clock_t stime = clock();
|
||||
for (int l = 1; l <= si.dimen(); l++)
|
||||
zc.multAndAdd(*(dg.ts[l-1]), res);
|
||||
stime = clock() - stime;
|
||||
std::cout << "\t\ttime for symmetry: "
|
||||
<< static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
|
||||
FGSTensor fold_res(res);
|
||||
const FGSTensor *mres = dg.rcont.get(si);
|
||||
fold_res.add(-1.0, *mres);
|
||||
double normtmp = fold_res.getData().getMax();
|
||||
std::cout << "\t\terror normMax: " << normtmp << '\n';
|
||||
if (normtmp > maxnorm)
|
||||
maxnorm = normtmp;
|
||||
}
|
||||
|
||||
return maxnorm < 1.0e-10;
|
||||
}
|
||||
|
||||
|
|
@ -470,12 +456,12 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
|
|||
utime = clock() - utime;
|
||||
|
||||
v.add(-1.0, res);
|
||||
printf("\ttime for folded contraction: %8.4g\n",
|
||||
((double) ctime)/CLOCKS_PER_SEC);
|
||||
printf("\ttime for unfolded power: %8.4g\n",
|
||||
((double) utime)/CLOCKS_PER_SEC);
|
||||
printf("\terror normMax: %10.6g\n", v.getMax());
|
||||
printf("\terror norm1: %10.6g\n", v.getNorm1());
|
||||
std::cout << "\ttime for folded contraction: "
|
||||
<< static_cast<double>(ctime)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\ttime for unfolded power: "
|
||||
<< static_cast<double>(utime)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\terror normMax: " << v.getMax() << '\n'
|
||||
<< "\terror norm1: " << v.getNorm1() << '\n';
|
||||
|
||||
delete f;
|
||||
delete x;
|
||||
|
|
@ -513,12 +499,12 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
|
|||
utime = clock() - utime;
|
||||
|
||||
v.add(-1.0, res);
|
||||
printf("\ttime for unfolded contraction: %8.4g\n",
|
||||
((double) ctime)/CLOCKS_PER_SEC);
|
||||
printf("\ttime for unfolded power: %8.4g\n",
|
||||
((double) utime)/CLOCKS_PER_SEC);
|
||||
printf("\terror normMax: %10.6g\n", v.getMax());
|
||||
printf("\terror norm1: %10.6g\n", v.getNorm1());
|
||||
std::cout << "\ttime for unfolded contraction: "
|
||||
<< static_cast<double>(ctime)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\ttime for unfolded power: "
|
||||
<< static_cast<double>(utime)/CLOCKS_PER_SEC << '\n'
|
||||
<< "\terror normMax: " << v.getMax() << '\n'
|
||||
<< "\terror norm1: " << v.getNorm1() << '\n';
|
||||
|
||||
delete u;
|
||||
delete x;
|
||||
|
|
@ -532,10 +518,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
|
|||
Factory fact;
|
||||
Vector *x = fact.makeVector(nv);
|
||||
|
||||
Vector out_ft(r); out_ft.zeros();
|
||||
Vector out_fh(r); out_fh.zeros();
|
||||
Vector out_ut(r); out_ut.zeros();
|
||||
Vector out_uh(r); out_uh.zeros();
|
||||
Vector out_ft(r);
|
||||
out_ft.zeros();
|
||||
Vector out_fh(r);
|
||||
out_fh.zeros();
|
||||
Vector out_ut(r);
|
||||
out_ut.zeros();
|
||||
Vector out_uh(r);
|
||||
out_uh.zeros();
|
||||
|
||||
UTensorPolynomial *up;
|
||||
{
|
||||
|
|
@ -544,14 +534,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
|
|||
clock_t ft_cl = clock();
|
||||
fp->evalTrad(out_ft, *x);
|
||||
ft_cl = clock() - ft_cl;
|
||||
printf("\ttime for folded power eval: %8.4g\n",
|
||||
((double) ft_cl)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for folded power eval: "
|
||||
<< static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
clock_t fh_cl = clock();
|
||||
fp->evalHorner(out_fh, *x);
|
||||
fh_cl = clock() - fh_cl;
|
||||
printf("\ttime for folded horner eval: %8.4g\n",
|
||||
((double) fh_cl)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for folded horner eval: "
|
||||
<< static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
up = new UTensorPolynomial(*fp);
|
||||
delete fp;
|
||||
|
|
@ -560,14 +550,14 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
|
|||
clock_t ut_cl = clock();
|
||||
up->evalTrad(out_ut, *x);
|
||||
ut_cl = clock() - ut_cl;
|
||||
printf("\ttime for unfolded power eval: %8.4g\n",
|
||||
((double) ut_cl)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for unfolded power eval: "
|
||||
<< static_cast<double>(ut_cl)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
clock_t uh_cl = clock();
|
||||
up->evalHorner(out_uh, *x);
|
||||
uh_cl = clock() - uh_cl;
|
||||
printf("\ttime for unfolded horner eval: %8.4g\n",
|
||||
((double) uh_cl)/CLOCKS_PER_SEC);
|
||||
std::cout << "\ttime for unfolded horner eval: "
|
||||
<< static_cast<double>(uh_cl)/CLOCKS_PER_SEC << '\n';
|
||||
|
||||
out_ft.add(-1.0, out_ut);
|
||||
double max_ft = out_ft.getMax();
|
||||
|
|
@ -576,9 +566,9 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
|
|||
out_uh.add(-1.0, out_ut);
|
||||
double max_uh = out_uh.getMax();
|
||||
|
||||
printf("\tfolded power error norm max: %10.6g\n", max_ft);
|
||||
printf("\tfolded horner error norm max: %10.6g\n", max_fh);
|
||||
printf("\tunfolded horner error norm max: %10.6g\n", max_uh);
|
||||
std::cout << "\tfolded power error norm max: " << max_ft << '\n'
|
||||
<< "\tfolded horner error norm max: " << max_fh << '\n'
|
||||
<< "\tunfolded horner error norm max: " << max_uh << '\n';
|
||||
|
||||
delete up;
|
||||
delete x;
|
||||
|
|
@ -599,7 +589,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3};
|
||||
IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
|
||||
IntSequence nvs{4, 2};
|
||||
return index_forward<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -615,7 +605,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3};
|
||||
IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
|
||||
IntSequence nvs{4, 2};
|
||||
return index_forward<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -631,7 +621,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
|
||||
IntSequence nvs{5, 2, 2};
|
||||
return index_forward<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -647,7 +637,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
|
||||
IntSequence nvs{5, 2, 2};
|
||||
return index_forward<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -663,7 +653,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{1, 1, 3};
|
||||
IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
|
||||
IntSequence nvs{3, 3, 2};
|
||||
return index_backward<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -679,7 +669,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
|
||||
IntSequence nvs{4, 2, 4};
|
||||
return index_backward<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -695,7 +685,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{1, 1, 3};
|
||||
IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
|
||||
IntSequence nvs{3, 3, 2};
|
||||
return index_backward<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -711,7 +701,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 4; nvs[1] = 2; nvs[2] = 4;
|
||||
IntSequence nvs{4, 2, 4};
|
||||
return index_backward<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -727,7 +717,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3};
|
||||
IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
|
||||
IntSequence nvs{4, 2};
|
||||
return index_offset<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -743,7 +733,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3};
|
||||
IntSequence nvs(2); nvs[0] = 4; nvs[1] = 2;
|
||||
IntSequence nvs{4, 2};
|
||||
return index_offset<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -759,7 +749,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
|
||||
IntSequence nvs{5, 2, 2};
|
||||
return index_offset<FGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -775,7 +765,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 3, 2};
|
||||
IntSequence nvs(3); nvs[0] = 5; nvs[1] = 2; nvs[2] = 2;
|
||||
IntSequence nvs{5, 2, 2};
|
||||
return index_offset<UGSTensor>(s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -805,7 +795,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{1, 2, 2};
|
||||
IntSequence nvs(3); nvs[0] = 3; nvs[1] = 3; nvs[2] = 2;
|
||||
IntSequence nvs{3, 3, 2};
|
||||
return gs_fold_unfold(5, s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -835,7 +825,7 @@ public:
|
|||
run() const override
|
||||
{
|
||||
Symmetry s{2, 1, 2};
|
||||
IntSequence nvs(3); nvs[0] = 6; nvs[1] = 2; nvs[2] = 6;
|
||||
IntSequence nvs{6, 2, 6};
|
||||
return gs_fold_unfold(5, s, nvs);
|
||||
}
|
||||
};
|
||||
|
|
@ -878,7 +868,7 @@ public:
|
|||
bool
|
||||
run() const override
|
||||
{
|
||||
IntSequence bnvs(2); bnvs[0] = 3; bnvs[1] = 2;
|
||||
IntSequence bnvs{3, 2};
|
||||
return dense_prod(Symmetry{1, 2}, bnvs, 2, 3, 2);
|
||||
}
|
||||
};
|
||||
|
|
@ -893,7 +883,7 @@ public:
|
|||
bool
|
||||
run() const override
|
||||
{
|
||||
IntSequence bnvs(2); bnvs[0] = 10; bnvs[1] = 7;
|
||||
IntSequence bnvs{10, 7};
|
||||
return dense_prod(Symmetry{2, 3}, bnvs, 3, 15, 10);
|
||||
}
|
||||
};
|
||||
|
|
@ -908,7 +898,7 @@ public:
|
|||
bool
|
||||
run() const override
|
||||
{
|
||||
IntSequence bnvs(2); bnvs[0] = 13; bnvs[1] = 11;
|
||||
IntSequence bnvs{13, 11};
|
||||
return dense_prod(Symmetry{3, 2}, bnvs, 3, 20, 20);
|
||||
}
|
||||
};
|
||||
|
|
@ -1116,86 +1106,83 @@ public:
|
|||
int
|
||||
main()
|
||||
{
|
||||
TestRunnable *all_tests[50];
|
||||
std::vector<std::unique_ptr<TestRunnable>> all_tests;
|
||||
// fill in vector of all tests
|
||||
int num_tests = 0;
|
||||
all_tests[num_tests++] = new SmallIndexForwardFold();
|
||||
all_tests[num_tests++] = new SmallIndexForwardUnfold();
|
||||
all_tests[num_tests++] = new IndexForwardFold();
|
||||
all_tests[num_tests++] = new IndexForwardUnfold();
|
||||
all_tests[num_tests++] = new SmallIndexBackwardFold();
|
||||
all_tests[num_tests++] = new IndexBackwardFold();
|
||||
all_tests[num_tests++] = new SmallIndexBackwardUnfold();
|
||||
all_tests[num_tests++] = new IndexBackwardUnfold();
|
||||
all_tests[num_tests++] = new SmallIndexOffsetFold();
|
||||
all_tests[num_tests++] = new SmallIndexOffsetUnfold();
|
||||
all_tests[num_tests++] = new IndexOffsetFold();
|
||||
all_tests[num_tests++] = new IndexOffsetUnfold();
|
||||
all_tests[num_tests++] = new SmallFoldUnfoldFS();
|
||||
all_tests[num_tests++] = new SmallFoldUnfoldGS();
|
||||
all_tests[num_tests++] = new FoldUnfoldFS();
|
||||
all_tests[num_tests++] = new FoldUnfoldGS();
|
||||
all_tests[num_tests++] = new SmallFoldUnfoldR();
|
||||
all_tests[num_tests++] = new FoldUnfoldR();
|
||||
all_tests[num_tests++] = new SmallDenseProd();
|
||||
all_tests[num_tests++] = new DenseProd();
|
||||
all_tests[num_tests++] = new BigDenseProd();
|
||||
all_tests[num_tests++] = new SmallFoldedMonomial();
|
||||
all_tests[num_tests++] = new FoldedMonomial();
|
||||
all_tests[num_tests++] = new SmallUnfoldedMonomial();
|
||||
all_tests[num_tests++] = new UnfoldedMonomial();
|
||||
all_tests[num_tests++] = new FoldedContractionSmall();
|
||||
all_tests[num_tests++] = new FoldedContractionBig();
|
||||
all_tests[num_tests++] = new UnfoldedContractionSmall();
|
||||
all_tests[num_tests++] = new UnfoldedContractionBig();
|
||||
all_tests[num_tests++] = new PolyEvalSmall();
|
||||
all_tests[num_tests++] = new PolyEvalBig();
|
||||
all_tests[num_tests++] = new FoldZContSmall();
|
||||
all_tests[num_tests++] = new FoldZCont();
|
||||
all_tests[num_tests++] = new UnfoldZContSmall();
|
||||
all_tests[num_tests++] = new UnfoldZCont();
|
||||
all_tests.push_back(std::make_unique<SmallIndexForwardFold>());
|
||||
all_tests.push_back(std::make_unique<SmallIndexForwardUnfold>());
|
||||
all_tests.push_back(std::make_unique<IndexForwardFold>());
|
||||
all_tests.push_back(std::make_unique<IndexForwardUnfold>());
|
||||
all_tests.push_back(std::make_unique<SmallIndexBackwardFold>());
|
||||
all_tests.push_back(std::make_unique<IndexBackwardFold>());
|
||||
all_tests.push_back(std::make_unique<SmallIndexBackwardUnfold>());
|
||||
all_tests.push_back(std::make_unique<IndexBackwardUnfold>());
|
||||
all_tests.push_back(std::make_unique<SmallIndexOffsetFold>());
|
||||
all_tests.push_back(std::make_unique<SmallIndexOffsetUnfold>());
|
||||
all_tests.push_back(std::make_unique<IndexOffsetFold>());
|
||||
all_tests.push_back(std::make_unique<IndexOffsetUnfold>());
|
||||
all_tests.push_back(std::make_unique<SmallFoldUnfoldFS>());
|
||||
all_tests.push_back(std::make_unique<SmallFoldUnfoldGS>());
|
||||
all_tests.push_back(std::make_unique<FoldUnfoldFS>());
|
||||
all_tests.push_back(std::make_unique<FoldUnfoldGS>());
|
||||
all_tests.push_back(std::make_unique<SmallFoldUnfoldR>());
|
||||
all_tests.push_back(std::make_unique<FoldUnfoldR>());
|
||||
all_tests.push_back(std::make_unique<SmallDenseProd>());
|
||||
all_tests.push_back(std::make_unique<DenseProd>());
|
||||
all_tests.push_back(std::make_unique<BigDenseProd>());
|
||||
all_tests.push_back(std::make_unique<SmallFoldedMonomial>());
|
||||
all_tests.push_back(std::make_unique<FoldedMonomial>());
|
||||
all_tests.push_back(std::make_unique<SmallUnfoldedMonomial>());
|
||||
all_tests.push_back(std::make_unique<UnfoldedMonomial>());
|
||||
all_tests.push_back(std::make_unique<FoldedContractionSmall>());
|
||||
all_tests.push_back(std::make_unique<FoldedContractionBig>());
|
||||
all_tests.push_back(std::make_unique<UnfoldedContractionSmall>());
|
||||
all_tests.push_back(std::make_unique<UnfoldedContractionBig>());
|
||||
all_tests.push_back(std::make_unique<PolyEvalSmall>());
|
||||
all_tests.push_back(std::make_unique<PolyEvalBig>());
|
||||
all_tests.push_back(std::make_unique<FoldZContSmall>());
|
||||
all_tests.push_back(std::make_unique<FoldZCont>());
|
||||
all_tests.push_back(std::make_unique<UnfoldZContSmall>());
|
||||
all_tests.push_back(std::make_unique<UnfoldZCont>());
|
||||
|
||||
// find maximum dimension and maximum nvar
|
||||
int dmax = 0;
|
||||
int nvmax = 0;
|
||||
for (int i = 0; i < num_tests; i++)
|
||||
for (const auto &test : all_tests)
|
||||
{
|
||||
if (dmax < all_tests[i]->dim)
|
||||
dmax = all_tests[i]->dim;
|
||||
if (nvmax < all_tests[i]->nvar)
|
||||
nvmax = all_tests[i]->nvar;
|
||||
if (dmax < test->dim)
|
||||
dmax = test->dim;
|
||||
if (nvmax < test->nvar)
|
||||
nvmax = test->nvar;
|
||||
}
|
||||
tls.init(dmax, nvmax); // initialize library
|
||||
|
||||
// launch the tests
|
||||
int success = 0;
|
||||
for (int i = 0; i < num_tests; i++)
|
||||
for (const auto &test : all_tests)
|
||||
{
|
||||
try
|
||||
{
|
||||
if (all_tests[i]->test())
|
||||
if (test->test())
|
||||
success++;
|
||||
}
|
||||
catch (const TLException &e)
|
||||
{
|
||||
printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
|
||||
std::cout << "Caught TL exception in <" << test->name << ">:\n";
|
||||
e.print();
|
||||
}
|
||||
catch (SylvException &e)
|
||||
{
|
||||
printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
|
||||
std::cout << "Caught Sylv exception in <" << test->name << ">:\n";
|
||||
e.printMessage();
|
||||
}
|
||||
}
|
||||
|
||||
printf("There were %d tests that failed out of %d tests run.\n",
|
||||
num_tests - success, num_tests);
|
||||
int nfailed = all_tests.size() - success;
|
||||
std::cout << "There were " << nfailed << " tests that failed out of "
|
||||
<< all_tests.size() << " tests run." << std::endl;
|
||||
|
||||
// destroy
|
||||
for (int i = 0; i < num_tests; i++)
|
||||
{
|
||||
delete all_tests[i];
|
||||
}
|
||||
|
||||
return 0;
|
||||
if (nfailed)
|
||||
return EXIT_FAILURE;
|
||||
else
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue