[pull] master from MRChemSoft:master

api/MWFunctions

@@ -27,7 +27,7 @@
 #include "trees/BoundingBox.h"
 #include "trees/FunctionTree.h"
 #include "trees/FunctionTreeVector.h"
-#include "utils/ComplexFunction.h"
+#include "utils/CompFunction.h"
 
 #include "core/InterpolatingBasis.h"
 #include "core/LegendreBasis.h"

api/mrcpp_declarations.h

@@ -34,34 +34,36 @@ 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;
+template <int D = 3> class CompFunction;
+class ComplexFunction;
 class OperatorNode;
 
 template <int D> class IdentityConvolution;
@@ -79,31 +81,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

cmake/compiler_flags/CXXFlags.cmake

@@ -13,7 +13,7 @@
 
 option(ENABLE_ARCH_FLAGS "Enable architecture-specific compiler flags" ON)
 
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED TRUE)
 set(CMAKE_CXX_EXTENSIONS FALSE)
 set(CMAKE_EXPORT_COMPILE_COMMANDS TRUE)

docs/mrcpp_api/mwfunctions.rst

@@ -165,7 +165,7 @@ Constructing an MRA
 
 An MRA is defined in two steps, first the computational domain is given by a
 ``BoundingBox`` (D is the dimension), e.g. for a total domain of
-:math:`[-32,32]^3` in three dimensions (eight root boxes of size :math:`[16]^3`
+:math:`[-32,32]^3` in three dimensions (eight root boxes of size :math:`[32]^3`
 each):
 
 .. code-block:: cpp

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

@@ -21,17 +21,14 @@ void setupNuclearPotential(double Z, FunctionTree<D> &V) {
 
     // Smoothing parameter
     auto c = 0.00435 * prec / std::pow(Z, 5);
-    auto u = [](double r) -> double {
-        return std::erf(r) / r +
-               1.0 / (3.0 * std::sqrt(mrcpp::pi)) * (std::exp(-r * r) + 16.0 * std::exp(-4.0 * r * r));
-    };
+    auto u = [](double r) -> double { return std::erf(r) / r + 1.0 / (3.0 * std::sqrt(mrcpp::pi)) * (std::exp(-r * r) + 16.0 * std::exp(-4.0 * r * r)); };
     auto f = [u, c, Z](const Coord<3> &r) -> double {
         auto x = std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
         return -1.0 * Z * u(x / c) / c;
     };
 
     // Projecting function
-    project<D>(prec, V, f);
+    project<D, double>(prec, V, f);
 
     print::footer(0, timer, 2);
     Printer::setPrintLevel(oldlevel);
@@ -48,7 +45,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

@@ -1,29 +1,26 @@
 #include "MRCPP/MWFunctions"
-#include <MRCPP/MWOperators>
-#include <MRCPP/Printer>
 #include "MRCPP/Plotter"
-#include <MRCPP/Timer>
-#include "operators/TimeEvolutionOperator.h"
 #include "functions/special_functions.h"
+#include "operators/TimeEvolutionOperator.h"
 #include "treebuilders/complex_apply.h"
-
-
+#include <MRCPP/MWOperators>
+#include <MRCPP/Printer>
+#include <MRCPP/Timer>
 
 const auto min_scale = 0;
 const auto max_depth = 25;
 
 const auto order = 4;
 const auto prec = 1.0e-7;
 
-int finest_scale = 10; //for time evolution operator construction (not recommended to use more than 10)
-int max_Jpower = 20;  //the amount of J integrals to be used in construction (20 should be enough)
+int finest_scale = 10; // for time evolution operator construction (not recommended to use more than 10)
+int max_Jpower = 20;   // the amount of J integrals to be used in construction (20 should be enough)
 
 // Time moments:
-double t1 = 0.001;         //initial time moment (not recommended to use more than 0.001)
-double delta_t = 0.001;    //time step (not recommended to use less than 0.001)
-double t2 = delta_t + t1;  //final time moment
+double t1 = 0.001;        // initial time moment (not recommended to use more than 0.001)
+double delta_t = 0.001;   // time step (not recommended to use less than 0.001)
+double t2 = delta_t + t1; // final time moment
 
-
 /**
  * @brief Exploring free-particle time evolution.
  * @details We check the time propagator.
@@ -41,17 +38,16 @@ double t2 = delta_t + t1;  //final time moment
  *   \psi(x, t) = \sqrt{\frac{\sigma}{4it + \sigma}} e^{-\frac{(x - x_0)^2}{4it + \sigma}}
  *   .
  * \f]
- * 
+ *
  */
-int main(int argc, char **argv)
-{
+int main(int argc, char **argv) {
     auto timer = mrcpp::Timer();
 
     // Initialize printing
     auto printlevel = 0;
     mrcpp::Printer::init(printlevel);
     mrcpp::print::environment(0);
-    
+
     // Initialize world in the unit cube [0,1]
     auto basis = mrcpp::LegendreBasis(order);
     auto world = mrcpp::BoundingBox<1>(min_scale);
@@ -74,87 +70,74 @@ int main(int argc, char **argv)
     double x0 = 0.5;
 
     // Functions f(x) = psi(x, t1) and g(x) = psi(x, t2)
-    auto Re_f = [sigma, x0, t=t1](const mrcpp::Coord<1> &r) -> double
-    {
-        return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).real();
-    };
-    auto Im_f = [sigma, x0, t=t1](const mrcpp::Coord<1> &r) -> double
-    {
-        return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).imag();
-    };
-    auto Re_g = [sigma, x0, t=t2](const mrcpp::Coord<1> &r) -> double
-    {
-        return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).real();
-    };
-    auto Im_g = [sigma, x0, t=t2](const mrcpp::Coord<1> &r) -> double
-    {
-        return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).imag();
-    };
+    auto Re_f = [sigma, x0, t = t1](const mrcpp::Coord<1> &r) -> double { return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).real(); };
+    auto Im_f = [sigma, x0, t = t1](const mrcpp::Coord<1> &r) -> double { return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).imag(); };
+    auto Re_g = [sigma, x0, t = t2](const mrcpp::Coord<1> &r) -> double { return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).real(); };
+    auto Im_g = [sigma, x0, t = t2](const mrcpp::Coord<1> &r) -> double { return mrcpp::free_particle_analytical_solution(r[0], x0, t, sigma).imag(); };
 
     // 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);
     mrcpp::FunctionTree<1> Im_fout_tree(MRA);
-    
+
     // Complex objects for use in apply()
-    mrcpp::ComplexObject< mrcpp::ConvolutionOperator<1> > E(ReExp, ImExp);
-    mrcpp::ComplexObject< mrcpp::FunctionTree<1> > input(Re_f_tree, Im_f_tree);
-    mrcpp::ComplexObject< mrcpp::FunctionTree<1> > output(Re_fout_tree, Im_fout_tree);
+    mrcpp::ComplexObject<mrcpp::ConvolutionOperator<1>> E(ReExp, ImExp);
+    mrcpp::ComplexObject<mrcpp::FunctionTree<1>> input(Re_f_tree, Im_f_tree);
+    mrcpp::ComplexObject<mrcpp::FunctionTree<1>> output(Re_fout_tree, Im_fout_tree);
 
     mrcpp::print::header(0, "Applying operator");
     mrcpp::print::footer(0, timer, 2);
 
     // Apply operator Exp(delta_t) f(x)
     mrcpp::apply(prec, output, E, input);
-    
+
     mrcpp::print::header(0, "Checking the result on analytical solution");
     mrcpp::print::footer(0, timer, 2);
 
     // Check g(x) = Exp(delta_t) f(x)
-    mrcpp::FunctionTree<1> Re_error(MRA);  // = Re_fout_tree - Re_g_tree
-    mrcpp::FunctionTree<1> Im_error(MRA);  // = Im_fout_tree - Im_g_tree
-    
+    mrcpp::FunctionTree<1> Re_error(MRA); // = Re_fout_tree - Re_g_tree
+    mrcpp::FunctionTree<1> Im_error(MRA); // = Im_fout_tree - Im_g_tree
+
     // Re_error = Re_fout_tree - Re_g_tree
     add(prec, Re_error, 1.0, Re_fout_tree, -1.0, Re_g_tree);
     auto Re_integral = Re_error.integrate();
     auto Re_sq_norm = Re_error.getSquareNorm();
     mrcpp::print::value(0, "Integral of    Re(Exp(delta_t) f(x) - g(x)) =", Re_integral);
     mrcpp::print::value(0, "Square norm of Re(Exp(delta_t) f(x) - g(x)) =", Re_sq_norm);
-    
+
     // Im_error = Im_fout_tree - Im_g_tree
     add(prec, Im_error, 1.0, Im_fout_tree, -1.0, Im_g_tree);
     auto Im_integral = Im_error.integrate();
     auto Im_sq_norm = Im_error.getSquareNorm();
     mrcpp::print::value(0, "Integral of    Im(Exp(delta_t) f(x) - g(x)) =", Im_integral);
     mrcpp::print::value(0, "Square norm of Im(Exp(delta_t) f(x) - g(x)) =", Im_sq_norm);
-        
+
     mrcpp::print::header(0, "Saving plots to files");
     mrcpp::print::footer(0, timer, 2);
 
     // Set plotting parameters
-    int nPts = 1000;                                     
-    mrcpp::Coord<1> o{0.0};                   
-    mrcpp::Coord<1> a{1.0};                   
-    mrcpp::Plotter<1> plot(o);                           
+    int nPts = 1000;
+    mrcpp::Coord<1> o{0.0};
+    mrcpp::Coord<1> a{1.0};
+    mrcpp::Plotter<1> plot(o);
     plot.setRange(a);
 
-    plot.linePlot({nPts}, Re_error, "Re_error");             // Write to file Re_error.line
-    plot.linePlot({nPts}, Im_error, "Im_error");             // Write to file Im_error.line
-    plot.linePlot({nPts}, Re_f_tree, "Re_f_tree");           // Write to file Re_f_tree.line
-    plot.linePlot({nPts}, Im_f_tree, "Im_f_tree");           // Write to file Im_f_tree.line
-    plot.linePlot({nPts}, Re_g_tree, "Re_g_tree");           // Write to file Re_g_tree.line
-    plot.linePlot({nPts}, Im_g_tree, "Im_g_tree");           // Write to file Im_g_tree.line
-    
+    plot.linePlot({nPts}, Re_error, "Re_error");   // Write to file Re_error.line
+    plot.linePlot({nPts}, Im_error, "Im_error");   // Write to file Im_error.line
+    plot.linePlot({nPts}, Re_f_tree, "Re_f_tree"); // Write to file Re_f_tree.line
+    plot.linePlot({nPts}, Im_f_tree, "Im_f_tree"); // Write to file Im_f_tree.line
+    plot.linePlot({nPts}, Re_g_tree, "Re_g_tree"); // Write to file Re_g_tree.line
+    plot.linePlot({nPts}, Im_g_tree, "Im_g_tree"); // Write to file Im_g_tree.line
+
     mrcpp::print::footer(0, timer, 2);
     return 0;
 }
-