Skip to content

[pull] master from MRChemSoft:master #52

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pull merged 5 commits into from
Feb 20, 2025
Merged

[pull] master from MRChemSoft:master #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c18197d
mpi reduce for no mpi
gitpeterwind Feb 19, 2025
8e1f167
macos fix for metric. remove norm_overlap
gitpeterwind Feb 20, 2025
d88880f
Merge branch 'master' of github.com:MRChemSoft/mrcpp into merge
gitpeterwind Feb 20, 2025
cba8256
cleanup
gitpeterwind Feb 20, 2025
088503c
Merge pull request #253 from gitpeterwind/merge
gitpeterwind Feb 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 12 additions & 12 deletions src/treebuilders/apply.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
* no coefs).
*
*/
template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, ComplexDouble metric[4][4], int maxIter, bool absPrec) {
template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {

for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
for (int ocomp = 0; ocomp < 4; ocomp++) {
Expand Down Expand Up @@ -250,7 +250,7 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
}

template <int D, typename T>
void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, ComplexDouble metric[4][4], int maxIter, bool absPrec) {
void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {

for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
for (int ocomp = 0; ocomp < 4; ocomp++) {
Expand Down Expand Up @@ -293,7 +293,7 @@ template <int D, typename T> void apply_far_field(double prec, FunctionTree<D, T
apply_on_unit_cell<D>(false, prec, out, oper, inp, maxIter, absPrec);
}

template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, ComplexDouble metric[4][4], int maxIter, bool absPrec) {
template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {

for (int icomp = 0; icomp < 4; icomp++) {
if (inp.Comp[icomp] != nullptr) {
Expand Down Expand Up @@ -338,7 +338,7 @@ template <int D, typename T> void apply_near_field(double prec, FunctionTree<D,
apply_on_unit_cell<D>(true, prec, out, oper, inp, maxIter, absPrec);
}

template <int D> void apply_near_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, ComplexDouble metric[4][4], int maxIter, bool absPrec) {
template <int D> void apply_near_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4], int maxIter, bool absPrec) {

for (int icomp = 0; icomp < 4; icomp++) {
if (inp.Comp[icomp] != nullptr) {
Expand Down Expand Up @@ -408,7 +408,7 @@ template <int D, typename T> void apply(FunctionTree<D, T> &out, DerivativeOpera
print::separator(10, ' ');
}

template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir, ComplexDouble metric[4][4]) {
template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir, const ComplexDouble (*metric)[4]) {
// TODO: sums and not only each components independently, when concrete examples with non diagonal metric are tested

for (int icomp = 0; icomp < inp.Ncomp(); icomp++) {
Expand Down Expand Up @@ -461,7 +461,7 @@ template <int D, typename T> FunctionTreeVector<D, T> gradient(DerivativeOperato
return out;
}

std::vector<CompFunction<3> *> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble metric[4][4]) {
std::vector<CompFunction<3> *> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, const ComplexDouble (*metric)[4]) {
std::vector<CompFunction<3> *> out;

for (int d = 0; d < 3; d++) {
Expand Down Expand Up @@ -525,7 +525,7 @@ template <int D, typename T> void divergence(FunctionTree<D, T> &out, Derivative
clear(tmp_vec, true);
}

template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, ComplexDouble metric[4][4]) {
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, const ComplexDouble (*metric)[4]) {
MSG_ABORT("not implemented");
}

Expand All @@ -534,16 +534,16 @@ template <int D, typename T> void divergence(FunctionTree<D, T> &out, Derivative
for (auto &t : inp) inp_vec.push_back({1.0, t});
divergence(out, oper, inp_vec);
}
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, ComplexDouble metric[4][4]) {
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, const ComplexDouble (*metric)[4]) {
MSG_ABORT("not implemented");
}

template void apply<1, double>(double prec, FunctionTree<1, double> &out, ConvolutionOperator<1> &oper, FunctionTree<1, double> &inp, int maxIter, bool absPrec);
template void apply<2, double>(double prec, FunctionTree<2, double> &out, ConvolutionOperator<2> &oper, FunctionTree<2, double> &inp, int maxIter, bool absPrec);
template void apply<3, double>(double prec, FunctionTree<3, double> &out, ConvolutionOperator<3> &oper, FunctionTree<3, double> &inp, int maxIter, bool absPrec);
template void apply<1>(double prec, CompFunction<1> &out, ConvolutionOperator<1> &oper, const CompFunction<1> &inp, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template void apply<2>(double prec, CompFunction<2> &out, ConvolutionOperator<2> &oper, const CompFunction<2> &inp, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template void apply<3>(double prec, CompFunction<3> &out, ConvolutionOperator<3> &oper, const CompFunction<3> &inp, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template void apply<1>(double prec, CompFunction<1> &out, ConvolutionOperator<1> &oper, const CompFunction<1> &inp, const ComplexDouble (*metric)[4], int maxIter = -1, bool absPrec = false);
template void apply<2>(double prec, CompFunction<2> &out, ConvolutionOperator<2> &oper, const CompFunction<2> &inp, const ComplexDouble (*metric)[4], int maxIter = -1, bool absPrec = false);
template void apply<3>(double prec, CompFunction<3> &out, ConvolutionOperator<3> &oper, const CompFunction<3> &inp, const ComplexDouble (*metric)[4], int maxIter = -1, bool absPrec = false);
template void
apply<1, double>(double prec, FunctionTree<1, double> &out, ConvolutionOperator<1> &oper, FunctionTree<1, double> &inp, FunctionTreeVector<1, double> &precTrees, int maxIter, bool absPrec);
template void
Expand Down Expand Up @@ -613,6 +613,6 @@ template FunctionTreeVector<1, ComplexDouble> gradient<1>(DerivativeOperator<1>
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, ComplexDouble metric[4][4] = nullptr);
template void apply(CompFunction<3> &out, DerivativeOperator<3> &oper, CompFunction<3> &inp, int dir = -1, const ComplexDouble (*metric)[4]);

} // namespace mrcpp
18 changes: 9 additions & 9 deletions src/treebuilders/apply.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,25 +36,25 @@ template <int D, typename T> class FunctionTree;
template <int D> class DerivativeOperator;
template <int D> class ConvolutionOperator;

const ComplexDouble defaultMetric [4][4] ={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
constexpr ComplexDouble defaultMetric [4][4] ={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};

template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, ConvolutionOperator<D> &oper, FunctionTree<D, T> &inp, int maxIter = -1, bool absPrec = false);
template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, ComplexDouble metric[4][4] = defaultMetric, int maxIter = -1, bool absPrec = false);
template <int D> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, const CompFunction<D> &inp, const ComplexDouble (*metric)[4] = defaultMetric, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, ConvolutionOperator<D> &oper, FunctionTree<D, T> &inp, FunctionTreeVector<D, T> &precTrees, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, FunctionTreeVector<D, T> *precTrees, ComplexDouble (*metric)[4] = nullptr, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply_far_field(double prec, FunctionTree<D, T> &out, ConvolutionOperator<D> &oper, FunctionTree<D, T> &inp, int maxIter = -1, bool absPrec = false);
template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template <int D> void apply_far_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4] = defaultMetric, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply_near_field(double prec, FunctionTree<D, T> &out, ConvolutionOperator<D> &oper, FunctionTree<D, T> &inp, int maxIter = -1, bool absPrec = false);
template <int D> void apply_near_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, ComplexDouble metric[4][4] = nullptr, int maxIter = -1, bool absPrec = false);
template <int D> void apply_near_field(double prec, CompFunction<D> &out, ConvolutionOperator<D> &oper, CompFunction<D> &inp, const ComplexDouble (*metric)[4] = defaultMetric, int maxIter = -1, bool absPrec = false);
template <int D, typename T> void apply(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, FunctionTree<D, T> &inp, int dir = -1);
template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir = -1, ComplexDouble metric[4][4] = nullptr);
template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir = -1, const ComplexDouble (*metric)[4] = defaultMetric);
template <int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> &inp);
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, ComplexDouble metric[4][4]);
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, const ComplexDouble (*metric)[4] = defaultMetric);
template <int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> &inp);
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, ComplexDouble metric[4][4] = nullptr);
template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, const ComplexDouble (*metric)[4] = defaultMetric);
template <int D, typename T> FunctionTreeVector<D, T> gradient(DerivativeOperator<D> &oper, FunctionTree<D, T> &inp);
// template <int D>
std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble metric[4][4] = nullptr);
std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, ComplexDouble (*metric)[4] = nullptr);
// clang-format on

} // namespace mrcpp
160 changes: 0 additions & 160 deletions src/utils/CompFunction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2551,166 +2551,6 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K
return S;
}

/** @brief Compute the overlap matrix of the absolute value of the functions S_ij = <|bra_i|||ket_j|>
*
*/
DoubleMatrix calc_norm_overlap_matrix(CompFunctionVector &BraKet) {
int N = BraKet.size();
DoubleMatrix S = DoubleMatrix::Zero(N, N);
DoubleMatrix Sreal = DoubleMatrix::Zero(N, N); // same as S, but stored as 4 blocks, rr,ri,ir,ii
MultiResolutionAnalysis<3> *mra = BraKet.vecMRA;

// 1) make union tree without coefficients
mrcpp::FunctionTree<3> refTree(*mra);
mrcpp::mpi::allreduce_Tree_noCoeff(refTree, BraKet, mpi::comm_wrk);

int sizecoeff = (1 << refTree.getDim()) * refTree.getKp1_d();
int sizecoeffW = ((1 << refTree.getDim()) - 1) * refTree.getKp1_d();

// get a list of all nodes in union grid, as defined by their indices
std::vector<double> scalefac;
std::vector<double *> coeffVec_ref;
std::vector<int> indexVec_ref; // serialIx of the nodes
std::vector<int> parindexVec_ref; // serialIx of the parent nodes
int max_ix; // largest index value (not used here)

refTree.makeCoeffVector(coeffVec_ref, indexVec_ref, parindexVec_ref, scalefac, max_ix, refTree);
int max_n = indexVec_ref.size();

// only used for serial case:
std::vector<std::vector<double *>> coeffVec(N);
std::map<int, std::vector<int>> node2orbVec; // for each node index, gives a vector with the indices of the orbitals using this node
std::vector<std::map<int, int>> orb2node(N); // for a given orbital and a given node, gives the node index in
// the orbital given the node index in the reference tree

bool serial = mrcpp::mpi::wrk_size == 1; // flag for serial/MPI switch
mrcpp::BankAccount nodesBraKet;

// In the serial case we store the coeff pointers in coeffVec. In the mpi case the coeff are stored in the bank
if (serial) {
// 2) make list of all coefficients, and their reference indices
// for different orbitals, indexVec will give the same index for the same node in space
std::vector<int> parindexVec; // serialIx of the parent nodes
std::vector<int> indexVec; // serialIx of the nodes
for (int j = 0; j < N; j++) {
// make vector with all coef pointers and their indices in the union grid
if (BraKet[j].hasReal()) {
BraKet[j].real().makeCoeffVector(coeffVec[j], indexVec, parindexVec, scalefac, max_ix, refTree);
// make a map that gives j from indexVec
int orb_node_ix = 0;
for (int ix : indexVec) {
orb2node[j][ix] = orb_node_ix++;
if (ix < 0) continue;
node2orbVec[ix].push_back(j);
}
}
if (BraKet[j].hasImag()) {
BraKet[j].imag().makeCoeffVector(coeffVec[j + N], indexVec, parindexVec, scalefac, max_ix, refTree);
// make a map that gives j from indexVec
int orb_node_ix = 0;
for (int ix : indexVec) {
orb2node[j + N][ix] = orb_node_ix++;
if (ix < 0) continue;
node2orbVec[ix].push_back(j + N);
}
}
}
} else { // MPI case
// 2) send own nodes to bank, identifying them through the serialIx of refTree
save_nodes(BraKet, refTree, nodesBraKet);
mrcpp::mpi::barrier(mrcpp::mpi::comm_wrk); // wait until everything is stored before fetching!
}

// 3) make dot product for all the nodes and accumulate into S

int ibank = 0;
#pragma omp parallel for schedule(dynamic) if (serial)
for (int n = 0; n < max_n; n++) {
if (n % mrcpp::mpi::wrk_size != mrcpp::mpi::wrk_rank) continue;
int csize;
int node_ix = indexVec_ref[n]; // SerialIx for this node in the reference tree
std::vector<int> orbVec; // identifies which orbitals use this node
if (serial and node2orbVec[node_ix].size() <= 0) continue;
if (parindexVec_ref[n] < 0)
csize = sizecoeff;
else
csize = sizecoeffW;
// In the serial case we copy the coeff coeffBlock. In the mpi case coeffBlock is provided by the bank
if (serial) {
int shift = sizecoeff - sizecoeffW; // to copy only wavelet part
if (parindexVec_ref[n] < 0) shift = 0;
DoubleMatrix 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, orbVec.size()) = coeffVec[j][orb_node_ix][k + shift];
orbVec.push_back(j);
}
if (orbVec.size() > 0) {
DoubleMatrix S_temp(orbVec.size(), orbVec.size());
coeffBlock = coeffBlock.cwiseAbs();
S_temp.noalias() = coeffBlock.transpose() * coeffBlock;
for (int i = 0; i < orbVec.size(); i++) {
for (int j = 0; j < orbVec.size(); j++) {
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;
double &Srealij = Sreal(orbVec[i], orbVec[j]);
double &Stempij = S_temp(i, j);
#pragma omp atomic
Srealij += Stempij;
}
}
}
} else { // MPI case
DoubleMatrix coeffBlock(csize, N);
nodesBraKet.get_nodeblock(indexVec_ref[n], coeffBlock.data(), orbVec);

if (orbVec.size() > 0) {
DoubleMatrix S_temp(orbVec.size(), orbVec.size());
coeffBlock.conservativeResize(Eigen::NoChange, orbVec.size());
coeffBlock = coeffBlock.cwiseAbs();
S_temp.noalias() = coeffBlock.transpose() * coeffBlock;
for (int i = 0; i < orbVec.size(); i++) {
for (int j = 0; j < orbVec.size(); j++) {
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;
Sreal(orbVec[i], orbVec[j]) += S_temp(i, j);
}
}
}
}
}

IntVector conjMat = IntVector::Zero(N);
for (int i = 0; i < N; i++) {
if (!mrcpp::mpi::my_func(i)) continue;
conjMat[i] = (BraKet[i].conjugate()) ? -1 : 1;
}
mrcpp::mpi::allreduce_vector(conjMat, mrcpp::mpi::comm_wrk);

for (int i = 0; i < N; i++) {
for (int j = 0; j <= i; j++) {
S(i, j) = Sreal(i, j) + conjMat[i] * conjMat[j] * Sreal(i + N, j + N) + conjMat[j] * Sreal(i, j + N) - conjMat[i] * Sreal(i + N, j);
S(j, i) = S(i, j);
}
}

// Assumes linearity: result is sum of all nodes contributions
mrcpp::mpi::allreduce_matrix(S, mrcpp::mpi::comm_wrk);
// multiply by CompFunction multiplicative factor
ComplexVector Fac = ComplexVector::Zero(N);
for (int i = 0; i < N; i++) {
if (!mrcpp::mpi::my_func(BraKet[i])) continue;
Fac[i] = BraKet[i].func_ptr->data.c1[0];
}
mrcpp::mpi::allreduce_vector(Fac, mrcpp::mpi::comm_wrk);
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) { S(i, j) *= std::norm(std::conj(Fac[i])) * std::norm(Fac[j]); }
}
return S;
}

/** @brief Orthogonalize the functions in Bra against all orbitals in Ket
*
*/
Expand Down
Loading
Loading