From 6ac7fef588a8ffcd9c1b8174e0a5f60aaac2b605 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Mon, 9 May 2016 17:04:34 -0400 Subject: [PATCH] fixed multiply --- Matrix.h | 19 ++++-- function1d.h | 149 +++++++++++++++++++++++++++++++++++------------ test_function.cc | 53 ++++++++++++++++- 3 files changed, 175 insertions(+), 46 deletions(-) diff --git a/Matrix.h b/Matrix.h index 90c2c58..9834e84 100644 --- a/Matrix.h +++ b/Matrix.h @@ -44,6 +44,7 @@ class VectorT { // is the matrix currently allocated? bool _allocated; + // this can serve as either a deep copy or as a zeros_like (in Numpy) friend VectorT copy(const VectorT& t, bool empty = false) { VectorT r; if (t._allocated) { @@ -94,9 +95,9 @@ class VectorT { for (auto i = 0; i < v.size(); i++) _p[i] = v[i]; } - VectorT(int k) { + VectorT(int k, T val = T(0)) { create(k); - for (auto i = 0; i < k; i++) _p[i] = T(0); + for (auto i = 0; i < k; i++) _p[i] = val; } void create(int d0) { @@ -176,6 +177,14 @@ class VectorT { return r; } + // point by point multiplication + VectorT operator*(const VectorT& t) const { + assert(_dim0 == t._dim0); + VectorT r = copy(t, true); + for (auto i = 0; i < size(); i++) r._p[i] = _p[i]*t._p[i]; + return r; + } + // addition VectorT operator+(const VectorT& t) const { return gaxpy_oop(1.0, t, 1.0); @@ -525,14 +534,12 @@ class MatrixT { // multiply with a VectorT VectorT operator* (const VectorT& v) const { int vsz = v.size(); - // VectorT r(vsz); - VectorT r; - r.create(vsz); + VectorT r(_dim0); assert(vsz == _dim1); for (int i = 0; i < _dim0; i++) { T s = T(0); for (int j = 0; j < _dim1; j++) { - s += _p[i*_dim0+j]*v[j]; + s += _p[i*_dim1+j]*v[j]; } r[i] = s; } diff --git a/function1d.h b/function1d.h index 1baf9b2..beced52 100644 --- a/function1d.h +++ b/function1d.h @@ -13,6 +13,28 @@ double normf(Vector v) { return std::sqrt(s); } +inline +void pack(int k, const Vector& s0, const Vector& s1, Vector& s) { + assert(s0.size() == k); + assert(s1.size() == k); + assert(s.size() == 2*k); + for (auto i = 0; i < k; i++) { + s[i] = s0[i]; + s[i+k] = s1[i]; + } +} + +inline +void unpack(int k, const Vector& s, Vector& s0, Vector& s1) { + assert(s0.size() == k); + assert(s1.size() == k); + assert(s.size() == 2*k); + for (auto i = 0; i < k; i++) { + s0[i] = s[i]; + s1[i] = s[i+k]; + } +} + class Function1D { private: enum class FunctionForm {RECONSTRUCTED, COMPRESSED, UNDEFINED}; @@ -33,6 +55,7 @@ class Function1D { real_matrix quad_phi; real_matrix quad_phiT; real_matrix quad_phiw; + real_matrix quad_phiwT; public: // I know that friends are evil, but I decided to have them anyway @@ -66,8 +89,8 @@ class Function1D { ~Function1D() {} - bool is_compressed() {return form == FunctionForm::COMPRESSED;} - bool is_reconstructed() {return form == FunctionForm::RECONSTRUCTED;} + bool is_compressed() const {return form == FunctionForm::COMPRESSED;} + bool is_reconstructed() const {return form == FunctionForm::RECONSTRUCTED;} void init_twoscale(int k) { hg = TwoScaleCoeffs::instance()->hg(k); @@ -79,9 +102,10 @@ class Function1D { quad_w = gauss_legendre_w(order); quad_npts = quad_w.size(); - quad_phi = zeros(quad_npts, k); - quad_phiT = zeros(k, quad_npts); - quad_phiw = zeros(quad_npts, k); + quad_phi = zeros(quad_npts, k); + quad_phiT = zeros(k, quad_npts); + quad_phiw = zeros(quad_npts, k); + quad_phiwT = zeros(k, quad_npts); for (auto i = 0; i < quad_npts; i++) { auto p = ScalingFunction::instance()->phi(quad_x[i], k); for (auto m = 0; m < k; m++) { @@ -90,6 +114,7 @@ class Function1D { quad_phiw(i,m) = quad_w[i]*p[m]; } } + quad_phiwT = transpose(quad_phiw); } Vector project_box(double (*f)(double), int n, int l) { @@ -110,16 +135,13 @@ class Function1D { Vector s0 = project_box(f, n+1, 2*l); Vector s1 = project_box(f, n+1, 2*l+1); Vector s(2*k); - for (auto i = 0; i < k; i++) { - s[i] = s0[i]; - s[i+k] = s1[i]; - } + pack(k, s0, s1, s); Vector d = hg*s; + stree[Key(n,l)] = Vector(); if (normf(Vector(d.slice(k,2*k-1))) < thresh || n >= maxlevel-1) { stree[Key(n+1,2*l)] = s0; stree[Key(n+1,2*l+1)] = s1; } else { - stree[Key(n,l)] = Vector(); refine(f, n+1, 2*l); refine(f, n+1, 2*l+1); } @@ -130,10 +152,7 @@ class Function1D { if (dp != dtree.end()) { Vector dd = dp->second; Vector d(2*k); - for (auto i = 0; i < k; i++) { - d[i] = ss[i]; - d[i+k] = dd[i]; - } + pack(k, ss, dd, d); auto s = hgT*d; Vector s0(k); Vector s1(k); @@ -165,10 +184,7 @@ class Function1D { s1 = compress_spawn(dtree_r, n+1, 2*l+1); } Vector s(2*k); - for (auto i = 0; i < k; i++) { - s[i] = s0[i]; - s[i+k] = s1[i]; - } + pack(k, s0, s1, s); Vector d = hg*s; auto sr = d.slice(0,k-1); auto dr = d.slice(k,2*k-1); @@ -220,23 +236,70 @@ class Function1D { return r; } - // void mul_helper(const CoeffTree& f, const CoeffTree& g, const Vector& fs, const Vector& gs, int n, int l) { - // auto fsin = fs; - // if (fsin.size() == 0) { - // const auto fp = f.find(Key(n,l)); - // - // } - // const auto gp = g.find(Key(n,l)); - // - // } + void mul_helper(CoeffTree& r, const CoeffTree& f, const CoeffTree& g, + const Vector& fsin, const Vector& gsin, int n, int l) { + auto mrefine = true; + auto fs = fsin; + if (fs.size() == 0) { + const auto fp = f.find(Key(n,l)); + assert(fp != f.end()); + fs = fp->second; + } + auto gs = gsin; + if (gs.size() == 0) { + const auto gp = g.find(Key(n,l)); + assert(gp != g.end()); + gs = gp->second; + } + if (fs.size() && gs.size()) { + if (mrefine) { + // refine to lower level for both f and g + Vector fd(2*k); + Vector gd(2*k); + // pack the coeffs together so that we can do two scale + pack(k, fs, Vector(k, 0.0), fd); + pack(k, gs, Vector(k, 0.0), gd); + auto fss = hgT*fd; auto gss = hgT*gd; + Vector fs0(k); Vector fs1(k); + Vector gs0(k); Vector gs1(k); + // unpack the coeffs on n+1 level + unpack(k, fss, fs0, fs1); + unpack(k, gss, gs0, gs1); + // convert to function values + auto scale = std::sqrt(std::pow(2.0, n+1)); + auto fs0vals = quad_phi * fs0; + auto fs1vals = quad_phi * fs1; + auto gs0vals = quad_phi * gs0; + auto gs1vals = quad_phi * gs1; + auto rs0 = fs0vals * gs0vals; + auto rs1 = fs1vals * gs1vals; + rs0 = (quad_phiwT * rs0); + rs1 = (quad_phiwT * rs1); + rs0.scale(scale); + rs1.scale(scale); + r[Key(n,l)] = Vector(); + r[Key(n+1,2*l)] = rs0; + r[Key(n+1,2*l+1)] = rs1; + } else { + // do multiply + auto rs = fs * gs; + r[Key(n,l)] = rs; + } + } else { + r[Key(n,l)] = Vector(); + mul_helper(r, f, g, fs, gs, n+1, 2*l); + mul_helper(r, f, g, fs, gs, n+1, 2*l+1); + } + } - // Function1D operator*(const Function1D& f) const { - // Function1D r(f.k, f.thresh, f.maxlevel, f.initiallevel); - // r.form = Function1D::FunctionForm::RECONSTRUCTED; - // assert(is_reconstructed()); - // assert(f.is_reconstructed()); - // r.mul_helper(stree, f.stree, Vector(), Vector(), 0, 0); - // } + Function1D operator*(const Function1D& g) const { + Function1D r(g.k, g.thresh, g.maxlevel, g.initiallevel); + r.form = Function1D::FunctionForm::RECONSTRUCTED; + assert(is_reconstructed()); + assert(g.is_reconstructed()); + r.mul_helper(r.stree, stree, g.stree, Vector(), Vector(), 0, 0); + return r; + } double eval(double x, int n, int l) { assert(n < maxlevel); @@ -274,18 +337,28 @@ class Function1D { printf(") %15.8e\n",normf(d)); } - void print_tree() { + void print_tree(bool docoeffs = false) { printf("sum coeffs:\n"); for (auto c : stree) { auto k = c.first; - auto s = c.second; - printf("[%d %d] %15.8e\n", k.n, k.l, normf(s)); + auto s = c.second; + if (docoeffs) { + printf("[%d %d] %15.8e\n", k.n, k.l, normf(s)); + print(s); + } else { + printf("[%d %d] %15.8e\n", k.n, k.l, normf(s)); + } } printf("diff coeffs:\n"); for (auto c : dtree) { auto k = c.first; auto d = c.second; - printf("[%d %d] %15.8e\n", k.n, k.l, normf(d)); + if (docoeffs) { + printf("[%d %d] %15.8e\n", k.n, k.l, normf(d)); + print(d); + } else { + printf("[%d %d] %15.8e\n", k.n, k.l, normf(d)); + } } } }; diff --git a/test_function.cc b/test_function.cc index 2470af3..f86d5f0 100644 --- a/test_function.cc +++ b/test_function.cc @@ -3,7 +3,15 @@ const double PI = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825; const int k = 4; const double thresh = 1e-8; -const int initiallevel = 2; +const int initiallevel = 1; + +double linear(double x) { + return x; +} + +double quadratic(double x) { + return x*x; +} double func1(double x) { auto c = 1.0; @@ -21,6 +29,10 @@ double sum_f12(double x) { return func1(x) + func2(x); } +double mul_f12(double x) { + return func1(x) * func2(x); +} + void test_function_point() { Function1D f(func1, k, thresh, 30, initiallevel); f.print_tree(); @@ -56,8 +68,45 @@ void test_add() { f12c.print_tree(); } +void test_mul() { + auto fr = Function1D(func1, k, thresh, 30, initiallevel); + printf("\nfunction fr:\n"); + fr.print_tree(); + auto gr = Function1D(func2, k, thresh, 30, initiallevel); + printf("\nfunction gr:\n"); + gr.print_tree(); + auto fgr = fr * gr; + printf("\nfunction f+g(r):\n"); + fgr.print_tree(); + auto f12r = Function1D(mul_f12, k, thresh, 30, initiallevel+1); + printf("\nfunction f12r:\n"); + f12r.print_tree(); + + auto x = 0.23111; + printf("x: %15.8e f: %15.8e func: %15.8e error: %15.8e\n", x, fgr(x), f12r(x), std::abs(fgr(x)-f12r(x))); +} + +void test_simple_mul() { + auto fr = Function1D(linear, k, thresh, 30, 0); + printf("\nfunction fr:\n"); + fr.print_tree(true); + auto gr = Function1D(linear, k, thresh, 30, 0); + printf("\nfunction gr:\n"); + gr.print_tree(true); + auto fgr = fr * gr; + printf("\nfunction f+g(r):\n"); + fgr.print_tree(true); + auto f12r = Function1D(quadratic, k, thresh, 30, 1); + printf("\nfunction f12r:\n"); + f12r.print_tree(true); + + auto x = 0.23111; + printf("x: %15.8e f: %15.8e func: %15.8e error: %15.8e\n", x, fgr(x), f12r(x), std::abs(fgr(x)-f12r(x))); +} + int main(int argc, char** argv) { //test_function_compress(); - test_add(); + //test_add(); + test_mul(); return 0; }