From d82d0d16fbefc341d12fed40faeb11aa261ce609 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Thu, 20 Feb 2025 16:02:13 +0100 Subject: [PATCH 01/31] multiply with empty CompFunctions --- src/utils/CompFunction.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 174655e8..43306a34 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -432,6 +432,7 @@ template class CompFunction<3>; template void deep_copy(CompFunction *out, const CompFunction &inp) { out->func_ptr->data = inp.func_ptr->data; out->alloc(inp.Ncomp()); + if (inp.getNNodes() == 0) return; for (int i = 0; i < inp.Ncomp(); i++) { if (inp.isreal()) { inp.CompD[i]->deep_copy(out->CompD[i]); @@ -448,6 +449,7 @@ template void deep_copy(CompFunction *out, const CompFunction &inp template void deep_copy(CompFunction &out, const CompFunction &inp) { out.func_ptr->data = inp.func_ptr->data; out.alloc(inp.Ncomp()); + if (inp.getNNodes() == 0) return; for (int i = 0; i < inp.Ncomp(); i++) { if (inp.isreal()) { inp.CompD[i]->deep_copy(out.CompD[i]); @@ -564,6 +566,10 @@ template void multiply(double prec, CompFunction &out, double coef, C out.func_ptr->data = inp_a.func_ptr->data; out.func_ptr->data.shared = share; // we don't inherit the shareness out.func_ptr->conj = false; // we don't inherit conjugaison + if (inp_a.getNNodes() == 0 or inp_b.getNNodes() == 0) { + if (!out_allocated) out.alloc(out.Ncomp()); + return; + } for (int comp = 0; comp < inp_a.Ncomp(); comp++) { out.func_ptr->data.c1[comp] = inp_a.func_ptr->data.c1[comp] * inp_b.func_ptr->data.c1[comp]; // we could put this is coef if everything is real? if (inp_a.isreal() and inp_b.isreal()) { From 94bfd3b44b4d288c88057203021b71e047752dfc Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 10 Mar 2025 13:20:57 +0100 Subject: [PATCH 02/31] build before project --- src/trees/NodeAllocator.cpp | 2 +- src/utils/CompFunction.cpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp index 5079459a..a71f37c8 100644 --- a/src/trees/NodeAllocator.cpp +++ b/src/trees/NodeAllocator.cpp @@ -130,7 +130,7 @@ template int NodeAllocator::alloc(int nNodes, bool coe // we require that the index for first child is a multiple of 2**D // so that we can find the sibling rank using rank=sIdx%(2**D) - if (sIdx % nNodes != 0) MSG_ERROR(" node allocate error"); + if (sIdx % nNodes != 0) std::cout << "Warning: recommended number of siblings is 2**D" << std::endl; // fill stack status auto &status = this->stackStatus; diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 43306a34..806bf605 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -782,6 +782,7 @@ template void project(CompFunction &out, RepresentableFunctionisreal = 1; out.func_ptr->iscomplex = 0; if (out.Ncomp() < 1) out.alloc(1); + if (need_to_project) build_grid(*out.CompD[0], f); if (need_to_project) mrcpp::project(prec, *out.CompD[0], f); mpi::share_function(out, 0, 132231, mpi::comm_share); } @@ -790,6 +791,7 @@ template void project(CompFunction &out, RepresentableFunctionisreal = 0; out.func_ptr->iscomplex = 1; if (out.Ncomp() < 1) out.alloc(1); + if (need_to_project) build_grid(*out.CompC[0], f); if (need_to_project) mrcpp::project(prec, *out.CompC[0], f); mpi::share_function(out, 0, 132231, mpi::comm_share); } From ef40e4e837317ee8d15da8be883b0ce19b641fd8 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 10 Mar 2025 14:03:47 +0100 Subject: [PATCH 03/31] documentation box sizes --- docs/mrcpp_api/mwfunctions.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/mrcpp_api/mwfunctions.rst b/docs/mrcpp_api/mwfunctions.rst index 527497f5..5b2b8c6f 100644 --- a/docs/mrcpp_api/mwfunctions.rst +++ b/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 From acc1894a20e1f11683ca74f8cb1e75adb00d5a0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Mon, 10 Mar 2025 14:39:45 +0100 Subject: [PATCH 04/31] Update build-test.yml --- .github/workflows/build-test.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 21a9034e..4c501180 100644 --- a/.github/workflows/build-test.yml +++ b/.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 From e7aa2970ebcad9d364029b262f277dad6ae676a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Mon, 10 Mar 2025 14:40:28 +0100 Subject: [PATCH 05/31] Update code-coverage.yml --- .github/workflows/code-coverage.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/code-coverage.yml b/.github/workflows/code-coverage.yml index 26ea24f0..f25bb30d 100644 --- a/.github/workflows/code-coverage.yml +++ b/.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} From 1a27762c2a07cd6a5383f885346600ba6a14fd52 Mon Sep 17 00:00:00 2001 From: Peter Wind Date: Mon, 10 Mar 2025 14:52:46 +0100 Subject: [PATCH 06/31] Update src/trees/NodeAllocator.cpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Roberto Di Remigio Eikås --- src/trees/NodeAllocator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp index a71f37c8..f4f72061 100644 --- a/src/trees/NodeAllocator.cpp +++ b/src/trees/NodeAllocator.cpp @@ -130,7 +130,7 @@ template int NodeAllocator::alloc(int nNodes, bool coe // we require that the index for first child is a multiple of 2**D // so that we can find the sibling rank using rank=sIdx%(2**D) - if (sIdx % nNodes != 0) std::cout << "Warning: recommended number of siblings is 2**D" << std::endl; + if (sIdx % nNodes != 0) MSG_WARN("Warning: recommended number of siblings is 2**D"); // fill stack status auto &status = this->stackStatus; From 499d0c9e0a5b867ab23b228d9aeb9c47af5fe2be Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Wed, 12 Mar 2025 13:28:29 +0100 Subject: [PATCH 07/31] complex dot product definition --- src/utils/CompFunction.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 806bf605..46e91940 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -707,7 +707,7 @@ template void multiply(CompFunction &out, FunctionTree = int bra^\dag(r) * ket(r) dr. * * Sum of component dots. - * Notice that the ComplexDouble dot(CompFunction bra, CompFunction ket) { @@ -725,14 +725,10 @@ template ComplexDouble dot(CompFunction bra, CompFunction ket) { } else { dotprod += mrcpp::dot(*bra.CompC[comp], *ket.CompC[comp]); } - dotprod *= bra.func_ptr->data.c1[comp] * ket.func_ptr->data.c1[comp]; + dotprod *= std::conj(bra.func_ptr->data.c1[comp]) * ket.func_ptr->data.c1[comp]; dotprodtot += dotprod; } - if (bra.isreal() and ket.isreal()) { - return dotprodtot.real(); - } else { - return dotprodtot; - } + return dotprodtot; } /** @brief Compute = int |bra^\dag(r)| * |ket(r)| dr. From 0e9dc7e6a490a0f304e54190c695217b10e31c64 Mon Sep 17 00:00:00 2001 From: Peter Wind Date: Wed, 12 Mar 2025 14:14:24 +0100 Subject: [PATCH 08/31] complex dot product definition (#257) --- src/utils/CompFunction.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 806bf605..46e91940 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -707,7 +707,7 @@ template void multiply(CompFunction &out, FunctionTree = int bra^\dag(r) * ket(r) dr. * * Sum of component dots. - * Notice that the ComplexDouble dot(CompFunction bra, CompFunction ket) { @@ -725,14 +725,10 @@ template ComplexDouble dot(CompFunction bra, CompFunction ket) { } else { dotprod += mrcpp::dot(*bra.CompC[comp], *ket.CompC[comp]); } - dotprod *= bra.func_ptr->data.c1[comp] * ket.func_ptr->data.c1[comp]; + dotprod *= std::conj(bra.func_ptr->data.c1[comp]) * ket.func_ptr->data.c1[comp]; dotprodtot += dotprod; } - if (bra.isreal() and ket.isreal()) { - return dotprodtot.real(); - } else { - return dotprodtot; - } + return dotprodtot; } /** @brief Compute = int |bra^\dag(r)| * |ket(r)| dr. From ea8bd98af2b59b1e686d9f6a0d813cebe8cfda78 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Sat, 22 Mar 2025 13:41:27 +0100 Subject: [PATCH 09/31] defaultmetric for gradient for CompFunctions --- src/treebuilders/apply.cpp | 3 --- src/treebuilders/apply.h | 3 +-- src/utils/CompFunction.cpp | 2 +- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp index 77d8674b..c45c818b 100644 --- a/src/treebuilders/apply.cpp +++ b/src/treebuilders/apply.cpp @@ -609,9 +609,6 @@ template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out, template void divergence<1, ComplexDouble>(FunctionTree<1, ComplexDouble> &out, DerivativeOperator<1> &oper, std::vector *> &inp); template void divergence<2, ComplexDouble>(FunctionTree<2, ComplexDouble> &out, DerivativeOperator<2> &oper, std::vector *> &inp); template void divergence<3, ComplexDouble>(FunctionTree<3, ComplexDouble> &out, DerivativeOperator<3> &oper, std::vector *> &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]); diff --git a/src/treebuilders/apply.h b/src/treebuilders/apply.h index e50b3b87..fa5a4366 100644 --- a/src/treebuilders/apply.h +++ b/src/treebuilders/apply.h @@ -53,8 +53,7 @@ template void divergence(CompFunction &out, DerivativeOpe template void divergence(FunctionTree &out, DerivativeOperator &oper, std::vector *> &inp); template void divergence(CompFunction &out, DerivativeOperator &oper, std::vector *> *inp, const ComplexDouble (*metric)[4] = defaultMetric); template FunctionTreeVector gradient(DerivativeOperator &oper, FunctionTree &inp); -// template -std::vector*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble (*metric)[4] = nullptr); +std::vector*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, const ComplexDouble (*metric)[4] = defaultMetric); // clang-format on } // namespace mrcpp diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 46e91940..4a85bec8 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -821,7 +821,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe int N = Phi.size(); int M = Psi.size(); for (int i = 0; i < M; i++) { - for (int j; j < 4; j++) delete Psi[i].CompD[j]; + for (int j = 0; j < 4; j++) delete Psi[i].CompD[j]; Psi[i].func_ptr->isreal = 0; Psi[i].func_ptr->iscomplex = 1; } From bfd2c35a470e98077939c1b7ca3695a79e6bcd5a Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 24 Mar 2025 11:08:10 +0100 Subject: [PATCH 10/31] gradient --- src/treebuilders/apply.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp index c45c818b..40d38c4c 100644 --- a/src/treebuilders/apply.cpp +++ b/src/treebuilders/apply.cpp @@ -469,7 +469,7 @@ std::vector *> 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; From f385d0c161172aba094d52467dacc25d266f9558 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Tue, 25 Mar 2025 08:51:58 +0100 Subject: [PATCH 11/31] Derive CompFunction --- src/treebuilders/apply.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp index 40d38c4c..20764e94 100644 --- a/src/treebuilders/apply.cpp +++ b/src/treebuilders/apply.cpp @@ -411,6 +411,7 @@ template void apply(FunctionTree &out, DerivativeOpera template void apply(CompFunction &out, DerivativeOperator &oper, CompFunction &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); for (int icomp = 0; icomp < inp.Ncomp(); icomp++) { for (int ocomp = 0; ocomp < 4; ocomp++) { if (std::norm(metric[icomp][ocomp]) > MachinePrec) { From ed65f174bb0af002f18aab67bcdef81e97894a4f Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Thu, 27 Mar 2025 13:07:55 +0100 Subject: [PATCH 12/31] apply CompFunction: copy metadata of inp into out --- src/treebuilders/apply.cpp | 5 ++++- src/utils/CompFunction.cpp | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp index 20764e94..bfcbf164 100644 --- a/src/treebuilders/apply.cpp +++ b/src/treebuilders/apply.cpp @@ -118,6 +118,7 @@ template void apply(double prec, FunctionTree &out, Co */ template void apply(double prec, CompFunction &out, ConvolutionOperator &oper, const CompFunction &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 void apply(double prec, FunctionTree &out, Co template void apply(double prec, CompFunction &out, ConvolutionOperator &oper, CompFunction &inp, FunctionTreeVector *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 void apply_far_field(double prec, FunctionTree void apply_far_field(double prec, CompFunction &out, ConvolutionOperator &oper, CompFunction &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,7 +414,7 @@ template void apply(FunctionTree &out, DerivativeOpera template void apply(CompFunction &out, DerivativeOperator &oper, CompFunction &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); + 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) { diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 4a85bec8..53256389 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -122,7 +122,10 @@ template * Returns a copy without defined trees. */ CompFunction CompFunction::paramCopy(bool alloc) const { - return CompFunction(func_ptr->data, alloc); + CompFunction out(func_ptr->data, alloc); + // we do not copy tree sizes: + for (int i = 0; i < 4; i++) out.func_ptr->data.Nchunks[i] = 0; + return out; } template void CompFunction::flushMRAData() { From e92280fe06bdd132f4d40ea28b54cf5070de3f75 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Thu, 16 Oct 2025 12:55:48 +0300 Subject: [PATCH 13/31] Allow CMake to be happy with Eigen3 version 5.0.0 --- cmake/MRCPPConfig.cmake.in | 2 +- external/upstream/fetch_eigen3.cmake | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/MRCPPConfig.cmake.in b/cmake/MRCPPConfig.cmake.in index f7a06bd3..adccf088 100644 --- a/cmake/MRCPPConfig.cmake.in +++ b/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) diff --git a/external/upstream/fetch_eigen3.cmake b/external/upstream/fetch_eigen3.cmake index be8f9e0b..68a6683c 100644 --- a/external/upstream/fetch_eigen3.cmake +++ b/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 ) From 1b49cf177944b756e524732cc1cea88afc124ccd Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Thu, 16 Oct 2025 12:56:21 +0300 Subject: [PATCH 14/31] Add missing include --- src/treebuilders/ProjectionCalculator.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/treebuilders/ProjectionCalculator.cpp b/src/treebuilders/ProjectionCalculator.cpp index 93173323..15ba014d 100644 --- a/src/treebuilders/ProjectionCalculator.cpp +++ b/src/treebuilders/ProjectionCalculator.cpp @@ -25,6 +25,7 @@ #include "ProjectionCalculator.h" #include "trees/MWNode.h" +#include using Eigen::MatrixXd; From fcde7f42b6ce7eb04733c03e63ca53cca72a666d Mon Sep 17 00:00:00 2001 From: Peter Wind Date: Thu, 16 Oct 2025 15:31:08 +0200 Subject: [PATCH 15/31] bug: tree was modified if saved (#259) * bug: tree was modified if saved * Make real or imaginary parts of a complex FunctionTree * make density from a CompFunction * define an explicit real and cplx project to help the compiler --------- Co-authored-by: gitpeterwind --- src/trees/FunctionTree.cpp | 32 ++++++++++++++++++------------- src/utils/CompFunction.cpp | 39 +++++++++++++++++++++++++++++++++++++- src/utils/CompFunction.h | 3 +++ 3 files changed, 60 insertions(+), 14 deletions(-) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index a4158169..72ca68e3 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -306,6 +306,7 @@ template void FunctionTree::saveTreeTXT(const std::str int ncoefs = 1; for (int d = 0; d < D; d++) ncoefs *= kp1; int Tdim = std::pow(2, D); + T values[ncoefs * Tdim]; int nout = this->endNodeTable.size(); out << Tdim * nout << std::endl; // could output only scaling coeff? @@ -331,7 +332,9 @@ template void FunctionTree::saveTreeTXT(const std::str std::array l; NodeIndex idx = this->endNodeTable[count]->getNodeIndex(); MWNode *node = &(this->getNode(idx, false)); - T *values = node->getCoefs(); + T *coefs = node->getCoefs(); + for (int i = 0; i < ncoefs * Tdim; i++) values[i] = coefs[i]; + node->attachCoefs(values); int n = idx.getScale(); node->mwTransform(Reconstruction); node->cvTransform(Forward); @@ -350,7 +353,8 @@ template void FunctionTree::saveTreeTXT(const std::str for (int i = 0; i < ncoefs; i++) out << values[cix * ncoefs + mapMRC[i]] << " "; out << std::endl; } - } + node->attachCoefs(coefs); // put back original coeff + } out.close(); } /** @brief Write the tree structure to disk, for later use @@ -358,9 +362,9 @@ template void FunctionTree::saveTreeTXT(const std::str */ template void FunctionTree::saveTree(const std::string &file) { Timer t1; + this->deleteGenerated(); auto &allocator = this->getNodeAllocator(); - std::stringstream fname; fname << file << ".tree"; std::fstream f; @@ -370,14 +374,12 @@ template void FunctionTree::saveTree(const std::string // Write size of tree int nChunks = allocator.getNChunksUsed(); f.write((char *)&nChunks, sizeof(int)); - // Write tree data, chunk by chunk for (int iChunk = 0; iChunk < nChunks; iChunk++) { - f.write((char *)allocator.getNodeChunk(iChunk), allocator.getNodeChunkSize()); - f.write((char *)allocator.getCoefChunk(iChunk), allocator.getCoefChunkSize()); + f.write((char *)allocator.getNodeChunk(iChunk), allocator.getNodeChunkSize()); + f.write((char *)allocator.getCoefChunk(iChunk), allocator.getCoefChunkSize()); } f.close(); - this->saveTreeTXT("MRC.dat"); print::time(10, "Time write", t1); } @@ -1142,19 +1144,22 @@ template void FunctionTree::deep_copy(FunctionTree FunctionTree *FunctionTree::Real() { FunctionTree *out = new FunctionTree(this->getMRA(), this->getName()); -#pragma omp parallel num_threads(mrcpp_get_num_threads()) + out->setZero(); + // TODO: why does the omp parallelism not always work here? + //#pragma omp parallel num_threads(mrcpp_get_num_threads()) { int nNodes = this->getNEndNodes(); -#pragma omp for schedule(guided) + //#pragma omp for schedule(guided) for (int n = 0; n < nNodes; n++) { MWNode &inp_node = *this->endNodeTable[n]; - MWNode out_node = out->getNode(out_node.getNodeIndex()); // Full copy + MWNode &out_node = out->getNode(inp_node.getNodeIndex(), true); double *out_coefs = out_node.getCoefs(); const T *inp_coefs = inp_node.getCoefs(); for (int i = 0; i < inp_node.getNCoefs(); i++) { out_coefs[i] = std::real(inp_coefs[i]); } out_node.calcNorms(); } } + out->resetEndNodeTable(); out->mwTransform(BottomUp); out->calcSquareNorm(); return out; @@ -1164,13 +1169,14 @@ template FunctionTree *FunctionTree::Real() */ template FunctionTree *FunctionTree::Imag() { FunctionTree *out = new FunctionTree(this->getMRA(), this->getName()); -#pragma omp parallel num_threads(mrcpp_get_num_threads()) + out->setZero(); + //#pragma omp parallel num_threads(mrcpp_get_num_threads()) { int nNodes = this->getNEndNodes(); -#pragma omp for schedule(guided) + //#pragma omp for schedule(guided) for (int n = 0; n < nNodes; n++) { MWNode &inp_node = *this->endNodeTable[n]; - MWNode out_node = out->getNode(out_node.getNodeIndex()); // Full copy + MWNode &out_node = out->getNode(inp_node.getNodeIndex(), true); double *out_coefs = out_node.getCoefs(); const T *inp_coefs = inp_node.getCoefs(); for (int i = 0; i < inp_node.getNCoefs(); i++) { out_coefs[i] = std::imag(inp_coefs[i]); } diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 53256389..a0c9db00 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -548,6 +548,25 @@ template void linear_combination(CompFunction &out, const std::vector } } +/** @brief out = conj(inp) * inp + * + * Note that output is always real + * + */ +template void make_density(CompFunction &out, CompFunction inp, double prec) { + multiply(prec, out, 1.0, inp, inp, -1, false, false, true); + if (out.iscomplex()) { + // copy onto real components + for (int i = 0; i < out.Ncomp(); i++) { + out.CompD[i] = out.CompC[i]->Real(); + delete out.CompC[i]; + } + out.func_ptr->isreal = 1; + out.func_ptr->iscomplex = 0; + } +} + + /** @brief out = inp_a * inp_b * */ @@ -766,7 +785,15 @@ void project(CompFunction<3> &out, std::function &r)> f, d mpi::share_function(out, 0, 123123, mpi::comm_share); } -// template +void project_real(CompFunction<3> &out, std::function &r)> f, double prec) { + bool need_to_project = not(out.isShared()) or mpi::share_master(); + out.func_ptr->isreal = 1; + out.func_ptr->iscomplex = 0; + if (out.Ncomp() < 1) out.alloc(1); + if (need_to_project) mrcpp::project<3>(prec, *out.CompD[0], f); + mpi::share_function(out, 0, 123123, mpi::comm_share); +} + void project(CompFunction<3> &out, std::function &r)> f, double prec) { bool need_to_project = not(out.isShared()) or mpi::share_master(); out.func_ptr->isreal = 0; @@ -776,6 +803,15 @@ void project(CompFunction<3> &out, std::function &r mpi::share_function(out, 0, 123123, mpi::comm_share); } +void project_cplx(CompFunction<3> &out, std::function &r)> f, double prec) { + bool need_to_project = not(out.isShared()) or mpi::share_master(); + out.func_ptr->isreal = 0; + out.func_ptr->iscomplex = 1; + if (out.Ncomp() < 1) out.alloc(1); + if (need_to_project) mrcpp::project<3>(prec, *out.CompC[0], f); + mpi::share_function(out, 0, 123123, mpi::comm_share); +} + template void project(CompFunction &out, RepresentableFunction &f, double prec) { bool need_to_project = not(out.isShared()) or mpi::share_master(); out.func_ptr->isreal = 1; @@ -2614,5 +2650,6 @@ template void add(CompFunction<3> &out, ComplexDouble a, CompFunction<3> inp_a, template void linear_combination(CompFunction<3> &out, const std::vector &c, std::vector> &inp, double prec, bool conjugate); template double node_norm_dot(CompFunction<3> bra, CompFunction<3> ket); template void orthogonalize(double prec, CompFunction<3> &Bra, CompFunction<3> &Ket); +template void make_density(CompFunction<3> &out, CompFunction<3> inp, double prec); } // namespace mrcpp diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 82209819..8537aff6 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -165,6 +165,7 @@ template void linear_combination(CompFunction &out, const std::vector template void multiply(CompFunction &out, CompFunction inp_a, CompFunction inp_b, double prec, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false); template void multiply(double prec, CompFunction &out, double coef, CompFunction inp_a, CompFunction inp_b, int maxIter = -1, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false); +template void make_density(CompFunction &out, CompFunction inp, double prec); template void multiply(CompFunction &out, CompFunction inp_a, CompFunction inp_b, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false); template void multiply(CompFunction &out, CompFunction &inp_a, RepresentableFunction &f, double prec, int nrefine = 0, bool conjugate = false); template void multiply(CompFunction &out, CompFunction &inp_a, RepresentableFunction &f, double prec, int nrefine = 0, bool conjugate = false); @@ -173,7 +174,9 @@ template void multiply(CompFunction &out, FunctionTree ComplexDouble dot(CompFunction bra, CompFunction ket); template double node_norm_dot(CompFunction bra, CompFunction ket); void project(CompFunction<3> &out, std::function &r)> f, double prec); +void project_real(CompFunction<3> &out, std::function &r)> f, double prec); //overload of project is not always recognized by the compiler void project(CompFunction<3> &out, std::function &r)> f, double prec); +void project_cplx(CompFunction<3> &out, std::function &r)> f, double prec); //overload of project is not always recognized by the compiler template void project(CompFunction &out, RepresentableFunction &f, double prec); template void project(CompFunction &out, RepresentableFunction &f, double prec); template void orthogonalize(double prec, CompFunction &Bra, CompFunction &Ket); From d53366e67c9b3ba61f3bf4b38426a1b4238defbc Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 24 Oct 2025 12:40:49 +0200 Subject: [PATCH 16/31] Copy CompFunction into complex type --- src/utils/CompFunction.cpp | 20 ++++++++++++++++++++ src/utils/CompFunction.h | 1 + 2 files changed, 21 insertions(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index a0c9db00..a43ad768 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -428,6 +428,25 @@ template class CompFunction<1>; template class CompFunction<2>; template class CompFunction<3>; +/** @brief Deep copy that changes type from real to complex + * + * Deep copy: makes an exact copy with type complex from a real input + */ +template void CopyToComplex(CompFunction &out, const CompFunction &inp) { + out.func_ptr->data = inp.func_ptr->data; + out.defcomplex(); + out.alloc(inp.Ncomp()); + if (inp.getNNodes() == 0) return; + for (int i = 0; i < inp.Ncomp(); i++) { + if (inp.isreal()) { + inp.CompD[i]->CopyTreeToComplex(out.CompC[i]); + } else { + inp.CompC[i]->deep_copy(out.CompC[i]); + } + } +} + + /** @brief Deep copy * * Deep copy: meta data is copied along with the content of each component. @@ -2644,6 +2663,7 @@ template void multiply(CompFunction<3> &out, FunctionTree<3, double> &inp_a, Rep template void multiply(CompFunction<3> &out, FunctionTree<3, ComplexDouble> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate); +template void CopyToComplex(CompFunction<3> &out, const CompFunction<3> &inp); template void deep_copy(CompFunction<3> *out, const CompFunction<3> &inp); template void deep_copy(CompFunction<3> &out, const CompFunction<3> &inp); template void add(CompFunction<3> &out, ComplexDouble a, CompFunction<3> inp_a, ComplexDouble b, CompFunction<3> inp_b, double prec, bool conjugate); diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 8537aff6..2d2c5732 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -158,6 +158,7 @@ template class CompFunction { std::shared_ptr> func_ptr; }; +template void CopyToComplex(CompFunction &out, const CompFunction &inp); template void deep_copy(CompFunction *out, const CompFunction &inp); template void deep_copy(CompFunction &out, const CompFunction &inp); template void add(CompFunction &out, ComplexDouble a, CompFunction inp_a, ComplexDouble b, CompFunction inp_b, double prec, bool conjugate = false); From 954ec4db8a7f0d02d781030b46159bd3554bdd41 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 24 Oct 2025 16:18:48 +0200 Subject: [PATCH 17/31] Copy OrbitalVector into complex type --- src/utils/CompFunction.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index a43ad768..bdec0f95 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -435,6 +435,7 @@ template class CompFunction<3>; template void CopyToComplex(CompFunction &out, const CompFunction &inp) { out.func_ptr->data = inp.func_ptr->data; out.defcomplex(); + out.func_ptr->data.isreal = 0; out.alloc(inp.Ncomp()); if (inp.getNNodes() == 0) return; for (int i = 0; i < inp.Ncomp(); i++) { From 079f8ccae16038001eca24e22badc84257ec63c5 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 24 Oct 2025 12:40:49 +0200 Subject: [PATCH 18/31] Copy CompFunction into complex type --- src/utils/CompFunction.cpp | 20 ++++++++++++++++++++ src/utils/CompFunction.h | 1 + 2 files changed, 21 insertions(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index a0c9db00..a43ad768 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -428,6 +428,25 @@ template class CompFunction<1>; template class CompFunction<2>; template class CompFunction<3>; +/** @brief Deep copy that changes type from real to complex + * + * Deep copy: makes an exact copy with type complex from a real input + */ +template void CopyToComplex(CompFunction &out, const CompFunction &inp) { + out.func_ptr->data = inp.func_ptr->data; + out.defcomplex(); + out.alloc(inp.Ncomp()); + if (inp.getNNodes() == 0) return; + for (int i = 0; i < inp.Ncomp(); i++) { + if (inp.isreal()) { + inp.CompD[i]->CopyTreeToComplex(out.CompC[i]); + } else { + inp.CompC[i]->deep_copy(out.CompC[i]); + } + } +} + + /** @brief Deep copy * * Deep copy: meta data is copied along with the content of each component. @@ -2644,6 +2663,7 @@ template void multiply(CompFunction<3> &out, FunctionTree<3, double> &inp_a, Rep template void multiply(CompFunction<3> &out, FunctionTree<3, ComplexDouble> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate); +template void CopyToComplex(CompFunction<3> &out, const CompFunction<3> &inp); template void deep_copy(CompFunction<3> *out, const CompFunction<3> &inp); template void deep_copy(CompFunction<3> &out, const CompFunction<3> &inp); template void add(CompFunction<3> &out, ComplexDouble a, CompFunction<3> inp_a, ComplexDouble b, CompFunction<3> inp_b, double prec, bool conjugate); diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 8537aff6..2d2c5732 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -158,6 +158,7 @@ template class CompFunction { std::shared_ptr> func_ptr; }; +template void CopyToComplex(CompFunction &out, const CompFunction &inp); template void deep_copy(CompFunction *out, const CompFunction &inp); template void deep_copy(CompFunction &out, const CompFunction &inp); template void add(CompFunction &out, ComplexDouble a, CompFunction inp_a, ComplexDouble b, CompFunction inp_b, double prec, bool conjugate = false); From 718584a3129785ebf67bc352eede61ba974ff8a1 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 24 Oct 2025 16:18:48 +0200 Subject: [PATCH 19/31] Copy OrbitalVector into complex type --- src/utils/CompFunction.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index a43ad768..bdec0f95 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -435,6 +435,7 @@ template class CompFunction<3>; template void CopyToComplex(CompFunction &out, const CompFunction &inp) { out.func_ptr->data = inp.func_ptr->data; out.defcomplex(); + out.func_ptr->data.isreal = 0; out.alloc(inp.Ncomp()); if (inp.getNNodes() == 0) return; for (int i = 0; i < inp.Ncomp(); i++) { From 90e92bda11806ca059b2c8dfc7744c356c7480f0 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 10 Nov 2025 10:57:24 +0100 Subject: [PATCH 20/31] make real density from complex CompFunction --- src/utils/CompFunction.cpp | 35 +++++++++++++++++++++++++---------- src/utils/CompFunction.h | 6 ++++-- 2 files changed, 29 insertions(+), 12 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index bdec0f95..c604c16a 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -336,7 +336,6 @@ template const FunctionTree &CompFunction::complex( /* for backwards compatibility */ template void CompFunction::setReal(FunctionTree *tree, int i) { func_ptr->isreal = 1; - // if (CompD[i] != nullptr) delete CompD[i]; CompD[i] = tree; if (tree != nullptr) { func_ptr->Ncomp = std::max(Ncomp(), i + 1); @@ -347,7 +346,6 @@ template void CompFunction::setReal(FunctionTree *tree, in template void CompFunction::setCplx(FunctionTree *tree, int i) { func_ptr->iscomplex = 1; - // if (CompC[i] != nullptr) delete CompC[i]; CompC[i] = tree; if (tree != nullptr) { func_ptr->Ncomp = std::max(Ncomp(), i + 1); @@ -363,23 +361,34 @@ template void CompFunction::setCplx(FunctionTree *t */ template void CompFunction::add(ComplexDouble c, CompFunction inp) { - if (Ncomp() < inp.Ncomp()) { - func_ptr->data = inp.func_ptr->data; - alloc(inp.Ncomp(), true); + if (Ncomp() > 0 and Ncomp() != inp.Ncomp()) { + MSG_ABORT("Cannot add CompFunction with different number of components"<<" "<isreal() and inp.isreal() and (c.imag() < MachineZero)) { + // everything is real write as real CompD[i]->add_inplace(c.real(), *inp.CompD[i]); } else { + // write as complex if (this->isreal()) { + // change type of output! CompD[i]->CopyTreeToComplex(CompC[i]); delete CompD[i]; CompD[i] = nullptr; - func_ptr->iscomplex = true; - func_ptr->isreal = false; + func_ptr->iscomplex = 1; + func_ptr->isreal = 0; + } + if (inp.iscomplex()) { + CompC[i]->add_inplace(c, *inp.CompC[i]); + } else { + // change type of input, only temporarily + inp.CompD[i]->CopyTreeToComplex(inp.CompC[i]); + CompC[i]->add_inplace(c, *inp.CompC[i]); + delete inp.CompC[i]; + inp.CompC[i] = nullptr; } - CompC[i]->add_inplace(c, *inp.CompC[i]); } } } @@ -580,6 +589,7 @@ template void make_density(CompFunction &out, CompFunction inp, do for (int i = 0; i < out.Ncomp(); i++) { out.CompD[i] = out.CompC[i]->Real(); delete out.CompC[i]; + out.CompC[i] = nullptr; } out.func_ptr->isreal = 1; out.func_ptr->iscomplex = 0; @@ -649,6 +659,8 @@ template void multiply(double prec, CompFunction &out, double coef, C out.func_ptr->isreal = 0; delete out.CompD[comp]; delete out.CompC[comp]; + out.CompD[comp] = nullptr; + out.CompC[comp] = nullptr; if (!out_allocated) out.alloc(out.Ncomp()); build_grid(*out.CompC[comp], *inp_a.CompC[comp]); build_grid(*out.CompC[comp], *inp_b.CompC[comp]); @@ -880,7 +892,10 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe int N = Phi.size(); int M = Psi.size(); for (int i = 0; i < M; i++) { - for (int j = 0; j < 4; j++) delete Psi[i].CompD[j]; + for (int j = 0; j < 4; j++) { + delete Psi[i].CompD[j]; + Psi[i].CompD[j] = nullptr; + } Psi[i].func_ptr->isreal = 0; Psi[i].func_ptr->iscomplex = 1; } diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 2d2c5732..33011bd9 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -115,8 +115,10 @@ template class CompFunction { int conj() const { return func_ptr->data.conj; } // soft conjugate int isreal() const { return func_ptr->data.isreal; } // T=double int iscomplex() const { return func_ptr->data.iscomplex; } // T=DoubleComplex - void defreal() { func_ptr->data.isreal = 1; } // define as real - void defcomplex() { func_ptr->data.iscomplex = 1; } // define as complex + void defreal() { func_ptr->data.isreal = 1; // define as real + func_ptr->data.iscomplex = 0;} + void defcomplex() { func_ptr->data.isreal = 0; // define as complex + func_ptr->data.iscomplex = 1;} int share() const { return func_ptr->data.shared; } int *Nchunks() const { return func_ptr->data.Nchunks; } // number of chunks of each component tree From 812777476a10253f7f9148bc175f1bf5bafa1f1a Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Thu, 13 Nov 2025 08:56:58 +0100 Subject: [PATCH 21/31] remove deprecated warning --- src/trees/NodeAllocator.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp index f4f72061..cd3cf9fc 100644 --- a/src/trees/NodeAllocator.cpp +++ b/src/trees/NodeAllocator.cpp @@ -128,10 +128,6 @@ template int NodeAllocator::alloc(int nNodes, bool coe // return value is index of first new node auto sIdx = this->topStack; - // we require that the index for first child is a multiple of 2**D - // so that we can find the sibling rank using rank=sIdx%(2**D) - if (sIdx % nNodes != 0) MSG_WARN("Warning: recommended number of siblings is 2**D"); - // fill stack status auto &status = this->stackStatus; for (int i = sIdx; i < sIdx + nNodes; i++) { From e8d15598b4a5669d3e3c70030ee3396438d11e9d Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 17 Nov 2025 11:33:12 +0100 Subject: [PATCH 22/31] allow Ncomp()=0 in add --- src/utils/CompFunction.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index c604c16a..7821df7b 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -361,7 +361,7 @@ template void CompFunction::setCplx(FunctionTree *t */ template void CompFunction::add(ComplexDouble c, CompFunction inp) { - if (Ncomp() > 0 and Ncomp() != inp.Ncomp()) { + if (Ncomp() > 0 and inp.Ncomp() > 0 and Ncomp() != inp.Ncomp()) { MSG_ABORT("Cannot add CompFunction with different number of components"<<" "< Date: Tue, 18 Nov 2025 14:39:35 +0100 Subject: [PATCH 23/31] To complex and back to real density (#266) * Copy CompFunction into complex type * Copy OrbitalVector into complex type * make real density from complex CompFunction * remove deprecated warning * allow Ncomp()=0 in add --- src/trees/NodeAllocator.cpp | 4 ---- src/utils/CompFunction.cpp | 35 +++++++++++++++++++++++++---------- src/utils/CompFunction.h | 6 ++++-- 3 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp index f4f72061..cd3cf9fc 100644 --- a/src/trees/NodeAllocator.cpp +++ b/src/trees/NodeAllocator.cpp @@ -128,10 +128,6 @@ template int NodeAllocator::alloc(int nNodes, bool coe // return value is index of first new node auto sIdx = this->topStack; - // we require that the index for first child is a multiple of 2**D - // so that we can find the sibling rank using rank=sIdx%(2**D) - if (sIdx % nNodes != 0) MSG_WARN("Warning: recommended number of siblings is 2**D"); - // fill stack status auto &status = this->stackStatus; for (int i = sIdx; i < sIdx + nNodes; i++) { diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index bdec0f95..7821df7b 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -336,7 +336,6 @@ template const FunctionTree &CompFunction::complex( /* for backwards compatibility */ template void CompFunction::setReal(FunctionTree *tree, int i) { func_ptr->isreal = 1; - // if (CompD[i] != nullptr) delete CompD[i]; CompD[i] = tree; if (tree != nullptr) { func_ptr->Ncomp = std::max(Ncomp(), i + 1); @@ -347,7 +346,6 @@ template void CompFunction::setReal(FunctionTree *tree, in template void CompFunction::setCplx(FunctionTree *tree, int i) { func_ptr->iscomplex = 1; - // if (CompC[i] != nullptr) delete CompC[i]; CompC[i] = tree; if (tree != nullptr) { func_ptr->Ncomp = std::max(Ncomp(), i + 1); @@ -363,23 +361,34 @@ template void CompFunction::setCplx(FunctionTree *t */ template void CompFunction::add(ComplexDouble c, CompFunction inp) { - if (Ncomp() < inp.Ncomp()) { - func_ptr->data = inp.func_ptr->data; - alloc(inp.Ncomp(), true); + if (Ncomp() > 0 and inp.Ncomp() > 0 and Ncomp() != inp.Ncomp()) { + MSG_ABORT("Cannot add CompFunction with different number of components"<<" "<isreal() and inp.isreal() and (c.imag() < MachineZero)) { + // everything is real write as real CompD[i]->add_inplace(c.real(), *inp.CompD[i]); } else { + // write as complex if (this->isreal()) { + // change type of output! CompD[i]->CopyTreeToComplex(CompC[i]); delete CompD[i]; CompD[i] = nullptr; - func_ptr->iscomplex = true; - func_ptr->isreal = false; + func_ptr->iscomplex = 1; + func_ptr->isreal = 0; + } + if (inp.iscomplex()) { + CompC[i]->add_inplace(c, *inp.CompC[i]); + } else { + // change type of input, only temporarily + inp.CompD[i]->CopyTreeToComplex(inp.CompC[i]); + CompC[i]->add_inplace(c, *inp.CompC[i]); + delete inp.CompC[i]; + inp.CompC[i] = nullptr; } - CompC[i]->add_inplace(c, *inp.CompC[i]); } } } @@ -580,6 +589,7 @@ template void make_density(CompFunction &out, CompFunction inp, do for (int i = 0; i < out.Ncomp(); i++) { out.CompD[i] = out.CompC[i]->Real(); delete out.CompC[i]; + out.CompC[i] = nullptr; } out.func_ptr->isreal = 1; out.func_ptr->iscomplex = 0; @@ -649,6 +659,8 @@ template void multiply(double prec, CompFunction &out, double coef, C out.func_ptr->isreal = 0; delete out.CompD[comp]; delete out.CompC[comp]; + out.CompD[comp] = nullptr; + out.CompC[comp] = nullptr; if (!out_allocated) out.alloc(out.Ncomp()); build_grid(*out.CompC[comp], *inp_a.CompC[comp]); build_grid(*out.CompC[comp], *inp_b.CompC[comp]); @@ -880,7 +892,10 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe int N = Phi.size(); int M = Psi.size(); for (int i = 0; i < M; i++) { - for (int j = 0; j < 4; j++) delete Psi[i].CompD[j]; + for (int j = 0; j < 4; j++) { + delete Psi[i].CompD[j]; + Psi[i].CompD[j] = nullptr; + } Psi[i].func_ptr->isreal = 0; Psi[i].func_ptr->iscomplex = 1; } diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 2d2c5732..33011bd9 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -115,8 +115,10 @@ template class CompFunction { int conj() const { return func_ptr->data.conj; } // soft conjugate int isreal() const { return func_ptr->data.isreal; } // T=double int iscomplex() const { return func_ptr->data.iscomplex; } // T=DoubleComplex - void defreal() { func_ptr->data.isreal = 1; } // define as real - void defcomplex() { func_ptr->data.iscomplex = 1; } // define as complex + void defreal() { func_ptr->data.isreal = 1; // define as real + func_ptr->data.iscomplex = 0;} + void defcomplex() { func_ptr->data.isreal = 0; // define as complex + func_ptr->data.iscomplex = 1;} int share() const { return func_ptr->data.shared; } int *Nchunks() const { return func_ptr->data.Nchunks; } // number of chunks of each component tree From 0ce5c5649a227668416c645dc85ba6eadbba3fff Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 24 Nov 2025 09:44:41 +0100 Subject: [PATCH 24/31] add mixed types CompFunction --- src/utils/CompFunction.cpp | 92 ++++++++++++++++++++++++++++---------- 1 file changed, 68 insertions(+), 24 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 7821df7b..10422a2d 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -360,35 +360,79 @@ template void CompFunction::setCplx(FunctionTree *t * */ template void CompFunction::add(ComplexDouble c, CompFunction inp) { - - if (Ncomp() > 0 and inp.Ncomp() > 0 and Ncomp() != inp.Ncomp()) { - MSG_ABORT("Cannot add CompFunction with different number of components"<<" "<isreal() and inp.isreal() and (c.imag() < MachineZero)) { - // everything is real write as real - CompD[i]->add_inplace(c.real(), *inp.CompD[i]); - } else { - // write as complex - if (this->isreal()) { - // change type of output! + if (inp.getSquareNorm() < MachineZero) { + // nothing to add + } else if(this->getSquareNorm() < MachineZero) { + //copy + func_ptr->data = inp.func_ptr->data; + alloc(inp.Ncomp(), true); + for (int i = 0; i < inp.Ncomp(); i++) { + if(inp.isreal()) CompD[i]->add_inplace(c.real(), *inp.CompD[i]); + if(inp.iscomplex()) CompC[i]->add_inplace(c, *inp.CompC[i]); + } + } else { + // we must check if the result is real or complex + if (inp.isreal() and this->isreal() and c.imag() < MachineZero) { + if ((std::norm(func_ptr->data.c1[0]-inp.func_ptr->data.c1[0]) < MachineZero) and + c.imag() < MachineZero){ + // real (possibly with complex c1) + for (int i = 0; i < inp.Ncomp(); i++) { + CompD[i]->add_inplace(c.real(), *inp.CompD[i]); + } + }else if (func_ptr->data.c1[0].imag() > MachineZero or inp.func_ptr->data.c1[0].imag() > MachineZero){ + MSG_ABORT("Not implemented"); + } else { + if(func_ptr->data.c1[0].real() > MachineZero){ + rescale(func_ptr->data.c1[0].real()); + for (int i = 1; i < Ncomp(); i++) { + if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero) + MSG_ABORT("different scaling not implemented"); + func_ptr->data.c1[i] = {1.0, 0.0}; + } + func_ptr->data.c1[0] = {1.0, 0.0}; + } + for (int i = 0; i < inp.Ncomp(); i++) { + CompD[i]->add_inplace(c.real()*inp.func_ptr->data.c1[i].real(), *inp.CompD[i]); + } + } + } else if (inp.isreal() and this->isreal() and c.imag() > MachineZero) { + MSG_ABORT("Not implemented"); + } else if( this->iscomplex() and inp.iscomplex()) { + if (std::norm(1.0 - func_ptr->data.c1[0]) > MachineZero) { + rescale(func_ptr->data.c1[0]); + for (int i = 1; i < Ncomp(); i++) { + if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero) + MSG_ABORT("different scaling not implemented"); + func_ptr->data.c1[i] = {1.0, 0.0}; + } + func_ptr->data.c1[0] = {1.0, 0.0}; + } + for (int i = 0; i < inp.Ncomp(); i++) { + CompC[i]->add_inplace(c*inp.func_ptr->data.c1[i], *inp.CompC[i]); + } + } else if( this->isreal()) { + // we set as complex + for (int i = 0; i < Ncomp(); i++) { CompD[i]->CopyTreeToComplex(CompC[i]); delete CompD[i]; CompD[i] = nullptr; - func_ptr->iscomplex = 1; - func_ptr->isreal = 0; } - if (inp.iscomplex()) { - CompC[i]->add_inplace(c, *inp.CompC[i]); - } else { - // change type of input, only temporarily - inp.CompD[i]->CopyTreeToComplex(inp.CompC[i]); - CompC[i]->add_inplace(c, *inp.CompC[i]); - delete inp.CompC[i]; - inp.CompC[i] = nullptr; + func_ptr->iscomplex = 1; + func_ptr->isreal = 0; + if (std::norm(1.0-func_ptr->data.c1[0]) > MachineZero) { + rescale(func_ptr->data.c1[0]); + for (int i = 1; i < Ncomp(); i++) { + if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero) + MSG_ABORT("different scaling not implemented"); + func_ptr->data.c1[i] = {1.0, 0.0}; + } + func_ptr->data.c1[0] = {1.0, 0.0}; + } + for (int i = 0; i < inp.Ncomp(); i++) { + CompC[i]->add_inplace(c*inp.func_ptr->data.c1[i], *inp.CompC[i]); } + } else { + MSG_ABORT("Not implemented"); } } } From 18b8c6aa9066d465496bf65aea2dcda1e1722bae Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 24 Nov 2025 14:37:58 +0100 Subject: [PATCH 25/31] do not copy c1 factor --- src/utils/CompFunction.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 10422a2d..cdf23b44 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -125,6 +125,8 @@ CompFunction CompFunction::paramCopy(bool alloc) const { CompFunction out(func_ptr->data, alloc); // we do not copy tree sizes: for (int i = 0; i < 4; i++) out.func_ptr->data.Nchunks[i] = 0; + // we do not copy the multiplicative scalar coefficient: + for (int i = 0; i < 4; i++) out.func_ptr->data.c1[i] = 1.0; return out; } From aa020958a9917924d78b52e8ec66c9edc0d29362 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Mon, 24 Nov 2025 15:47:38 +0100 Subject: [PATCH 26/31] reverse do not copy c1! --- src/utils/CompFunction.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index cdf23b44..10422a2d 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -125,8 +125,6 @@ CompFunction CompFunction::paramCopy(bool alloc) const { CompFunction out(func_ptr->data, alloc); // we do not copy tree sizes: for (int i = 0; i < 4; i++) out.func_ptr->data.Nchunks[i] = 0; - // we do not copy the multiplicative scalar coefficient: - for (int i = 0; i < 4; i++) out.func_ptr->data.c1[i] = 1.0; return out; } From ede5dffd862b36337ba9bbcb4b45a8b22c4001a0 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Wed, 26 Nov 2025 15:16:52 +0100 Subject: [PATCH 27/31] Complex rotation --- src/utils/CompFunction.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 10422a2d..205c6336 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -2042,7 +2042,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) { if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and BraKet[orbVec[j]].func_ptr->data.n1[0] != 0) continue; - S_omp(orbVec[i], orbVec[j]) += S_temp(i, j); + S(orbVec[i], orbVec[j]) += S_temp(i, j); } } } From 6da082c5364f39c79f232ca7f52e692eae98ed76 Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 5 Dec 2025 14:32:54 +0100 Subject: [PATCH 28/31] CompFunc factor, test abs(fac) and crop --- src/utils/CompFunction.cpp | 8 ++++---- src/utils/CompFunction.h | 4 +++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 205c6336..13bd4618 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -382,7 +382,7 @@ template void CompFunction::add(ComplexDouble c, CompFunction inp) }else if (func_ptr->data.c1[0].imag() > MachineZero or inp.func_ptr->data.c1[0].imag() > MachineZero){ MSG_ABORT("Not implemented"); } else { - if(func_ptr->data.c1[0].real() > MachineZero){ + if(std::abs(1.0 - func_ptr->data.c1[0].real()) > MachineZero){ rescale(func_ptr->data.c1[0].real()); for (int i = 1; i < Ncomp(); i++) { if (std::norm(func_ptr->data.c1[0]-func_ptr->data.c1[i]) > MachineZero) @@ -437,14 +437,14 @@ template void CompFunction::add(ComplexDouble c, CompFunction inp) } } -template int CompFunction::crop(double prec) { +template int CompFunction::crop(double prec, bool absPrec) { if (prec < 0.0) return 0; int nChunksremoved = 0; for (int i = 0; i < Ncomp(); i++) { if (isreal()) { - nChunksremoved += CompD[i]->crop(prec, 1.0, false); + nChunksremoved += CompD[i]->crop(prec, 1.0, absPrec); } else { - nChunksremoved += CompC[i]->crop(prec, 1.0, false); + nChunksremoved += CompC[i]->crop(prec, 1.0, absPrec); } } return nChunksremoved; diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h index 33011bd9..3b63a63b 100644 --- a/src/utils/CompFunction.h +++ b/src/utils/CompFunction.h @@ -121,6 +121,8 @@ template class CompFunction { func_ptr->data.iscomplex = 1;} int share() const { return func_ptr->data.shared; } int *Nchunks() const { return func_ptr->data.Nchunks; } // number of chunks of each component tree + ComplexDouble getFac() const { return func_ptr->data.c1[0]; } // returns the overall multiplicative factor + void setFac(ComplexDouble fac) { func_ptr->data.c1[0] = fac; } // sets the overall multiplicative factor CompFunction paramCopy(bool alloc = false) const; ComplexDouble integrate() const; @@ -134,7 +136,7 @@ template class CompFunction { const int getRank() const { return func_ptr->rank; }; void add(ComplexDouble c, CompFunction inp); - int crop(double prec); + int crop(double prec, bool absPrec = true); void rescale(ComplexDouble c); void free(); int getSizeNodes() const; From 66a35fa3aa2a2bfd31a1ed18a19514b2fa4108cd Mon Sep 17 00:00:00 2001 From: Bin Gao Date: Thu, 11 Dec 2025 14:35:02 +0100 Subject: [PATCH 29/31] remove unused prec and GAUSS_EXP_PREC in GaussExp --- src/functions/GaussExp.cpp | 2 +- src/functions/GaussExp.h | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/functions/GaussExp.cpp b/src/functions/GaussExp.cpp index a57fe670..21432686 100644 --- a/src/functions/GaussExp.cpp +++ b/src/functions/GaussExp.cpp @@ -41,7 +41,7 @@ namespace mrcpp { template double GaussExp::defaultScreening = 10.0; -template GaussExp::GaussExp(int nTerms, double prec) { +template GaussExp::GaussExp(int nTerms) { for (int i = 0; i < nTerms; i++) { this->funcs.push_back(nullptr); } } diff --git a/src/functions/GaussExp.h b/src/functions/GaussExp.h index a4315e38..38221321 100644 --- a/src/functions/GaussExp.h +++ b/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 class GaussExp : public RepresentableFunction { public: - GaussExp(int nTerms = 0, double prec = GAUSS_EXP_PREC); + GaussExp(int nTerms = 0); GaussExp(const GaussExp &gExp); GaussExp &operator=(const GaussExp &gExp); ~GaussExp() override; From 14ebd2dc1bb984673270be68f189779dc2530f6c Mon Sep 17 00:00:00 2001 From: gitpeterwind Date: Fri, 2 Jan 2026 10:07:08 +0100 Subject: [PATCH 30/31] bug fix for omp in serial case (race condition) --- src/utils/CompFunction.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 13bd4618..306aebbf 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -1031,6 +1031,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree // 4a) make a dense contiguous matrix with the coefficient from all the orbitals using node n std::vector orbjVec; // to remember which orbital correspond to each orbVec.size(); + if (node2orbVec.count(node_ix) == 0) continue; if (node2orbVec[node_ix].size() <= 0) continue; csize = sizecoeffW; if (parindexVec_ref[n] < 0) csize = sizecoeff; // for root nodes we include scaling coeff @@ -1038,6 +1039,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe int shift = sizecoeff - sizecoeffW; // to copy only wavelet part if (parindexVec_ref[n] < 0) shift = 0; ComplexMatrix coeffBlock(csize, node2orbVec[node_ix].size()); + for (int j : node2orbVec[node_ix]) { // loop over indices of the orbitals using this node int orb_node_ix = orb2node[j][node_ix]; for (int k = 0; k < csize; k++) coeffBlock(k, orbjVec.size()) = coeffVec[j][orb_node_ix][k + shift]; @@ -1319,6 +1321,7 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree // 4a) make a dense contiguous matrix with the coefficient from all the orbitals using node n std::vector orbjVec; // to remember which orbital correspond to each orbVec.size(); + if (node2orbVec.count(node_ix) == 0) continue; if (node2orbVec[node_ix].size() <= 0) continue; csize = sizecoeffW; if (parindexVec_ref[n] < 0) csize = sizecoeff; // for root nodes we include scaling coeff @@ -1680,6 +1683,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f for (int n = 0; n < max_n; n++) { MWNode node(*(refNodes[n]), false); int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree + if (node2orbVec.count(node_ix) == 0) continue; // 3a) make values for f at this node // 3a1) get coordinates of quadrature points for this node @@ -2001,6 +2005,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) { int csize; int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree std::vector orbVec; // identifies which orbitals use this node + if (serial and node2orbVec.count(node_ix) == 0) continue; if (serial and node2orbVec[node_ix].size() <= 0) continue; if (parindexVec_ref[n] < 0) csize = sizecoeff; @@ -2149,6 +2154,7 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) { int csize; int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree std::vector orbVec; // identifies which orbitals use this node + if (serial and node2orbVec.count(node_ix) == 0) continue; if (serial and node2orbVec[node_ix].size() <= 0) continue; if (parindexVec_ref[n] < 0) csize = sizecoeff; @@ -2361,6 +2367,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect csize = sizecoeffW; if (serial) { int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree + if (node2orbVecBra.count(node_ix) == 0) continue; + if (node2orbVecKet.count(node_ix) == 0) continue; int shift = sizecoeff - sizecoeffW; // to copy only wavelet part ComplexMatrix coeffBlockBra(csize, node2orbVecBra[node_ix].size()); ComplexMatrix coeffBlockKet(csize, node2orbVecKet[node_ix].size()); @@ -2587,6 +2595,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K csize = sizecoeffW; if (serial) { int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree + if (node2orbVecBra.count(node_ix) == 0) continue; + if (node2orbVecKet.count(node_ix) == 0) continue; int shift = sizecoeff - sizecoeffW; // to copy only wavelet part DoubleMatrix coeffBlockBra(csize, node2orbVecBra[node_ix].size()); DoubleMatrix coeffBlockKet(csize, node2orbVecKet[node_ix].size()); From cec3cc96d578cdf77e8e9222cb641ab3e9ff7236 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niklas=20G=C3=B6llmann?= <153195865+msnik1999@users.noreply.github.com> Date: Tue, 6 Jan 2026 16:29:30 +0100 Subject: [PATCH 31/31] Warning removal (#274) * removed all warnings; tested with clang21 on macOS 26 on M4 processor * fixed failing tests * included suggestion; fix for compiling mrchem+mrcpp on mac --- src/core/CrossCorrelation.cpp | 8 ++++---- src/core/MWFilter.cpp | 8 ++++---- src/functions/Gaussian.cpp | 1 - src/treebuilders/ABGVCalculator.cpp | 5 ++--- src/treebuilders/ConvolutionCalculator.cpp | 5 ++--- src/treebuilders/DerivativeCalculator.cpp | 14 ++++--------- src/treebuilders/grid.cpp | 1 - src/treebuilders/multiply.cpp | 8 ++++---- src/trees/FunctionNode.cpp | 4 ++-- src/trees/FunctionTree.cpp | 24 ++++------------------ src/trees/HilbertPath.cpp | 15 ++++++++++++++ src/trees/HilbertPath.h | 8 ++++---- src/trees/MWNode.cpp | 14 ++++--------- src/trees/OperatorNode.cpp | 4 ---- src/utils/Bank.cpp | 8 ++++++-- src/utils/Bank.h | 4 +++- src/utils/CompFunction.cpp | 21 ++++++------------- src/utils/tree_utils.cpp | 22 ++++++++++---------- tests/factory_functions.h | 1 - tests/functions/polynomial.cpp | 1 - tests/operators/helmholtz_operator.cpp | 1 - tests/operators/poisson_operator.cpp | 2 -- tests/treebuilders/map.cpp | 1 - tests/trees/function_tree.cpp | 5 ++--- tests/trees/node_index.cpp | 1 - 25 files changed, 77 insertions(+), 109 deletions(-) diff --git a/src/core/CrossCorrelation.cpp b/src/core/CrossCorrelation.cpp index 8b04cf18..788491ef 100644 --- a/src/core/CrossCorrelation.cpp +++ b/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 dL(2 * K); + std::vector 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(dL.data()), sizeof(double) * 2 * K); + R_fis.read(reinterpret_cast(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; diff --git a/src/core/MWFilter.cpp b/src/core/MWFilter.cpp index 8ada4f2c..19ed535c 100644 --- a/src/core/MWFilter.cpp +++ b/src/core/MWFilter.cpp @@ -206,14 +206,14 @@ void MWFilter::generateBlocks() { int K = this->order + 1; - double dH[K]; - double dG[K]; + std::vector dH(K); + std::vector 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(dH.data()), sizeof(double) * K); + G_fis.read(reinterpret_cast(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 diff --git a/src/functions/Gaussian.cpp b/src/functions/Gaussian.cpp index 6dbfa7c5..0b8e3214 100644 --- a/src/functions/Gaussian.cpp +++ b/src/functions/Gaussian.cpp @@ -194,7 +194,6 @@ template GaussExp Gaussian::periodify(const std::array auto needed_cells_vec = std::vector(); 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])); } diff --git a/src/treebuilders/ABGVCalculator.cpp b/src/treebuilders/ABGVCalculator.cpp index 10252f46..376068b5 100644 --- a/src/treebuilders/ABGVCalculator.cpp +++ b/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 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 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(); diff --git a/src/treebuilders/ConvolutionCalculator.cpp b/src/treebuilders/ConvolutionCalculator.cpp index 497fe0dd..8cfcac1e 100644 --- a/src/treebuilders/ConvolutionCalculator.cpp +++ b/src/treebuilders/ConvolutionCalculator.cpp @@ -143,7 +143,6 @@ template MWNodeVector *ConvolutionCalculator::ma auto *band = new MWNodeVector; 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 void ConvolutionCalculator::calcNode(MWNodeoper->getOperatorRoot(); if (manipulateOperator and this->oper->getOperatorRoot() < 0) o_depth = gNode.getDepth(); - T tmpCoefs[gNode.getNCoefs()]; - OperatorState os(gNode, tmpCoefs); + std::vector tmpCoefs(gNode.getNCoefs()); + OperatorState os(gNode, tmpCoefs.data()); this->operStat.incrementGNodeCounters(gNode); // Get all nodes in f within the bandwith of O in g diff --git a/src/treebuilders/DerivativeCalculator.cpp b/src/treebuilders/DerivativeCalculator.cpp index b298d1b6..f4936c01 100644 --- a/src/treebuilders/DerivativeCalculator.cpp +++ b/src/treebuilders/DerivativeCalculator.cpp @@ -90,8 +90,8 @@ template void DerivativeCalculator::calcNode(MWNodeoper->getMaxBandWidth() > 1) MSG_ABORT("Only implemented for zero bw"); outNode.zeroCoefs(); int nComp = (1 << D); - T tmpCoefs[outNode.getNCoefs()]; - OperatorState os(outNode, tmpCoefs); + std::vector tmpCoefs(outNode.getNCoefs()); + OperatorState os(outNode, tmpCoefs.data()); os.setFNode(inpNode); os.setFIndex(inpNode.nodeIndex); @@ -116,8 +116,8 @@ template void DerivativeCalculator::calcNode(MWNode os(gNode, tmpCoefs); + std::vector tmpCoefs(gNode.getNCoefs()); + OperatorState 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 void DerivativeCalculator::applyOperator_bw0( // cout<<" applyOperator "< &gNode = *os.gNode; MWNode &fNode = *os.fNode; - const NodeIndex &fIdx = *os.fIdx; - const NodeIndex &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 void DerivativeCalculator::applyOperator(Oper const NodeIndex &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 void DerivativeCalculator::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(oNode.getCoefs()) + oIdx * os.kp1_2; } else { diff --git a/src/treebuilders/grid.cpp b/src/treebuilders/grid.cpp index 0e7fb968..eaaf3696 100644 --- a/src/treebuilders/grid.cpp +++ b/src/treebuilders/grid.cpp @@ -113,7 +113,6 @@ template void build_grid(FunctionTree &out, const GaussExp &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); diff --git a/src/treebuilders/multiply.cpp b/src/treebuilders/multiply.cpp index 4e046126..46587ffa 100644 --- a/src/treebuilders/multiply.cpp +++ b/src/treebuilders/multiply.cpp @@ -334,8 +334,8 @@ template double node_norm_dot(FunctionTree &bra, Funct double result = 0.0; int ncoef = bra.getKp1_d() * bra.getTDim(); - T valA[ncoef]; - T valB[ncoef]; + std::vector valA(ncoef); + std::vector valB(ncoef); int nNodes = bra.getNEndNodes(); for (int n = 0; n < nNodes; n++) { @@ -345,8 +345,8 @@ template double node_norm_dot(FunctionTree &bra, Funct // convert to interpolating coef, take abs, convert back FunctionNode *mwNode = static_cast *>(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 diff --git a/src/trees/FunctionNode.cpp b/src/trees/FunctionNode.cpp index ff23fb39..102e5bd2 100644 --- a/src/trees/FunctionNode.cpp +++ b/src/trees/FunctionNode.cpp @@ -130,7 +130,7 @@ template T FunctionNode::integrateInterpolating() cons int qOrder = this->getKp1(); getQuadratureCache(qc); const VectorXd &weights = qc.getWeights(qOrder); - double sqWeights[qOrder]; + std::vector sqWeights(qOrder); for (int i = 0; i < qOrder; i++) sqWeights[i] = std::sqrt(weights[i]); int kp1_p[D]; @@ -170,7 +170,7 @@ template T FunctionNode::integrateValues() const { this->getCoefs(coefs); int ncoefs = coefs.size(); int ncoefChild = ncoefs / (1 << D); - T cc[ncoefChild]; + std::vector 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++) diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp index 72ca68e3..99557843 100644 --- a/src/trees/FunctionTree.cpp +++ b/src/trees/FunctionTree.cpp @@ -290,9 +290,6 @@ template void FunctionTree::loadTreeTXT(const std::str * @param[in] file: File name */ template void FunctionTree::saveTreeTXT(const std::string &fname) { - int nRoots = this->getRootBox().size(); - MWNode **roots = this->getRootBox().getNodes(); - std::ofstream out(fname); out << std::setprecision(14); out << D << std::endl; @@ -306,7 +303,7 @@ template void FunctionTree::saveTreeTXT(const std::str int ncoefs = 1; for (int d = 0; d < D; d++) ncoefs *= kp1; int Tdim = std::pow(2, D); - T values[ncoefs * Tdim]; + std::vector values(ncoefs * Tdim); int nout = this->endNodeTable.size(); out << Tdim * nout << std::endl; // could output only scaling coeff? @@ -334,7 +331,7 @@ template void FunctionTree::saveTreeTXT(const std::str MWNode *node = &(this->getNode(idx, false)); T *coefs = node->getCoefs(); for (int i = 0; i < ncoefs * Tdim; i++) values[i] = coefs[i]; - node->attachCoefs(values); + node->attachCoefs(values.data()); int n = idx.getScale(); node->mwTransform(Reconstruction); node->cvTransform(Forward); @@ -443,7 +440,6 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct // traverse tree, and treat end nodes only std::vector *> stack; // node from this for (int i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i))); - int basis = getMRA().getScalingBasis().getScalingType(); double result = 0.0; while (stack.size() > 0) { FunctionNode<3> *Node = stack.back(); @@ -456,9 +452,9 @@ template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunct double *coefs = Node->getCoefs(); // save position of coeff, but do not use them! // The data in fmat is not organized so that two consecutive points are stored after each other in memory, so needs to copy before mwtransform, cannot use memory adress directly. int nc = fmat.cols(); - double cc[nc]; + std::vector cc(nc); for (int i = 0; i < nc; i++) cc[i] = fmat(0, i); - Node->attachCoefs(cc); + Node->attachCoefs(cc.data()); result += Node->integrateValues(); Node->attachCoefs(coefs); // put back original coeff } @@ -880,8 +876,6 @@ void FunctionTree::makeCoeffVector(std::vector &coefs, indices.clear(); parent_indices.clear(); max_index = 0; - int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d(); - int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d(); std::vector *> refstack; // nodes from refTree std::vector *> thisstack; // nodes from this Tree for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { @@ -1096,7 +1090,6 @@ template <> int FunctionTree<3, double>::saveNodesAndRmCoeff() { for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); } while (stack.size() > stack_p) { MWNode<3, double> *Node = stack[stack_p++]; - int id = 0; NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs()); for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); } } @@ -1119,7 +1112,6 @@ template <> int FunctionTree<3, ComplexDouble>::saveNodesAndRmCoeff() { for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); } while (stack.size() > stack_p) { MWNode<3, ComplexDouble> *Node = stack[stack_p++]; - int id = 0; NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs()); for (int i = 0; i < Node->getNChildren(); i++) { stack.push_back(Node->children[i]); } } @@ -1195,7 +1187,6 @@ template FunctionTree *FunctionTree::Imag() template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, ComplexDouble> *&outTree) { delete outTree; - double ref = 0.0; outTree = new FunctionTree<3, ComplexDouble>(this->getMRA()); std::vector *> instack; // node from this std::vector *> outstack; // node from outTree @@ -1204,7 +1195,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp instack.push_back(this->getRootBox().getNodes()[rIdx]); outstack.push_back(outTree->getRootBox().getNodes()[rIdx]); } - int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk()); int ncoefs = this->getNodeAllocator().getNCoefs(); while (instack.size() > 0) { // inNode and outNode are the same node in space, but on different trees @@ -1235,7 +1225,6 @@ template <> void FunctionTree<3, double>::CopyTreeToComplex(FunctionTree<3, Comp template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, ComplexDouble> *&outTree) { delete outTree; - double ref = 0.0; outTree = new FunctionTree<2, ComplexDouble>(this->getMRA()); std::vector *> instack; // node from this std::vector *> outstack; // node from outTree @@ -1244,7 +1233,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp instack.push_back(this->getRootBox().getNodes()[rIdx]); outstack.push_back(outTree->getRootBox().getNodes()[rIdx]); } - int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk()); int ncoefs = this->getNodeAllocator().getNCoefs(); while (instack.size() > 0) { // inNode and outNode are the same node in space, but on different trees @@ -1275,7 +1263,6 @@ template <> void FunctionTree<2, double>::CopyTreeToComplex(FunctionTree<2, Comp template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, ComplexDouble> *&outTree) { delete outTree; - double ref = 0.0; outTree = new FunctionTree<1, ComplexDouble>(this->getMRA()); std::vector *> instack; // node from this std::vector *> outstack; // node from outTree @@ -1284,7 +1271,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp instack.push_back(this->getRootBox().getNodes()[rIdx]); outstack.push_back(outTree->getRootBox().getNodes()[rIdx]); } - int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk()); int ncoefs = this->getNodeAllocator().getNCoefs(); while (instack.size() > 0) { // inNode and outNode are the same node in space, but on different trees @@ -1316,7 +1302,6 @@ template <> void FunctionTree<1, double>::CopyTreeToComplex(FunctionTree<1, Comp // for testing template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double> *&outTree) { delete outTree; - double ref = 0.0; // FunctionTree<3, double>* inTree = this; outTree = new FunctionTree<3, double>(this->getMRA()); std::vector *> instack; // node from this @@ -1326,7 +1311,6 @@ template <> void FunctionTree<3, double>::CopyTreeToReal(FunctionTree<3, double> instack.push_back(this->getRootBox().getNodes()[rIdx]); outstack.push_back(outTree->getRootBox().getNodes()[rIdx]); } - int nNodes = std::min(this->getNNodes(), this->getNodeAllocator().getMaxNodesPerChunk()); int ncoefs = this->getNodeAllocator().getNCoefs(); while (instack.size() > 0) { // inNode and outNode are the same node in space, but on different trees diff --git a/src/trees/HilbertPath.cpp b/src/trees/HilbertPath.cpp index b052064a..ac252985 100644 --- a/src/trees/HilbertPath.cpp +++ b/src/trees/HilbertPath.cpp @@ -114,6 +114,21 @@ const int HilbertPath<3>::hTable[12][8] = { {4,5,7,6,3,2,0,1} }; +template +short int HilbertPath::getChildPath(int hIdx) const { + return pTable[this->path][hIdx]; +} + +template +int HilbertPath::getZIndex(int hIdx) const { + return zTable[this->path][hIdx]; +} + +template +int HilbertPath::getHIndex(int zIdx) const { + return hTable[this->path][zIdx]; +} + // clang-format on template class HilbertPath<1>; diff --git a/src/trees/HilbertPath.h b/src/trees/HilbertPath.h index 519e7ac7..981c5354 100644 --- a/src/trees/HilbertPath.h +++ b/src/trees/HilbertPath.h @@ -42,10 +42,10 @@ template class HilbertPath final { } short int getPath() const { return this->path; } - short int getChildPath(int hIdx) const { return this->pTable[this->path][hIdx]; } + short int getChildPath(int hIdx) const; - int getZIndex(int hIdx) const { return this->zTable[this->path][hIdx]; } - int getHIndex(int zIdx) const { return this->hTable[this->path][zIdx]; } + int getZIndex(int hIdx) const; + int getHIndex(int zIdx) const; private: short int path{0}; @@ -54,4 +54,4 @@ template class HilbertPath final { static const int hTable[][8]; }; -} // namespace mrcpp +} // namespace mrcpp \ No newline at end of file diff --git a/src/trees/MWNode.cpp b/src/trees/MWNode.cpp index 2d521b46..9f40bed2 100644 --- a/src/trees/MWNode.cpp +++ b/src/trees/MWNode.cpp @@ -346,13 +346,9 @@ template void MWNode::giveChildrenCoefs(bool overwrite * copied/summed in the correct child node. */ template void MWNode::giveChildCoefs(int cIdx, bool overwrite) { - MWNode node_i = *this; - node_i.mwTransform(Reconstruction); - int kp1_d = this->getKp1_d(); - int nChildren = this->getTDim(); if (this->children[cIdx] == nullptr) MSG_ABORT("Child does not exist!"); MWNode &child = getMWChild(cIdx); @@ -457,8 +453,8 @@ template void MWNode::cvTransform(int operation, bool auto sb = this->getMWTree().getMRA().getScalingBasis(); const MatrixXd &S = sb.getCVMap(operation); - T o_vec[nCoefs]; - T *out_vec = o_vec; + std::vector o_vec(nCoefs); + T *out_vec = o_vec.data(); T *in_vec = this->coefs; int nChildren = this->getTDim(); @@ -566,8 +562,8 @@ template void MWNode::mwTransform(int operation) { const MWFilter &filter = getMWTree().getMRA().getFilter(); double overwrite = 0.0; - T o_vec[nCoefs]; - T *out_vec = o_vec; + std::vector o_vec(nCoefs); + T *out_vec = o_vec.data(); T *in_vec = this->coefs; for (int i = 0; i < D; i++) { @@ -1259,7 +1255,6 @@ template std::ostream &MWNode::print(std::ostream &o) * i.e. *not* same normalization as a squareNorm */ template void MWNode::setMaxSquareNorm() { - auto n = this->getScale(); this->maxWSquareNorm = calcScaledWSquareNorm(); this->maxSquareNorm = calcScaledSquareNorm(); @@ -1275,7 +1270,6 @@ template void MWNode::setMaxSquareNorm() { /** @brief recursively reset maxSquaredNorm and maxWSquareNorm of parent and descendants to value -1 */ template void MWNode::resetMaxSquareNorm() { - auto n = this->getScale(); this->maxSquareNorm = -1.0; this->maxWSquareNorm = -1.0; if (not this->isEndNode()) { diff --git a/src/trees/OperatorNode.cpp b/src/trees/OperatorNode.cpp index 37f576ea..6253ba51 100644 --- a/src/trees/OperatorNode.cpp +++ b/src/trees/OperatorNode.cpp @@ -94,10 +94,6 @@ double OperatorNode::calcComponentNorm(int i) const { * */ MatrixXd OperatorNode::getComponent(int i) { - int depth = getDepth(); - double prec = getOperTree().getNormPrecision(); - double thrs = std::max(MachinePrec, prec / (8.0 * (1 << depth))); - VectorXd coef_vec; this->getCoefs(coef_vec); diff --git a/src/utils/Bank.cpp b/src/utils/Bank.cpp index f8c111a5..94152c06 100644 --- a/src/utils/Bank.cpp +++ b/src/utils/Bank.cpp @@ -8,8 +8,10 @@ namespace mrcpp { using namespace Eigen; using namespace std; -int metadata_block[3]; // can add more metadata in future -int const size_metadata = 3; +#ifdef MRCPP_HAS_MPI + int metadata_block[3]; // can add more metadata in future + int const size_metadata = 3; +#endif Bank::~Bank() { // delete all data and accounts @@ -57,7 +59,9 @@ std::map *> get_orbid2block; // to get blo std::map mem; +#ifdef MRCPP_HAS_MPI int const MIN_SCALE = -999; // Smaller than smallest scale +#endif int naccounts = 0; void Bank::open() { diff --git a/src/utils/Bank.h b/src/utils/Bank.h index 69719c53..d0e032c3 100644 --- a/src/utils/Bank.h +++ b/src/utils/Bank.h @@ -77,7 +77,10 @@ class Bank { void clear_bank(); void remove_account(int account); // remove the content and the account + #ifdef MRCPP_HAS_MPI long long totcurrentsize = 0ll; // number of kB used by all accounts + long long maxsize = 0; // max total deposited data size (without containers) + #endif std::vector accounts; // open bank accounts std::map *> get_deposits; // gives deposits of an account std::map *> get_id2ix; @@ -85,7 +88,6 @@ class Bank { std::map *> get_queue; // gives deposits of an account std::map> *> get_readytasks; // used by task manager std::map currentsize; // total deposited data size (without containers) - long long maxsize = 0; // max total deposited data size (without containers) }; class BankAccount { diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp index 306aebbf..250796c2 100644 --- a/src/utils/CompFunction.cpp +++ b/src/utils/CompFunction.cpp @@ -1517,7 +1517,6 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, double prec) { void save_nodes(CompFunctionVector &Phi, FunctionTree<3> &refTree, BankAccount &account, int sizes) { int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d(); int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d(); - int max_nNodes = refTree.getNNodes(); std::vector coeffVec; std::vector coeffVec_cplx; std::vector scalefac; @@ -1594,7 +1593,6 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f // refine_grid(refTree, f); //to test mpi::allreduce_Tree_noCoeff(refTree, Phi, mpi::comm_wrk); - int kp1 = refTree.getKp1(); int kp1_d = refTree.getKp1_d(); int nCoefs = refTree.getTDim() * kp1_d; @@ -1688,7 +1686,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f // 3a) make values for f at this node // 3a1) get coordinates of quadrature points for this node Eigen::MatrixXd pts; // Eigen::Zero(D, nCoefs); - double fval[nCoefs]; + std::vector fval(nCoefs); Coord r; double *originalCoef = nullptr; MWNode<3> *Fnode = nullptr; @@ -1709,7 +1707,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f } else { originalCoef = Fnode->getCoefs(); for (int j = 0; j < nCoefs; j++) fval[j] = originalCoef[j]; - Fnode->attachCoefs(fval); // note that each thread has its own copy + Fnode->attachCoefs(fval.data()); // note that each thread has its own copy Fnode->mwTransform(Reconstruction); Fnode->cvTransform(Forward); } @@ -1755,8 +1753,6 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f } } else { // MPI - int count1 = 0; - int count2 = 0; TaskManager tasks(max_n); for (int nn = 0; nn < max_n; nn++) { int n = tasks.next_task(); @@ -1766,7 +1762,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f // 3a1) get coordinates of quadrature points for this node Eigen::MatrixXd pts; // Eigen::Zero(D, nCoefs); node.getExpandedChildPts(pts); // TODO: use getPrimitiveChildPts (less cache). - double fval[nCoefs]; + std::vector fval(nCoefs); Coord r; MWNode Fnode(*(refNodes[n]), false); if (Func == nullptr) { @@ -1776,17 +1772,15 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f } } else { int nIdx = Func->real().getIx(node.getNodeIndex()); - count1++; if (nIdx < 0) { // use the function f instead of Func - count2++; for (int j = 0; j < nCoefs; j++) { for (int d = 0; d < D; d++) r[d] = pts(d, j); fval[j] = f.evalf(r); } } else { - Func->real().getNodeCoeff(nIdx, fval); // fetch coef from Bank - Fnode.attachCoefs(fval); + Func->real().getNodeCoeff(nIdx, fval.data()); // fetch coef from Bank + Fnode.attachCoefs(fval.data()); Fnode.mwTransform(Reconstruction); Fnode.cvTransform(Forward); } @@ -1994,7 +1988,6 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) { } // 3) make dot product for all the nodes and accumulate into S - int ibank = 0; #pragma omp parallel if (serial) { ComplexMatrix S_omp = ComplexMatrix::Zero(N, N); // copy for each thread @@ -2143,7 +2136,6 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) { } // 3) make dot product for all the nodes and accumulate into S - int ibank = 0; #pragma omp parallel if (serial) { ComplexMatrix S_omp = ComplexMatrix::Zero(N, N); // copy for each thread @@ -2349,7 +2341,6 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect int totsiz = 0; int totget = 0; int mxtotsiz = 0; - int ibank = 0; // the omp crashes sometime for unknown reasons? #pragma omp parallel if (serial) { @@ -2578,7 +2569,6 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K int totsiz = 0; int totget = 0; int mxtotsiz = 0; - int ibank = 0; #pragma omp parallel if (serial) { DoubleMatrix S_omp = DoubleMatrix::Zero(N, M); // copy for each thread @@ -2729,6 +2719,7 @@ template ComplexDouble dot(CompFunction<3> bra, CompFunction<3> ket); template void project(CompFunction<3> &out, RepresentableFunction<3, double> &f, double prec); template void project(CompFunction<3> &out, RepresentableFunction<3, ComplexDouble> &f, double prec); template void multiply(CompFunction<3> &out, CompFunction<3> inp_a, CompFunction<3> inp_b, double prec, bool absPrec, bool useMaxNorms, bool conjugate); +template void multiply(double prec, CompFunction<3> &out, double coef, CompFunction<3> inp_a, CompFunction<3> inp_b, int maxIter = -1, bool absPrec = false, bool useMaxNorms = false, bool conjugate = false); template void multiply(CompFunction<3> &out, FunctionTree<3, double> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, FunctionTree<3, ComplexDouble> &inp_a, RepresentableFunction<3, ComplexDouble> &f, double prec, int nrefine = 0, bool conjugate); template void multiply(CompFunction<3> &out, CompFunction<3> &inp_a, RepresentableFunction<3, double> &f, double prec, int nrefine = 0, bool conjugate); diff --git a/src/utils/tree_utils.cpp b/src/utils/tree_utils.cpp index 333544f6..b01c0334 100644 --- a/src/utils/tree_utils.cpp +++ b/src/utils/tree_utils.cpp @@ -118,8 +118,8 @@ template void tree_utils::mw_transform(const MWTree &t int kp1_dm1 = math_utils::ipow(kp1, D - 1); const MWFilter &filter = tree.getMRA().getFilter(); double overwrite = 0.0; - T tmpcoeff[kp1_d * tDim]; - T tmpcoeff2[kp1_d * tDim]; + std::vector tmpcoeff(kp1_d * tDim); + std::vector tmpcoeff2(kp1_d * tDim); int ftlim = tDim; int ftlim2 = tDim; int ftlim3 = tDim; @@ -135,7 +135,7 @@ template void tree_utils::mw_transform(const MWTree &t int i = 0; int mask = 1; for (int gt = 0; gt < tDim; gt++) { - T *out = tmpcoeff + gt * kp1_d; + T *out = tmpcoeff.data() + gt * kp1_d; for (int ft = 0; ft < ftlim; ft++) { // Operate in direction i only if the bits along other // directions are identical. The bit of the direction we @@ -155,13 +155,13 @@ template void tree_utils::mw_transform(const MWTree &t i++; mask = 2; // 1 << i; for (int gt = 0; gt < tDim; gt++) { - T *out = tmpcoeff2 + gt * kp1_d; + T *out = tmpcoeff2.data() + gt * kp1_d; for (int ft = 0; ft < ftlim2; ft++) { // Operate in direction i only if the bits along other // directions are identical. The bit of the direction we // operate on determines the appropriate filter/operator if ((gt | mask) == (ft | mask)) { - T *in = tmpcoeff + ft * kp1_d; + T *in = tmpcoeff.data() + ft * kp1_d; int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1); const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation); @@ -184,7 +184,7 @@ template void tree_utils::mw_transform(const MWTree &t // directions are identical. The bit of the direction we // operate on determines the appropriate filter/operator if ((gt | mask) == (ft | mask)) { - T *in = tmpcoeff2 + ft * kp1_d; + T *in = tmpcoeff2.data() + ft * kp1_d; int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1); const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation); @@ -201,8 +201,8 @@ template void tree_utils::mw_transform(const MWTree &t if (D < 3) { T *out; - if (D == 1) out = tmpcoeff; - if (D == 2) out = tmpcoeff2; + if (D == 1) out = tmpcoeff.data(); + if (D == 2) out = tmpcoeff2.data(); if (b_overwrite) { for (int j = 0; j < tDim; j++) { for (int i = 0; i < kp1_d; i++) { coeff_out[i + j * stride] = out[i + j * kp1_d]; } @@ -234,7 +234,7 @@ template void tree_utils::mw_transform_back(MWTree<3, T> &tree, T * int kp1_dm1 = math_utils::ipow(kp1, 2); const MWFilter &filter = tree.getMRA().getFilter(); double overwrite = 0.0; - T tmpcoeff[kp1_d * tDim]; + std::vector tmpcoeff(kp1_d * tDim); int ftlim = tDim; int ftlim2 = tDim; @@ -262,7 +262,7 @@ template void tree_utils::mw_transform_back(MWTree<3, T> &tree, T * i++; mask = 2; // 1 << i; for (int gt = 0; gt < tDim; gt++) { - T *out = tmpcoeff + gt * kp1_d; + T *out = tmpcoeff.data() + gt * kp1_d; for (int ft = 0; ft < ftlim2; ft++) { // Operate in direction i only if the bits along other // directions are identical. The bit of the direction we @@ -288,7 +288,7 @@ template void tree_utils::mw_transform_back(MWTree<3, T> &tree, T * // directions are identical. The bit of the direction we // operate on determines the appropriate filter/operator if ((gt | mask) == (ft | mask)) { - T *in = tmpcoeff + ft * kp1_d; + T *in = tmpcoeff.data() + ft * kp1_d; int filter_index = 2 * ((gt >> i) & 1) + ((ft >> i) & 1); const Eigen::MatrixXd &oper = filter.getSubFilter(filter_index, operation); diff --git a/tests/factory_functions.h b/tests/factory_functions.h index d1c14a98..57112326 100644 --- a/tests/factory_functions.h +++ b/tests/factory_functions.h @@ -129,7 +129,6 @@ void initialize(mrcpp::GaussFunc **func) { double beta = 1.0e4; double alpha = std::pow(beta / mrcpp::pi, D / 2.0); double pos_data[3] = {-0.2, 0.5, 1.0}; - const auto pow = std::array{}; auto pos = mrcpp::details::convert_to_std_array(pos_data); diff --git a/tests/functions/polynomial.cpp b/tests/functions/polynomial.cpp index 48b09cc5..e643d163 100644 --- a/tests/functions/polynomial.cpp +++ b/tests/functions/polynomial.cpp @@ -223,7 +223,6 @@ SCENARIO("Polynomials can be added and multiplied", "[poly_arithmetics], [polyno Polynomial R; R = P * Q; THEN("The coefficients of R are known") { - VectorXd &cr = R.getCoefs(); REQUIRE(R.getCoefs()[0] == Catch::Approx(0.0)); REQUIRE(R.getCoefs()[1] == Catch::Approx(0.0)); REQUIRE(R.getCoefs()[2] == Catch::Approx(1.0)); diff --git a/tests/operators/helmholtz_operator.cpp b/tests/operators/helmholtz_operator.cpp index 8f570691..ce669e27 100644 --- a/tests/operators/helmholtz_operator.cpp +++ b/tests/operators/helmholtz_operator.cpp @@ -52,7 +52,6 @@ TEST_CASE("Helmholtz' kernel", "[init_helmholtz], [helmholtz_operator], [mw_oper const double exp_prec = 1.0e-4; const double proj_prec = 1.0e-3; const double ccc_prec = 1.0e-3; - const double band_prec = 1.0e-3; const int n = -3; const int k = 5; diff --git a/tests/operators/poisson_operator.cpp b/tests/operators/poisson_operator.cpp index 23bb22a0..e6ad04f6 100644 --- a/tests/operators/poisson_operator.cpp +++ b/tests/operators/poisson_operator.cpp @@ -50,7 +50,6 @@ TEST_CASE("Initialize Poisson operator", "[init_poisson], [poisson_operator], [m const double exp_prec = 1.0e-4; const double proj_prec = 1.0e-3; const double ccc_prec = 1.0e-3; - const double band_prec = 1.0e-3; const int n = -3; const int k = 5; @@ -161,7 +160,6 @@ TEST_CASE("Apply Periodic Poisson' operator", "[apply_periodic_Poisson], [poisso // 2.0*pi periodic in all dirs // UPDATE ME auto scaling_factor = std::array{pi, pi, pi}; - auto periodic = true; auto corner = std::array{-1, -1, -1}; auto boxes = std::array{2, 2, 2}; diff --git a/tests/treebuilders/map.cpp b/tests/treebuilders/map.cpp index 6be2153f..64ced1e9 100644 --- a/tests/treebuilders/map.cpp +++ b/tests/treebuilders/map.cpp @@ -80,7 +80,6 @@ template void testMapping() { const double ref_int = ref_tree.integrate(); const double ref_norm = ref_tree.getSquareNorm(); - const double inp_int = inp_tree.integrate(); const double inp_norm = inp_tree.getSquareNorm(); FMap fmap = [](double val) { return val * val; }; diff --git a/tests/trees/function_tree.cpp b/tests/trees/function_tree.cpp index d9492342..99e5e709 100644 --- a/tests/trees/function_tree.cpp +++ b/tests/trees/function_tree.cpp @@ -69,8 +69,7 @@ SCENARIO("Generating FunctionTree nodes", "[function_tree_generating], [function } template void testGeneratedNodes() { - const int depth = 3; - + unsigned int depth = 3; Coord r; if (r.size() >= 1) r[0] = -0.3; if (r.size() >= 2) r[1] = 0.6; @@ -86,7 +85,6 @@ template void testGeneratedNodes() { WHEN("a non-existing node is fetched") { MWNode &node = tree.getNode(r, depth); - THEN("there will be allocated GenNodes") { REQUIRE(tree.getNGenNodes() > 0); @@ -95,6 +93,7 @@ template void testGeneratedNodes() { THEN("there will be no GenNodes") { REQUIRE(tree.getNGenNodes() == 0); } } } + (void)&node; // Clean up the fetched node } finalize(&mra); } diff --git a/tests/trees/node_index.cpp b/tests/trees/node_index.cpp index 791c9439..78d684a6 100644 --- a/tests/trees/node_index.cpp +++ b/tests/trees/node_index.cpp @@ -72,7 +72,6 @@ template void testConstructors() { } SECTION("Parent constructor") { - int i = D; auto cIdx = nIdx->parent(); REQUIRE(cIdx.getScale() == (nIdx->getScale() - 1)); }