[pull] master from MRChemSoft:master

.github/workflows/build-test.yml

@@ -27,10 +27,10 @@ jobs:
         os: [ubuntu-latest]  #, macos-12]
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Cache conda
-      uses: actions/cache@v1
+      uses: actions/cache@v4
       env:
         CACHE_NUMBER: 0  # Increase this value to reset cache if .github/mrcpp-gha.yml has not changed
       with:
@@ -45,7 +45,7 @@ jobs:
           activate-environment: mrcpp-gha
           environment-file: .github/mrcpp-gha.yml
           channel-priority: true
-          python-version: 3.6
+          python-version: 3.11
           use-only-tar-bz2: true  # IMPORTANT: This needs to be set for caching to work properly!
 
     - name: Configure

.github/workflows/code-coverage.yml

@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Set up environment
       uses: conda-incubator/setup-miniconda@v3
@@ -26,7 +26,7 @@ jobs:
           activate-environment: mrcpp-codecov
           environment-file: .github/mrcpp-codecov.yml
           channel-priority: true
-          python-version: 3.6
+          python-version: 3.11
 
     - name: Configure
       shell: bash -l {0}

cmake/MRCPPConfig.cmake.in

@@ -94,7 +94,7 @@ if(NOT TARGET ${PN}::mrcpp)
   get_filename_component(_fext ${${PN}_LIBRARY} EXT)
   include("${CMAKE_CURRENT_LIST_DIR}/${PN}Targets.cmake")
 
-  find_package(Eigen3 3.3 CONFIG REQUIRED)
+  find_package(Eigen3 CONFIG REQUIRED)
 
   get_target_property(MRCPP_HAS_OMP MRCPP::mrcpp MRCPP_HAS_OMP)
   if(MRCPP_HAS_OMP)

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:`[32]^3`
+:math:`[-16,16]^3` in three dimensions (eight root boxes of size :math:`[16]^3`
 each):
 
 .. code-block:: cpp

external/upstream/fetch_eigen3.cmake

@@ -1,4 +1,4 @@
-find_package(Eigen3 3.4 CONFIG QUIET
+find_package(Eigen3 CONFIG QUIET
   NO_CMAKE_PATH
   NO_CMAKE_PACKAGE_REGISTRY
   )

src/core/CrossCorrelation.cpp

@@ -106,11 +106,11 @@ void CrossCorrelation::readCCCBin() {
     int K = this->order + 1;
     this->Left = MatrixXd::Zero(K * K, 2 * K);
     this->Right = MatrixXd::Zero(K * K, 2 * K);
-    double dL[2 * K];
-    double dR[2 * K];
+    std::vector<double> dL(2 * K);
+    std::vector<double> dR(2 * K);
     for (int i = 0; i < K * K; i++) {
-        L_fis.read((char *)dL, sizeof(double) * 2 * K);
-        R_fis.read((char *)dR, sizeof(double) * 2 * K);
+        L_fis.read(reinterpret_cast<char *>(dL.data()), sizeof(double) * 2 * K);
+        R_fis.read(reinterpret_cast<char *>(dR.data()), sizeof(double) * 2 * K);
         for (int j = 0; j < 2 * K; j++) {
             if (std::abs(dL[j]) < MachinePrec) dL[j] = 0.0;
             if (std::abs(dR[j]) < MachinePrec) dR[j] = 0.0;

src/core/MWFilter.cpp

@@ -206,14 +206,14 @@ void MWFilter::generateBlocks() {
 
     int K = this->order + 1;
 
-    double dH[K];
-    double dG[K];
+    std::vector<double> dH(K);
+    std::vector<double> dG(K);
     /* read H0 and G0 from disk */
     this->G0 = Eigen::MatrixXd::Zero(K, K);
     this->H0 = Eigen::MatrixXd::Zero(K, K);
     for (int i = 0; i < K; i++) {
-        H_fis.read((char *)dH, sizeof(double) * K);
-        G_fis.read((char *)dG, sizeof(double) * K);
+        H_fis.read(reinterpret_cast<char *>(dH.data()), sizeof(double) * K);
+        G_fis.read(reinterpret_cast<char *>(dG.data()), sizeof(double) * K);
         for (int j = 0; j < K; j++) {
             this->G0(i, j) = dG[j]; // G0
             this->H0(i, j) = dH[j]; // H0

src/functions/GaussExp.cpp

@@ -41,7 +41,7 @@ namespace mrcpp {
 
 template <int D> double GaussExp<D>::defaultScreening = 10.0;
 
-template <int D> GaussExp<D>::GaussExp(int nTerms, double prec) {
+template <int D> GaussExp<D>::GaussExp(int nTerms) {
     for (int i = 0; i < nTerms; i++) { this->funcs.push_back(nullptr); }
 }
 

src/functions/GaussExp.h

@@ -35,8 +35,6 @@
 
 namespace mrcpp {
 
-#define GAUSS_EXP_PREC 1.e-10
-
 /** @class GaussExp
  *
  * @brief Gaussian expansion in D dimensions
@@ -53,7 +51,7 @@ namespace mrcpp {
 
 template <int D> class GaussExp : public RepresentableFunction<D, double> {
 public:
-    GaussExp(int nTerms = 0, double prec = GAUSS_EXP_PREC);
+    GaussExp(int nTerms = 0);
     GaussExp(const GaussExp<D> &gExp);
     GaussExp &operator=(const GaussExp<D> &gExp);
     ~GaussExp() override;

src/functions/Gaussian.cpp

@@ -194,7 +194,6 @@ template <int D> GaussExp<D> Gaussian<D>::periodify(const std::array<double, D>
         auto needed_cells_vec = std::vector<int>();
         for (auto i = 0; i < D; i++) {
             auto upper_bound = pos[i] + x_std;
-            auto lower_bound = pos[i] - x_std;
             // number of cells upp and down relative to the center of the Gaussian
             needed_cells_vec.push_back(std::ceil(upper_bound / period[i]));
         }

src/treebuilders/ABGVCalculator.cpp

@@ -48,7 +48,7 @@ ABGVCalculator::ABGVCalculator(const ScalingBasis &basis, double a, double b)
 
 void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {
     int kp1 = basis.getQuadratureOrder();
-    double sqrtCoef[kp1];
+    std::vector<double> sqrtCoef(kp1);
     for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }
 
     switch (basis.getScalingType()) {
@@ -74,7 +74,7 @@ void ABGVCalculator::calcValueVectors(const ScalingBasis &basis) {
 
 void ABGVCalculator::calcKMatrix(const ScalingBasis &basis) {
     int kp1 = basis.getQuadratureOrder();
-    double sqrtCoef[kp1];
+    std::vector<double> sqrtCoef(kp1);
     for (int i = 0; i < kp1; i++) { sqrtCoef[i] = std::sqrt(2.0 * i + 1.0); }
     getQuadratureCache(qCache);
     const VectorXd &roots = qCache.getRoots(kp1);
@@ -105,7 +105,6 @@ void ABGVCalculator::calcNode(MWNode<2> &node) {
     node.zeroCoefs();
 
     const auto &idx = node.getNodeIndex();
-    int l = idx[1] - idx[0];
     int np1 = idx.getScale() + 1;
     int kp1 = node.getKp1();
     int kp1_d = node.getKp1_d();

src/treebuilders/ConvolutionCalculator.cpp

@@ -143,7 +143,6 @@ template <int D, typename T> MWNodeVector<D, T> *ConvolutionCalculator<D, T>::ma
     auto *band = new MWNodeVector<D, T>;
 
     int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
-    int g_depth = gNode.getDepth();
     int width = this->oper->getMaxBandWidth(o_depth);
 
     bool periodic = gNode.getMWTree().isPeriodic();
@@ -228,8 +227,8 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::calcNode(MWNode<D
 
     int o_depth = gNode.getScale() - this->oper->getOperatorRoot();
     if (manipulateOperator and this->oper->getOperatorRoot() < 0) o_depth = gNode.getDepth();
-    T tmpCoefs[gNode.getNCoefs()];
-    OperatorState<D, T> os(gNode, tmpCoefs);
+    std::vector<T> tmpCoefs(gNode.getNCoefs());
+    OperatorState<D, T> os(gNode, tmpCoefs.data());
     this->operStat.incrementGNodeCounters(gNode);
 
     // Get all nodes in f within the bandwith of O in g

src/treebuilders/DerivativeCalculator.cpp

@@ -90,8 +90,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
     // if (this->oper->getMaxBandWidth() > 1) MSG_ABORT("Only implemented for zero bw");
     outNode.zeroCoefs();
     int nComp = (1 << D);
-    T tmpCoefs[outNode.getNCoefs()];
-    OperatorState<D, T> os(outNode, tmpCoefs);
+    std::vector<T> tmpCoefs(outNode.getNCoefs());
+    OperatorState<D, T> os(outNode, tmpCoefs.data());
 
     os.setFNode(inpNode);
     os.setFIndex(inpNode.nodeIndex);
@@ -116,8 +116,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
     gNode.zeroCoefs();
 
     int nComp = (1 << D);
-    T tmpCoefs[gNode.getNCoefs()];
-    OperatorState<D, T> os(gNode, tmpCoefs);
+    std::vector<T> tmpCoefs(gNode.getNCoefs());
+    OperatorState<D, T> os(gNode, tmpCoefs.data());
     this->operStat.incrementGNodeCounters(gNode);
 
     // Get all nodes in f within the bandwith of O in g
@@ -182,11 +182,8 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator_bw0(
     // cout<<" applyOperator "<<endl;
     MWNode<D, T> &gNode = *os.gNode;
     MWNode<D, T> &fNode = *os.fNode;
-    const NodeIndex<D> &fIdx = *os.fIdx;
-    const NodeIndex<D> &gIdx = gNode.getNodeIndex();
     int depth = gNode.getDepth();
 
-    double oNorm = 1.0;
     double **oData = os.getOperData();
 
     for (int d = 0; d < D; d++) {
@@ -218,7 +215,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper
     const NodeIndex<D> &gIdx = gNode.getNodeIndex();
     int depth = gNode.getDepth();
 
-    double oNorm = 1.0;
     double **oData = os.getOperData();
 
     for (int d = 0; d < D; d++) {
@@ -236,8 +232,6 @@ template <int D, typename T> void DerivativeCalculator<D, T>::applyOperator(Oper
 
         const OperatorNode &oNode = oTree.getNode(depth, oTransl);
         int oIdx = os.getOperIndex(d);
-        double ocn = oNode.getComponentNorm(oIdx);
-        oNorm *= ocn;
         if (this->applyDir == d) {
             oData[d] = const_cast<double *>(oNode.getCoefs()) + oIdx * os.kp1_2;
         } else {

src/treebuilders/ProjectionCalculator.cpp

@@ -25,6 +25,7 @@
 
 #include "ProjectionCalculator.h"
 #include "trees/MWNode.h"
+#include <cassert>
 
 using Eigen::MatrixXd;
 

src/treebuilders/apply.cpp

@@ -118,6 +118,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
  */
 template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -252,6 +253,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
 template <int D, typename T>
 void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -295,6 +297,7 @@ template <int D, typename T> void apply_far_field(double prec, FunctionTree<D, T
 
 template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {
 
+    out = inp.paramCopy(true);
     for (int icomp = 0; icomp < 4; icomp++) {
         if (inp.Comp[icomp] != nullptr) {
             for (int ocomp = 0; ocomp < 4; ocomp++) {
@@ -411,6 +414,7 @@ template <int D, typename T> void apply(FunctionTree<D, T> &out, DerivativeOpera
 template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir, const ComplexDouble (*metric)[4]) {
     // TODO: sums and not only each components independently, when concrete examples with non diagonal metric are tested
 
+    out = inp.paramCopy(true); // note that this will copy the factor of inp (inp.func_ptr->data.c1)
     for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
         for (int ocomp = 0; ocomp < 4; ocomp++) {
             if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
@@ -469,7 +473,7 @@ std::vector<CompFunction<3> *> gradient(DerivativeOperator<3> &oper, CompFunctio
         for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
             for (int ocomp = 0; ocomp < 4; ocomp++) {
                 if (std::norm(metric[icomp][ocomp]) > MachinePrec) {
-                    grad_d->func_ptr->Ncomp = ocomp;
+                    grad_d->func_ptr->Ncomp = ocomp + 1;
                     if (inp.isreal()) {
                         grad_d->func_ptr->isreal = 1;
                         grad_d->func_ptr->iscomplex = 0;
@@ -609,9 +613,6 @@ template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out,
 template void divergence<1, ComplexDouble>(FunctionTree<1, ComplexDouble> &out, DerivativeOperator<1> &oper, std::vector<FunctionTree<1, ComplexDouble> *> &inp);
 template void divergence<2, ComplexDouble>(FunctionTree<2, ComplexDouble> &out, DerivativeOperator<2> &oper, std::vector<FunctionTree<2, ComplexDouble> *> &inp);
 template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out, DerivativeOperator<3> &oper, std::vector<FunctionTree<3, ComplexDouble> *> &inp);
-template FunctionTreeVector<1, ComplexDouble> gradient<1>(DerivativeOperator<1> &oper, FunctionTree<1, ComplexDouble> &inp);
-template FunctionTreeVector<2, ComplexDouble> gradient<2>(DerivativeOperator<2> &oper, FunctionTree<2, ComplexDouble> &inp);
-template FunctionTreeVector<3, ComplexDouble> gradient<3>(DerivativeOperator<3> &oper, FunctionTree<3, ComplexDouble> &inp);
 
 template void apply(CompFunction<3> &out, DerivativeOperator<3> &oper, CompFunction<3> &inp, int dir = -1, const ComplexDouble (*metric)[4]);
 

src/treebuilders/apply.h

@@ -53,8 +53,7 @@ template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOpe
 template <int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> &inp);
 template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, const ComplexDouble (*metric)[4]  = defaultMetric);
 template <int D, typename T> FunctionTreeVector<D, T> gradient(DerivativeOperator<D> &oper, FunctionTree<D, T> &inp);
-// template <int D>
-std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble  (*metric)[4] = nullptr);
+std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, const ComplexDouble  (*metric)[4] = defaultMetric);
 // clang-format on
 
 } // namespace mrcpp

src/treebuilders/grid.cpp

@@ -113,7 +113,6 @@ template <int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, i
             builder.build(out, calculator, adaptor, maxIter);
         }
     } else {
-        auto period = out.getMRA().getWorldBox().getScalingFactors();
         for (auto i = 0; i < inp.size(); i++) {
             auto *gauss = inp.getFunc(i).copy();
             build_grid(out, *gauss, maxIter);

src/treebuilders/multiply.cpp

@@ -334,8 +334,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct
 
     double result = 0.0;
     int ncoef = bra.getKp1_d() * bra.getTDim();
-    T valA[ncoef];
-    T valB[ncoef];
+    std::vector<T> valA(ncoef);
+    std::vector<T> valB(ncoef);
     int nNodes = bra.getNEndNodes();
 
     for (int n = 0; n < nNodes; n++) {
@@ -345,8 +345,8 @@ template <int D, typename T> double node_norm_dot(FunctionTree<D, T> &bra, Funct
             // convert to interpolating coef, take abs, convert back
             FunctionNode<D, T> *mwNode = static_cast<FunctionNode<D, T> *>(ket.findNode(idx));
             if (mwNode == nullptr) MSG_ABORT("Trees must have same grid");
-            node.getAbsCoefs(valA);
-            mwNode->getAbsCoefs(valB);
+            node.getAbsCoefs(valA.data());
+            mwNode->getAbsCoefs(valB.data());
             for (int i = 0; i < ncoef; i++) result += std::norm(valA[i] * valB[i]);
         } else {
             // approximate by product of node norms

src/trees/FunctionNode.cpp

@@ -130,7 +130,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateInterpolating() cons
     int qOrder = this->getKp1();
     getQuadratureCache(qc);
     const VectorXd &weights = qc.getWeights(qOrder);
-    double sqWeights[qOrder];
+    std::vector<double> sqWeights(qOrder);
     for (int i = 0; i < qOrder; i++) sqWeights[i] = std::sqrt(weights[i]);
 
     int kp1_p[D];
@@ -170,7 +170,7 @@ template <int D, typename T> T FunctionNode<D, T>::integrateValues() const {
     this->getCoefs(coefs);
     int ncoefs = coefs.size();
     int ncoefChild = ncoefs / (1 << D);
-    T cc[ncoefChild];
+    std::vector<T> cc(ncoefChild);
     // factorize out the children
     for (int i = 0; i < ncoefChild; i++) cc[i] = coefs[i];
     for (int j = 1; j < (1 << D); j++)