template based native complex trees

api/mrcpp_declarations.h

@@ -34,34 +34,34 @@ namespace mrcpp {
 
 class Timer;
 class Printer;
-template <int D> class Plotter;
+template <int D, typename T = double> class Plotter;
 
 template <int D> class Gaussian;
 template <int D> class GaussFunc;
 template <int D> class GaussPoly;
 template <int D> class GaussExp;
 
 template <int D> class BoundingBox;
-template <int D> class NodeBox;
+template <int D, typename T = double> class NodeBox;
 template <int D> class NodeIndex;
 template <int D> class NodeIndexComp;
 
-class SharedMemory;
+template <typename T = double> class SharedMemory;
 class ScalingBasis;
 class LegendreBasis;
 class InterpolatingBasis;
 
-template <int D> class RepresentableFunction;
+template <int D, typename T = double> class RepresentableFunction;
 template <int D> class MultiResolutionAnalysis;
 
-template <int D> class MWTree;
-template <int D> class FunctionTree;
+template <int D, typename T = double> class MWTree;
+template <int D, typename T = double> class FunctionTree;
 class OperatorTree;
 
-template <int D> class NodeAllocator;
+template <int D, typename T = double> class NodeAllocator;
 
-template <int D> class MWNode;
-template <int D> class FunctionNode;
+template <int D, typename T = double> class MWNode;
+template <int D, typename T = double> class FunctionNode;
 class OperatorNode;
 
 template <int D> class IdentityConvolution;
@@ -79,31 +79,30 @@ template <int D> class DerivativeKernel;
 class PoissonKernel;
 class HelmholtzKernel;
 
-template <int D> class TreeBuilder;
-template <int D> class TreeCalculator;
-template <int D> class DefaultCalculator;
-template <int D> class ProjectionCalculator;
-template <int D> class AdditionCalculator;
-template <int D> class MultiplicationCalculator;
-template <int D> class ConvolutionCalculator;
-template <int D> class DerivativeCalculator;
+template <int D, typename T = double> class TreeBuilder;
+template <int D, typename T = double> class TreeCalculator;
+template <int D, typename T = double> class DefaultCalculator;
+template <int D, typename T = double> class ProjectionCalculator;
+template <int D, typename T = double> class AdditionCalculator;
+template <int D, typename T = double> class MultiplicationCalculator;
+template <int D, typename T = double> class ConvolutionCalculator;
+template <int D, typename T = double> class DerivativeCalculator;
 class CrossCorrelationCalculator;
 
-template <int D> class TreeAdaptor;
-template <int D> class AnalyticAdaptor;
-template <int D> class WaveletAdaptor;
-template <int D> class CopyAdaptor;
+template <int D, typename T = double> class TreeAdaptor;
+template <int D, typename T = double> class AnalyticAdaptor;
+template <int D, typename T = double> class WaveletAdaptor;
+template <int D, typename T = double> class CopyAdaptor;
 
-template <int D> class TreeIterator;
-template <int D> class IteratorNode;
+template <int D, typename T = double> class TreeIterator;
+template <int D, typename T = double> class IteratorNode;
 
 class BandWidth;
-template <int D> class OperatorState;
+template <int D, typename T = double> class OperatorState;
 
 template <int D> using Coord = std::array<double, D>;
-template <int D> using MWNodeVector = std::vector<MWNode<D> *>;
+template <int D, typename T = double> using MWNodeVector = std::vector<MWNode<D, T> *>;
 
-template <typename T, typename U> using FMap_ = std::function<T(U)>;
-typedef FMap_<double, double> FMap;
+template <typename T = double, typename U = double> using FMap = std::function<T(U)>;
 
 } // namespace mrcpp

examples/derivative.cpp

@@ -51,8 +51,8 @@ int main(int argc, char **argv) {
     mrcpp::FunctionTree<D> err_tree(MRA);
 
     // Projecting functions
-    mrcpp::project<D>(prec, f_tree, f);
-    mrcpp::project<D>(prec, df_tree, df);
+    mrcpp::project<D, double>(prec, f_tree, f);
+    mrcpp::project<D, double>(prec, df_tree, df);
 
     // Applying derivative operator
     mrcpp::apply(dg_tree, D_00, f_tree, 0);

examples/mpi_matrix.cpp

@@ -54,7 +54,7 @@ int main(int argc, char **argv) {
         };
         mrcpp::FunctionTree<3> *tree = new mrcpp::FunctionTree<3>(MRA);
         if (i % wsize == wrank) {
-            mrcpp::project<3>(prec, *tree, f);
+	  mrcpp::project<3, double>(prec, *tree, f);
             tree->normalize();
         }
         f_vec.push_back(std::make_tuple(1.0, tree));

examples/mpi_send_tree.cpp

@@ -55,7 +55,7 @@ int main(int argc, char **argv) {
     mrcpp::FunctionTree<D> f_tree(MRA);
 
     // Only rank 0 projects the function
-    if (wrank == 0) mrcpp::project<D>(prec, f_tree, f);
+    if (wrank == 0) mrcpp::project<D, double>(prec, f_tree, f);
 
     { // Print data before send
         auto integral = f_tree.integrate();

examples/mpi_shared_tree.cpp

@@ -63,12 +63,12 @@ int main(int argc, char **argv) {
     };
 
     // Initialize a shared memory tree, max 100MB
-    auto shared_mem = new mrcpp::SharedMemory(scomm, 100);
+    auto shared_mem = new mrcpp::SharedMemory<double>(scomm, 100);
     mrcpp::FunctionTree<D> f_tree(MRA, shared_mem);
 
     // Only first rank projects
     auto frank = 0;
-    if (srank == frank) mrcpp::project<D>(prec, f_tree, f);
+    if (srank == frank) mrcpp::project<D, double>(prec, f_tree, f);
     mrcpp::share_tree(f_tree, frank, 0, scomm);
 
     { // Print data after share

examples/projection.cpp

@@ -37,7 +37,7 @@ int main(int argc, char **argv) {
 
     // Projecting function
     mrcpp::FunctionTree<D> f_tree(MRA);
-    mrcpp::project<D>(prec, f_tree, f, -1);
+    mrcpp::project<D, double>(prec, f_tree, f, -1);
     auto integral = f_tree.integrate();
 
     mrcpp::print::header(0, "Projecting analytic function");

examples/scf.cpp

@@ -31,7 +31,7 @@ void setupNuclearPotential(double Z, FunctionTree<D> &V) {
     };
 
     // Projecting function
-    project<D>(prec, V, f);
+    project<D, double>(prec, V, f);
 
     print::footer(0, timer, 2);
     Printer::setPrintLevel(oldlevel);
@@ -48,7 +48,7 @@ void setupInitialGuess(FunctionTree<D> &phi) {
     };
 
     // Projecting and normalizing function
-    project<D>(prec, phi, f);
+    project<D, double>(prec, phi, f);
     phi.normalize();
 
     print::footer(0, timer, 2);

examples/schrodinger_semigroup1d.cpp

@@ -93,13 +93,13 @@ int main(int argc, char **argv)
 
     // Projecting functions
     mrcpp::FunctionTree<1> Re_f_tree(MRA);
-    mrcpp::project<1>(prec, Re_f_tree, Re_f);
+    mrcpp::project<1, double>(prec, Re_f_tree, Re_f);
     mrcpp::FunctionTree<1> Im_f_tree(MRA);
-    mrcpp::project<1>(prec, Im_f_tree, Im_f);
+    mrcpp::project<1, double>(prec, Im_f_tree, Im_f);
     mrcpp::FunctionTree<1> Re_g_tree(MRA);
-    mrcpp::project<1>(prec, Re_g_tree, Re_g);
+    mrcpp::project<1, double>(prec, Re_g_tree, Re_g);
     mrcpp::FunctionTree<1> Im_g_tree(MRA);
-    mrcpp::project<1>(prec, Im_g_tree, Im_g);
+    mrcpp::project<1, double>(prec, Im_g_tree, Im_g);
 
     // Output function trees
     mrcpp::FunctionTree<1> Re_fout_tree(MRA);

examples/tree_cleaner.cpp

@@ -9,6 +9,7 @@ const auto order = 7;
 const auto prec = 1.0e-5;
 
 const auto D = 3;
+
 int main(int argc, char **argv) {
     auto timer = mrcpp::Timer();
 
@@ -42,14 +43,14 @@ int main(int argc, char **argv) {
     auto iter = 0;
     auto n_nodes = 1;
     while (n_nodes > 0) {
-        mrcpp::project<D>(-1.0, f_tree, f);         // Projecting on fixed grid
+      mrcpp::project<D, double>(-1.0, f_tree, f);         // Projecting on fixed grid
         n_nodes = mrcpp::refine_grid(f_tree, prec); // Refine grid
         mrcpp::clear_grid(f_tree);                  // Clear MW coefs
         printout(0, " iter " << std::setw(3) << iter++ << std::setw(45));
         printout(0, " n_nodes " << std::setw(5) << n_nodes << std::endl);
     }
     // Projecting on final converged grid
-    mrcpp::project<D>(-1.0, f_tree, f);
+    mrcpp::project<D, double>(-1.0, f_tree, f);
 
     auto integral = f_tree.integrate();
     auto sq_norm = f_tree.getSquareNorm();

src/functions/AnalyticFunction.h

@@ -32,29 +32,29 @@
 
 namespace mrcpp {
 
-template <int D> class AnalyticFunction : public RepresentableFunction<D> {
+template <int D, typename T = double> class AnalyticFunction : public RepresentableFunction<D, T> {
 public:
     AnalyticFunction() = default;
     ~AnalyticFunction() override = default;
 
-    AnalyticFunction(std::function<double(const Coord<D> &r)> f, const double *a = nullptr, const double *b = nullptr)
-            : RepresentableFunction<D>(a, b)
+    AnalyticFunction(std::function<T(const Coord<D> &r)> f, const double *a = nullptr, const double *b = nullptr)
+            : RepresentableFunction<D, T>(a, b)
             , func(f) {}
-    AnalyticFunction(std::function<double(const Coord<D> &r)> f,
+    AnalyticFunction(std::function<T(const Coord<D> &r)> f,
                      const std::vector<double> &a,
                      const std::vector<double> &b)
             : AnalyticFunction(f, a.data(), b.data()) {}
 
-    void set(std::function<double(const Coord<D> &r)> f) { this->func = f; }
+    void set(std::function<T(const Coord<D> &r)> f) { this->func = f; }
 
-    double evalf(const Coord<D> &r) const override {
-        double val = 0.0;
+    T evalf(const Coord<D> &r) const override {
+        T val = 0.0;
         if (not this->outOfBounds(r)) val = this->func(r);
         return val;
     }
 
 protected:
-    std::function<double(const Coord<D> &r)> func;
+    std::function<T(const Coord<D> &r)> func;
 };
 
 } // namespace mrcpp

src/functions/BoysFunction.cpp

@@ -32,7 +32,7 @@
 namespace mrcpp {
 
 BoysFunction::BoysFunction(int n, double p)
-        : RepresentableFunction<1>()
+  : RepresentableFunction<1, double>()
         , order(n)
         , prec(p)
         , MRA(BoundingBox<1>(), InterpolatingBasis(13)) {}
@@ -50,8 +50,8 @@ double BoysFunction::evalf(const Coord<1> &r) const {
         return std::exp(-xt_2) * t_2n;
     };
 
-    FunctionTree<1> tree(this->MRA);
-    mrcpp::project<1>(this->prec, tree, f);
+    FunctionTree<1, double> tree(this->MRA);
+    mrcpp::project<1, double>(this->prec, tree, f);
     double result = tree.integrate();
 
     Printer::setPrintLevel(oldlevel);

src/functions/BoysFunction.h

@@ -30,7 +30,7 @@
 
 namespace mrcpp {
 
-class BoysFunction final : public RepresentableFunction<1> {
+  class BoysFunction final : public RepresentableFunction<1, double> {
 public:
     BoysFunction(int n, double prec = 1.0e-10);
 

src/functions/GaussExp.h

@@ -51,7 +51,7 @@ namespace mrcpp {
  *
  */
 
-template <int D> class GaussExp : public RepresentableFunction<D> {
+    template <int D> class GaussExp : public RepresentableFunction<D, double> {
 public:
     GaussExp(int nTerms = 0, double prec = GAUSS_EXP_PREC);
     GaussExp(const GaussExp<D> &gExp);

src/functions/Gaussian.h

@@ -40,7 +40,7 @@
 
 namespace mrcpp {
 
-template <int D> class Gaussian : public RepresentableFunction<D> {
+    template <int D> class Gaussian : public RepresentableFunction<D, double> {
 public:
     Gaussian(double a, double c, const Coord<D> &r, const std::array<int, D> &p);
     Gaussian(const std::array<double, D> &a, double c, const Coord<D> &r, const std::array<int, D> &p);

src/functions/Polynomial.cpp

@@ -45,7 +45,7 @@ namespace mrcpp {
 /** Construct polynomial of order zero with given size and bounds.
  * Includes default constructor. */
 Polynomial::Polynomial(int k, const double *a, const double *b)
-        : RepresentableFunction<1>(a, b) {
+  : RepresentableFunction<1, double>(a, b) {
     assert(k >= 0);
     this->N = 1.0;
     this->L = 0.0;

src/functions/Polynomial.h

@@ -44,7 +44,7 @@
 
 namespace mrcpp {
 
-class Polynomial : public RepresentableFunction<1> {
+  class Polynomial : public RepresentableFunction<1, double> {
 public:
     Polynomial(int k = 0, const double *a = nullptr, const double *b = nullptr);
     Polynomial(int k, const std::vector<double> &a, const std::vector<double> &b)

src/functions/RepresentableFunction.cpp

@@ -38,7 +38,7 @@
 
 namespace mrcpp {
 
-template <int D> RepresentableFunction<D>::RepresentableFunction(const double *a, const double *b) {
+template <int D, typename T> RepresentableFunction<D, T>::RepresentableFunction(const double *a, const double *b) {
     if (a == nullptr or b == nullptr) {
         this->bounded = false;
         this->A = nullptr;
@@ -56,7 +56,7 @@ template <int D> RepresentableFunction<D>::RepresentableFunction(const double *a
 }
 
 /** Constructs a new function with same bounds as the input function */
-template <int D> RepresentableFunction<D>::RepresentableFunction(const RepresentableFunction<D> &func) {
+template <int D, typename T> RepresentableFunction<D, T>::RepresentableFunction(const RepresentableFunction<D, T> &func) {
     if (func.isBounded()) {
         this->bounded = true;
         this->A = new double[D];
@@ -74,11 +74,11 @@ template <int D> RepresentableFunction<D>::RepresentableFunction(const Represent
 
 /** Copies function, not bounds. Use copy constructor if you want an
  * identical function. */
-template <int D> RepresentableFunction<D> &RepresentableFunction<D>::operator=(const RepresentableFunction<D> &func) {
+template <int D, typename T> RepresentableFunction<D, T> &RepresentableFunction<D, T>::operator=(const RepresentableFunction<D, T> &func) {
     return *this;
 }
 
-template <int D> RepresentableFunction<D>::~RepresentableFunction() {
+template <int D, typename T> RepresentableFunction<D, T>::~RepresentableFunction() {
     if (this->isBounded()) {
         delete[] this->A;
         delete[] this->B;
@@ -87,7 +87,7 @@ template <int D> RepresentableFunction<D>::~RepresentableFunction() {
     this->B = nullptr;
 }
 
-template <int D> void RepresentableFunction<D>::setBounds(const double *a, const double *b) {
+template <int D, typename T> void RepresentableFunction<D, T>::setBounds(const double *a, const double *b) {
     if (a == nullptr or b == nullptr) { MSG_ERROR("Invalid arguments"); }
     if (not isBounded()) {
         this->bounded = true;
@@ -101,7 +101,7 @@ template <int D> void RepresentableFunction<D>::setBounds(const double *a, const
     }
 }
 
-template <int D> bool RepresentableFunction<D>::outOfBounds(const Coord<D> &r) const {
+template <int D, typename T> bool RepresentableFunction<D, T>::outOfBounds(const Coord<D> &r) const {
     if (not isBounded()) { return false; }
     for (int d = 0; d < D; d++) {
         if (r[d] < getLowerBound(d)) return true;
@@ -110,8 +110,11 @@ template <int D> bool RepresentableFunction<D>::outOfBounds(const Coord<D> &r) c
     return false;
 }
 
-template class RepresentableFunction<1>;
-template class RepresentableFunction<2>;
-template class RepresentableFunction<3>;
+template class RepresentableFunction<1, double>;
+template class RepresentableFunction<2, double>;
+template class RepresentableFunction<3, double>;
+template class RepresentableFunction<1, ComplexDouble>;
+template class RepresentableFunction<2, ComplexDouble>;
+template class RepresentableFunction<3, ComplexDouble>;
 
 } // namespace mrcpp