Skip to content

Commit

Permalink
extract precondition function
Browse files Browse the repository at this point in the history
  • Loading branch information
maki49 committed Nov 14, 2024
1 parent a2fdb95 commit a8bb37f
Show file tree
Hide file tree
Showing 7 changed files with 113 additions and 55 deletions.
3 changes: 2 additions & 1 deletion python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,9 @@ class PyDiagoDavSubspace
std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
};

hsolver::PreOP<std::complex<double>> pre_op(precond_vec, hsolver::transfunc::qe_pw<double>);
obj = std::make_unique<hsolver::Diago_DavSubspace<std::complex<double>, base_device::DEVICE_CPU>>(
precond_vec,
hsolver::bind_pre_op(pre_op),
nband,
nbasis,
dav_ndim,
Expand Down
51 changes: 4 additions & 47 deletions source/module_hsolver/diago_dav_subspace.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,16 @@
using namespace hsolver;

template <typename T, typename Device>
Diago_DavSubspace<T, Device>::Diago_DavSubspace(const std::vector<Real>& precondition_in,
Diago_DavSubspace<T, Device>::Diago_DavSubspace(PreFunc<T>&& precondition_in,
const int& nband_in,
const int& nbasis_in,
const int& david_ndim_in,
const double& diag_thr_in,
const int& diag_nmax_in,
const bool& need_subspace_in,
const diag_comm_info& diag_comm_in)
: precondition(precondition_in), n_band(nband_in), dim(nbasis_in), nbase_x(nband_in * david_ndim_in),
diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in)
: precondition(std::forward<PreFunc<T>>(precondition_in)), n_band(nband_in), dim(nbasis_in), nbase_x(nband_in* david_ndim_in),
diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in)
{
this->device = base_device::get_device_type<Device>(this->ctx);

Expand Down Expand Up @@ -55,14 +55,6 @@ Diago_DavSubspace<T, Device>::Diago_DavSubspace(const std::vector<Real>& precond
resmem_complex_op()(this->ctx, this->vcc, this->nbase_x * this->nbase_x, "DAV::vcc");
setmem_complex_op()(this->ctx, this->vcc, 0, this->nbase_x * this->nbase_x);
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#if defined(__CUDA) || defined(__ROCM)
if (this->device == base_device::GpuDevice)
{
resmem_real_op()(this->ctx, this->d_precondition, nbasis_in);
// syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition.data(), nbasis_in);
}
#endif
}

template <typename T, typename Device>
Expand All @@ -74,13 +66,6 @@ Diago_DavSubspace<T, Device>::~Diago_DavSubspace()
delmem_complex_op()(this->ctx, this->hcc);
delmem_complex_op()(this->ctx, this->scc);
delmem_complex_op()(this->ctx, this->vcc);

#if defined(__CUDA) || defined(__ROCM)
if (this->device == base_device::GpuDevice)
{
delmem_real_op()(this->ctx, this->d_precondition);
}
#endif
}

template <typename T, typename Device>
Expand Down Expand Up @@ -334,35 +319,7 @@ void Diago_DavSubspace<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
this->dim);

// "precondition!!!"
std::vector<Real> pre(this->dim, 0.0);
for (int m = 0; m < notconv; m++)
{
for (size_t i = 0; i < this->dim; i++)
{
// pre[i] = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]);
double x = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]);
pre[i] = 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0)));
}
#if defined(__CUDA) || defined(__ROCM)
if (this->device == base_device::GpuDevice)
{
syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, pre.data(), this->dim);
vector_div_vector_op<T, Device>()(this->ctx,
this->dim,
psi_iter + (nbase + m) * this->dim,
psi_iter + (nbase + m) * this->dim,
this->d_precondition);
}
else
#endif
{
vector_div_vector_op<T, Device>()(this->ctx,
this->dim,
psi_iter + (nbase + m) * this->dim,
psi_iter + (nbase + m) * this->dim,
pre.data());
}
}
this->precondition(psi_iter + nbase * this->dim, eigenvalue_iter->data(), this->dim, notconv);

// "normalize!!!" in order to improve numerical stability of subspace diagonalization
std::vector<Real> psi_norm(notconv, 0.0);
Expand Down
9 changes: 5 additions & 4 deletions source/module_hsolver/diago_dav_subspace.h
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include <vector>
#include <functional>
#include "module_hsolver/precondition_funcs.h"

namespace hsolver
{
Expand All @@ -24,7 +25,7 @@ class Diago_DavSubspace
using Real = typename GetTypeReal<T>::type;

public:
Diago_DavSubspace(const std::vector<Real>& precondition_in,
Diago_DavSubspace(PreFunc<T>&& precondition_in, /// pass in a function, lambda or PreOP object
const int& nband_in,
const int& nbasis_in,
const int& david_ndim_in,
Expand Down Expand Up @@ -67,9 +68,9 @@ class Diago_DavSubspace
/// the maximum dimension of the reduced basis set
const int nbase_x = 0;

/// precondition for diag
const std::vector<Real>& precondition;
Real* d_precondition = nullptr;
/// The precondition operation, can be a function, lambda or PreOP object
const PreFunc<T> precondition;
// note that lambdas can only passed by value

/// record for how many bands not have convergence eigenvalues
int notconv = 0;
Expand Down
3 changes: 2 additions & 1 deletion source/module_hsolver/hsolver_pw.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -503,7 +503,8 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
};
bool scf = this->calculation_type == "nscf" ? false : true;

Diago_DavSubspace<T, Device> dav_subspace(pre_condition,
PreOP<T, Device> pre_op(pre_condition, transfunc::qe_pw<Real>);
Diago_DavSubspace<T, Device> dav_subspace(bind_pre_op(pre_op),
psi.get_nbands(),
psi.get_k_first() ? psi.get_current_nbas()
: psi.get_nk() * psi.get_nbasis(),
Expand Down
96 changes: 96 additions & 0 deletions source/module_hsolver/precondition_funcs.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
#include <functional>
#include "module_base/module_device/types.h"
#include "module_hsolver/kernels/math_kernel_op.h"
#include <iostream> // for debugging
#include <vector>
namespace hsolver
{
/// @brief Transforming a single value,
namespace transfunc
{
template <typename T> T none(const T& x) { return x; }
template <typename T> T qe_pw(const T& x) { return 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0))); }
}

template <typename T>
using Real = typename GetTypeReal<T>::type;

/// @brief to be called in the iterative eigensolver.
/// fixed parameters: object vector, eigenvalue, leading dimension, number of vectors
template <typename T>
using PreFunc = const std::function<void(T*, const Real<T>*, const size_t&, const size_t&)>;
// using PreFunc = std::function<void(T*, const Real<T>*, const int&, const int&)>;

/// type1: Divide transfunc(precon_vec - eigen_subspace[m]) for each vector[m]
///$X \to (A-\lambda I)^{-1} X$
// There may be other types of operation than this one.
template <typename T, typename Device = base_device::DEVICE_CPU>
void div_trans_prevec_minus_eigen(T* ptr, const Real<T>* eig, const size_t& dim, const size_t& nvec,
const Real<T>* const pre, Real<T>* const d_pre = nullptr, const std::function<Real<T>(const Real<T>&)>& transfunc = transfunc::none<Real<T>>)
{
using syncmem_var_h2d_op = base_device::memory::synchronize_memory_op<Real<T>, Device, base_device::DEVICE_CPU>;
std::vector<Real<T>> pre_trans(dim, 0.0);
const auto device = base_device::get_device_type<Device>({});

for (int m = 0; m < nvec; m++)
{
T* const ptr_m = ptr + m * dim;
for (size_t i = 0; i < dim; i++) { pre_trans[i] = transfunc(pre[i] - eig[m]); }
std::cout << std::endl;
#if defined(__CUDA) || defined(__ROCM)
if (device == base_device::GpuDevice)
{
assert(d_pre);
syncmem_var_h2d_op()({}, {}, d_pre, pre_trans.data(), dim);
vector_div_vector_op<T, Device>()({}, dim, ptr_m, ptr_m, d_pre);
}
else
#endif
{
vector_div_vector_op<T, Device>()({}, dim, ptr_m, ptr_m, pre_trans.data());
}
}
}

/// @brief A operator-like class of precondition function
/// to encapsulate the pre-allocation of memory on different devices before starting the iterative eigensolver.
/// One can pass the operatr() function of this class, or other custom lambdas/functions to eigensolvers.
template <typename T, typename Device = base_device::DEVICE_CPU>
struct PreOP
{
PreOP(const std::vector<Real<T>>& prevec, const std::function<Real<T>(const Real<T>&)>& transfunc = transfunc::none)
: PreOP<T, Device>(prevec.data(), prevec.size(), transfunc) {}
PreOP(const Real<T>* const prevec, const int& dim, const std::function<Real<T>(const Real<T>&)>& transfunc = transfunc::none)
: prevec_(prevec), dim_(dim), transfunc_(transfunc),
dev_(base_device::get_device_type<Device>({}))
{
#if defined(__CUDA) || defined(__ROCM)
if (this->dev_ == base_device::GpuDevice) { resmem_real_op<T, Device>()({}, this->d_prevec_, dim_); }
#endif
}
PreOP(const PreOP& other) = delete;
~PreOP() {
#if defined(__CUDA) || defined(__ROCM)
if (this->dev_ == base_device::GpuDevice) { delmem_real_op<T, Device>()({}, this->d_precondition); }
#endif
}
void operator()(T* ptr, const Real<T>* eig, const size_t& dim, const size_t& nvec) const
{
assert(dim <= dim_);
div_trans_prevec_minus_eigen<T, Device>(ptr, eig, dim, nvec, prevec_, d_prevec_, transfunc_);
}
private:
const Real<T>* const prevec_;
const int dim_;
Real<T>* d_prevec_;
const std::function<Real<T>(const Real<T>&)> transfunc_;
const base_device::AbacusDevice_t dev_;
};

/// @brief Bind a PreOP object to a function
template <typename T, typename Device>
PreFunc<T> bind_pre_op(const PreOP<T, Device>& pre_op)
{
return std::bind(&PreOP<T, Device>::operator(), &pre_op, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
}
}
4 changes: 3 additions & 1 deletion source/module_lr/hsolver_lrtd.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,8 @@ namespace LR
auto hpsi_func = [&hm](T* psi_in, T* hpsi, const int ld_psi, const int nvec) {hm.hPsi(psi_in, hpsi, ld_psi, nvec);};
auto spsi_func = [&hm](const T* psi_in, T* spsi, const int ld_psi, const int nvec)
{ std::memcpy(spsi, psi_in, sizeof(T) * ld_psi * nvec); };
auto pre_func = [&precondition](T* ptr, const Real<T>* eig, const int& ld, const int& nvec)->void
{ hsolver::div_trans_prevec_minus_eigen(ptr, eig, ld, nvec, precondition.data()); };

if (method == "dav")
{
Expand All @@ -88,7 +90,7 @@ namespace LR
}
else if (method == "dav_subspace") //need refactor
{
hsolver::Diago_DavSubspace<T> dav_subspace(precondition,
hsolver::Diago_DavSubspace<T> dav_subspace(pre_func,
nband,
dim,
PARAM.inp.pw_diag_ndim,
Expand Down
2 changes: 1 addition & 1 deletion tests/integrate/291_NO_KP_LR/result.ref
Original file line number Diff line number Diff line change
@@ -1 +1 @@
totexcitationenergyref 0.784274
totexcitationenergyref 0.786881

0 comments on commit a8bb37f

Please sign in to comment.