diff --git a/python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp b/python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp
index 8ef539902e..682b14941c 100644
--- a/python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp
+++ b/python/pyabacus/src/hsolver/py_diago_dav_subspace.hpp
@@ -130,8 +130,10 @@ class PyDiagoDavSubspace
             std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
         };
 
+        hsolver::PreOP<std::complex<double>, base_device::DEVICE_CPU, hsolver::fvec::DivTransMinusEigKernel<std::complex<double>, base_device::DEVICE_CPU>>
+            pre_op(precond_vec, hsolver::fvec::div_trans_prevec_minus_eigen<std::complex<double>>, hsolver::fval::qe_pw<double>);
         obj = std::make_unique<hsolver::Diago_DavSubspace<std::complex<double>, base_device::DEVICE_CPU>>(
-            precond_vec, 
+            pre_op.get(),
             nband, 
             nbasis, 
             dav_ndim, 
diff --git a/python/pyabacus/src/hsolver/py_diago_david.hpp b/python/pyabacus/src/hsolver/py_diago_david.hpp
index 2008eb6b85..bf2789f32d 100644
--- a/python/pyabacus/src/hsolver/py_diago_david.hpp
+++ b/python/pyabacus/src/hsolver/py_diago_david.hpp
@@ -137,8 +137,9 @@ class PyDiagoDavid
             syncmem_op()(this->ctx, this->ctx, spsi_out, psi_in, static_cast<size_t>(nbands * nrow));
         };
 
+        hsolver::PreOP<std::complex<double>> pre_op(precond_vec);
         obj = std::make_unique<hsolver::DiagoDavid<std::complex<double>, base_device::DEVICE_CPU>>(
-            precond_vec.data(), 
+            pre_op.get(),
             nband, 
             nbasis, 
             dav_ndim, 
diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp
index d11d7093f1..037740813d 100644
--- a/source/module_hsolver/diago_dav_subspace.cpp
+++ b/source/module_hsolver/diago_dav_subspace.cpp
@@ -13,7 +13,7 @@
 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&& precondition_in,
                                                 const int& nband_in,
                                                 const int& nbasis_in,
                                                 const int& david_ndim_in,
@@ -21,8 +21,8 @@ Diago_DavSubspace<T, Device>::Diago_DavSubspace(const std::vector<Real>& precond
                                                 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(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);
 
@@ -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>
@@ -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>
@@ -334,35 +319,8 @@ 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());
-        }
-    }
+    auto* start_ptr = psi_iter + nbase * this->dim;
+    this->precondition(start_ptr, start_ptr, eigenvalue_iter->data(), this->dim, notconv);
 
     // "normalize!!!" in order to improve numerical stability of subspace diagonalization
     std::vector<Real> psi_norm(notconv, 0.0);
diff --git a/source/module_hsolver/diago_dav_subspace.h b/source/module_hsolver/diago_dav_subspace.h
index d93c252602..85488af29b 100644
--- a/source/module_hsolver/diago_dav_subspace.h
+++ b/source/module_hsolver/diago_dav_subspace.h
@@ -10,6 +10,7 @@
 
 #include <vector>
 #include <functional>
+#include "module_hsolver/precondition_funcs.h"
 
 namespace hsolver
 {
@@ -23,8 +24,9 @@ class Diago_DavSubspace
     // otherwise return the real type of T(complex<float>, complex<double>)
     using Real = typename GetTypeReal<T>::type;
 
-  public:
-    Diago_DavSubspace(const std::vector<Real>& precondition_in,
+    using PreFunc = fvec::DivTransMinusEig<T>;
+public:
+    Diago_DavSubspace(PreFunc&& precondition_in,    /// pass in a function, lambda or PreOP object
                       const int& nband_in,
                       const int& nbasis_in,
                       const int& david_ndim_in,
@@ -67,9 +69,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 precondition;
+    // note that lambdas can only passed by value
 
     /// record for how many bands not have convergence eigenvalues
     int notconv = 0;
diff --git a/source/module_hsolver/diago_david.cpp b/source/module_hsolver/diago_david.cpp
index 8f404da135..f6da73b799 100644
--- a/source/module_hsolver/diago_david.cpp
+++ b/source/module_hsolver/diago_david.cpp
@@ -32,16 +32,16 @@ using namespace hsolver;
  * @note Auxiliary memory is allocated in the constructor and deallocated in the destructor.
  */
 template <typename T, typename Device>
-DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in, 
+DiagoDavid<T, Device>::DiagoDavid(PreFunc&& precondition_in,
                                   const int nband_in,
                                   const int dim_in,
                                   const int david_ndim_in,
                                   const bool use_paw_in,
                                   const diag_comm_info& diag_comm_in)
-    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in)
+    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in* nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in),
+    precondition(precondition_in)
 {
     this->device = base_device::get_device_type<Device>(this->ctx);
-    this->precondition = precondition_in;
 
     this->one = &one_;
     this->zero = &zero_;
@@ -110,15 +110,6 @@ DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in,
     // lagrange_matrix(nband, nband); // for orthogonalization
     resmem_complex_op()(this->ctx, this->lagrange_matrix, nband * nband);
     setmem_complex_op()(this->ctx, this->lagrange_matrix, 0, nband * nband);
-
-#if defined(__CUDA) || defined(__ROCM)
-    // device precondition array
-    if (this->device == base_device::GpuDevice)
-    {
-        resmem_var_op()(this->ctx, this->d_precondition, dim);
-        syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition, dim);
-    }
-#endif
 }
 
 /**
@@ -139,13 +130,6 @@ DiagoDavid<T, Device>::~DiagoDavid()
     delmem_complex_op()(this->ctx, this->vcc);
     delmem_complex_op()(this->ctx, this->lagrange_matrix);
     base_device::memory::delete_memory_op<Real, base_device::DEVICE_CPU>()(this->cpu_ctx, this->eigenvalue);
-    // If the device is a GPU device, free the d_precondition array.
-#if defined(__CUDA) || defined(__ROCM)
-    if (this->device == base_device::GpuDevice)
-    {
-        delmem_var_op()(this->ctx, this->d_precondition);
-    }
-#endif
 }
 
 template <typename T, typename Device>
@@ -499,40 +483,15 @@ void DiagoDavid<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
                               dim // LDC: if(N) max(1, m)
     );
     //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
     // Preconditioning
     // basis[nbase] = T * basis[nbase] = T * (H - lambda * S) * psi
     // where T, the preconditioner, is an approximate inverse of H
     //          T is a diagonal stored in array `precondition`
     // to do preconditioning, divide each column of basis by the corresponding element of precondition
-    for (int m = 0; m < notconv; m++)
-    {
-        //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-        if (this->device == base_device::GpuDevice)
-        {
-#if defined(__CUDA) || defined(__ROCM)
-            vector_div_vector_op<T, Device>()(this->ctx,
-                                                   dim,
-                                                   basis + dim*(nbase + m),
-                                                   basis + dim*(nbase + m),
-                                                   this->d_precondition);
-#endif
-        }
-        else
-        {
-            vector_div_vector_op<T, Device>()(this->ctx,
-                                                   dim,
-                                                   basis + dim*(nbase + m),
-                                                   basis + dim*(nbase + m),
-                                                   this->precondition);
-        }
-        //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-        // for (int ig = 0; ig < dim; ig++)
-        // {
-        //     ppsi[ig] /= this->precondition[ig];
-        // }
-        //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-    }
+    //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+    auto* start_ptr = basis + dim * nbase;
+    this->precondition(start_ptr, start_ptr, dim, notconv);
+    //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
     // there is a nbase to nbase + notconv band orthogonalise
     // plan for SchmidtOrth
diff --git a/source/module_hsolver/diago_david.h b/source/module_hsolver/diago_david.h
index 62ffc2654f..322274a760 100644
--- a/source/module_hsolver/diago_david.h
+++ b/source/module_hsolver/diago_david.h
@@ -6,7 +6,7 @@
 #include "module_base/module_device/memory_op.h"// base_device::memory
 
 #include "module_hsolver/diag_comm_info.h"
-
+#include "module_hsolver/precondition_funcs.h"
 #include <vector>
 #include <functional>
 
@@ -21,10 +21,11 @@ class DiagoDavid
     // return T if T is real type(float, double), 
     // otherwise return the real type of T(complex<float>, complex<double>)
     using Real = typename GetTypeReal<T>::type;
-  
-  public:
 
-    DiagoDavid(const Real* precondition_in,
+    using PreFunc = fvec::Div<T>;
+public:
+
+    DiagoDavid(PreFunc&& precondition_in,
                const int nband_in,
                const int dim_in,
                const int david_ndim_in,
@@ -102,8 +103,7 @@ class DiagoDavid
     int notconv = 0;
 
     /// precondition for diag, diagonal approximation of matrix A(i.e. Hamilt)
-    const Real* precondition = nullptr;
-    Real* d_precondition = nullptr;
+    const PreFunc precondition;
 
     /// eigenvalue results
     Real* eigenvalue = nullptr;
diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
index e9d64e8455..012fdb7123 100644
--- a/source/module_hsolver/hsolver_pw.cpp
+++ b/source/module_hsolver/hsolver_pw.cpp
@@ -503,7 +503,9 @@ 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,
+        // const auto pre_op = make_pre_op(pre_condition, fvec::div_trans_prevec_minus_eigen<T, Device>, fval::qe_pw<Real>);
+        const PreOP<T, Device, fvec::DivTransMinusEigKernel<T, Device>> pre_op(pre_condition, fvec::div_trans_prevec_minus_eigen<T, Device>, fval::qe_pw<Real>);
+        Diago_DavSubspace<T, Device> dav_subspace(pre_op.get(),
                                                   psi.get_nbands(),
                                                   psi.get_k_first() ? psi.get_current_nbas()
                                                                     : psi.get_nk() * psi.get_nbasis(),
@@ -572,7 +574,9 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
             ModuleBase::timer::tick("David", "spsi_func");
         };
 
-        DiagoDavid<T, Device> david(pre_condition.data(),
+        // const auto pre_op = make_pre_op(pre_condition, fvec::div_prevec<T, Device>);
+        const PreOP<T, Device> pre_op(pre_condition);
+        DiagoDavid<T, Device> david(pre_op.get(),
                                     nband,
                                     dim,
                                     PARAM.inp.pw_diag_ndim,
diff --git a/source/module_hsolver/precondition_funcs.h b/source/module_hsolver/precondition_funcs.h
new file mode 100644
index 0000000000..58834a5572
--- /dev/null
+++ b/source/module_hsolver/precondition_funcs.h
@@ -0,0 +1,146 @@
+#pragma once
+#include <vector>
+#include <functional>
+#include "module_base/module_device/types.h"
+#include "module_base/module_device/memory_op.h"
+#include "module_hsolver/kernels/math_kernel_op.h"
+
+/// @brief Preconditioner Function Library
+/// Users can add other types of operation than the following ones at one's need.
+namespace hsolver
+{
+    template <typename T>
+    using Real = typename GetTypeReal<T>::type;
+
+    /// @brief Transform a single value
+    namespace fval
+    {
+        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))); }
+    }
+
+    /// @brief Transform vectors
+    namespace fvec
+    {
+        ///---------------------------------------------------------------------------------------------
+        /// type 1: directly divide each vector by the precondition vector
+        ///---------------------------------------------------------------------------------------------
+        template <typename T, typename Device = base_device::DEVICE_CPU>
+        void div_prevec(T* const dst, const T* const src, const size_t& dim, const size_t& nvec,
+            const Real<T>* const pre)
+        {
+            Device* ctx = {};
+            for (size_t m = 0; m < nvec; m++)
+            {
+                const size_t offset = m * dim;
+                vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, pre);
+            }
+        }
+        /// Intereface to be called in the eigensolver
+        template <typename T>
+        using Div = std::function<void(T* const, const T* const, const size_t&, const size_t&)>;
+        // Kernel function full of dependence
+        template <typename T, typename Device = base_device::DEVICE_CPU>
+        using DivKernel = std::function<decltype(div_prevec<T, Device>)>;
+
+        ///---------------------------------------------------------------------------------------------
+        ///type2: Divide transval(precon_vec - eigen_subspace[m]) for each vector[m]
+        ///$X \to (A-\lambda I)^{-1} X$
+        ///---------------------------------------------------------------------------------------------
+        template <typename T, typename Device = base_device::DEVICE_CPU>
+        void div_trans_prevec_minus_eigen(T* const dst, const T* const src, 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>&)>& transval = fval::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);
+            Device* ctx = {};
+            base_device::DEVICE_CPU* cpu_ctx = {};
+            const auto device = base_device::get_device_type<Device>(ctx);
+            for (size_t m = 0; m < nvec; m++)
+            {
+                const size_t offset = m * dim;
+                for (size_t i = 0; i < dim; i++) { pre_trans[i] = transval(pre[i] - eig[m]); }
+#if defined(__CUDA) || defined(__ROCM)
+                if (device == base_device::GpuDevice)
+                {
+                    assert(d_pre);
+                    syncmem_var_h2d_op()(ctx, cpu_ctx, d_pre, pre_trans.data(), dim);
+                    vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, d_pre);
+                }
+                else
+#endif
+                {
+                    vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, pre_trans.data());
+                }
+            }
+        }
+        /// Intereface to be called in the eigensolver
+        template <typename T>
+        using DivTransMinusEig = std::function<void(T* const, const T* const, const Real<T>*, const size_t&, const size_t&)>;
+        /// Kernel function full of dependence
+        template <typename T, typename Device = base_device::DEVICE_CPU>
+        using DivTransMinusEigKernel = std::function<decltype(div_trans_prevec_minus_eigen<T, Device>)>;
+    }
+
+    /// @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 use `.get()` interface to get the function to be called by the eigensovler,
+    ///     or pass a custom lambdas/function to replace the one returned by `.get()`.
+    template <typename T, typename Device = base_device::DEVICE_CPU, typename Kernel_t = fvec::DivKernel<T, Device>>
+    struct PreOP
+    {
+        using FVal_t = std::function<Real<T>(const Real<T>&)>;  //single-value transformer
+        using resmem_real_op = base_device::memory::resize_memory_op<Real<T>, Device>;
+        using delmem_real_op = base_device::memory::delete_memory_op<Real<T>, Device>;
+        PreOP(const std::vector<Real<T>>& prevec,
+            const Kernel_t& transvec = fvec::div_prevec<T, Device>,
+            const FVal_t& transval = fval::none<Real<T>>)
+            : PreOP<T, Device, Kernel_t>(prevec.data(), prevec.size(), transvec, transval) {}
+        PreOP(const Real<T>* const prevec, const int& dim,
+            const Kernel_t& transvec = fvec::div_prevec<T, Device>,
+            const FVal_t& transval = fval::none<Real<T>>)
+            : prevec_(prevec), dim_(dim), transvec_(transvec), transval_(transval)
+        {
+            using syncmem_var_h2d_op = base_device::memory::synchronize_memory_op<Real<T>, Device, base_device::DEVICE_CPU>;
+            this->dev_ = base_device::get_device_type<Device>(this->ctx_);
+#if defined(__CUDA) || defined(__ROCM)
+            if (this->dev_ == base_device::GpuDevice) {
+                resmem_real_op()(this->ctx_, this->d_prevec_, dim_);
+                base_device::DEVICE_CPU* cpu_ctx = {};
+                syncmem_var_h2d_op()(this->ctx_, cpu_ctx, this->d_prevec_, prevec_, dim_);
+            }
+#endif
+        }
+        PreOP(const PreOP&) = delete;
+        ~PreOP() {
+#if defined(__CUDA) || defined(__ROCM)
+            if (this->dev_ == base_device::GpuDevice) { delmem_real_op()(this->ctx_, this->d_prevec_); }
+#endif
+        }
+
+        /// @brief Bind a PreOP object to a function
+        template<typename U = Kernel_t, typename std::enable_if<std::is_same<U, fvec::DivKernel<T, Device>>::value, bool>::type = 0>
+        fvec::Div<T> get() const
+        {
+            return std::bind(PreOP<T, Device, Kernel_t>::transvec_,
+                std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4,
+                this->dev_ == base_device::GpuDevice ? this->d_prevec_ : this->prevec_);
+        }
+        template<typename U = Kernel_t, typename std::enable_if<std::is_same<U, fvec::DivTransMinusEigKernel<T, Device>>::value, bool>::type = 0>
+        fvec::DivTransMinusEig<T> get() const
+        {
+            return std::bind(PreOP<T, Device, Kernel_t>::transvec_,
+                std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5,
+                this->prevec_, this->d_prevec_, this->transval_);
+        }
+
+    private:
+        const Real<T>* const prevec_ = nullptr;
+        const int dim_;
+        Real<T>* d_prevec_ = nullptr;
+        const Kernel_t transvec_;
+        const FVal_t transval_;
+        base_device::AbacusDevice_t dev_;
+        Device* ctx_ = {};
+    };
+}
diff --git a/source/module_hsolver/test/diago_david_float_test.cpp b/source/module_hsolver/test/diago_david_float_test.cpp
index da52b7272b..f1569564ca 100644
--- a/source/module_hsolver/test/diago_david_float_test.cpp
+++ b/source/module_hsolver/test/diago_david_float_test.cpp
@@ -92,8 +92,9 @@ class DiagoDavPrepare
 
 		const int dim = phi.get_current_nbas() ;
 		const int nband = phi.get_nbands();
-		const int ld_psi =phi.get_nbasis();
-		hsolver::DiagoDavid<std::complex<float>> dav(precondition, nband, dim, order, false, comm_info);
+        const int ld_psi = phi.get_nbasis();
+        const hsolver::PreOP<std::complex<float>> pre_op(precondition, dim);
+        hsolver::DiagoDavid<std::complex<float>> dav(pre_op.get(), nband, dim, order, false, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<float>>::PW_DIAG_NMAX = maxiter;
 		hsolver::DiagoIterAssist<std::complex<float>>::PW_DIAG_THR = eps;
diff --git a/source/module_hsolver/test/diago_david_real_test.cpp b/source/module_hsolver/test/diago_david_real_test.cpp
index e451d45fbe..3432fcc3bf 100644
--- a/source/module_hsolver/test/diago_david_real_test.cpp
+++ b/source/module_hsolver/test/diago_david_real_test.cpp
@@ -92,7 +92,8 @@ class DiagoDavPrepare
 		const int dim = phi.get_current_nbas();
         const int nband = phi.get_nbands();
         const int ld_psi = phi.get_nbasis();
-        hsolver::DiagoDavid<double> dav(precondition, nband, dim, order, false, comm_info);
+        const hsolver::PreOP<double> pre_op(precondition, dim);
+        hsolver::DiagoDavid<double> dav(pre_op.get(), nband, dim, order, false, comm_info);
 
         hsolver::DiagoIterAssist<double>::PW_DIAG_NMAX = maxiter;
         hsolver::DiagoIterAssist<double>::PW_DIAG_THR = eps;
diff --git a/source/module_hsolver/test/diago_david_test.cpp b/source/module_hsolver/test/diago_david_test.cpp
index 4911d40a4f..4e8f59359e 100644
--- a/source/module_hsolver/test/diago_david_test.cpp
+++ b/source/module_hsolver/test/diago_david_test.cpp
@@ -91,8 +91,10 @@ class DiagoDavPrepare
 
 		const int dim = phi.get_current_nbas();
 		const int nband = phi.get_nbands();
-		const int ld_psi = phi.get_nbasis();
-		hsolver::DiagoDavid<std::complex<double>> dav(precondition, nband, dim, order, false, comm_info);
+        const int ld_psi = phi.get_nbasis();
+        const auto pre_func = [&precondition](std::complex<double>* const dst, const std::complex<double>* const src, const int& ld, const int& nvec)->void
+            {  hsolver::fvec::div_prevec(dst, src, ld, nvec, precondition); };
+        hsolver::DiagoDavid<std::complex<double>> dav(pre_func, nband, dim, order, false, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<double>>::PW_DIAG_NMAX = maxiter;
 		hsolver::DiagoIterAssist<std::complex<double>>::PW_DIAG_THR = eps;
diff --git a/source/module_hsolver/test/hsolver_pw_sup.h b/source/module_hsolver/test/hsolver_pw_sup.h
index c70025a2c2..c2f4d98c9d 100644
--- a/source/module_hsolver/test/hsolver_pw_sup.h
+++ b/source/module_hsolver/test/hsolver_pw_sup.h
@@ -1,4 +1,5 @@
 #include "module_basis/module_pw/pw_basis_k.h"
+#include "module_hsolver/precondition_funcs.h"
 
 namespace ModulePW {
 
@@ -121,15 +122,15 @@ template class DiagoCG<std::complex<float>, base_device::DEVICE_CPU>;
 template class DiagoCG<std::complex<double>, base_device::DEVICE_CPU>;
 
 template <typename T, typename Device>
-DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in,
+DiagoDavid<T, Device>::DiagoDavid(PreFunc&& precondition_in,
                                   const int nband_in,
                                   const int dim_in,
                                   const int david_ndim_in,
                                   const bool use_paw_in,
                                   const diag_comm_info& diag_comm_in)
-    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in) {
+    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in* nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in),
+    precondition(std::forward<PreFunc>(precondition_in)) {
     this->device = base_device::get_device_type<Device>(this->ctx);
-    this->precondition = precondition_in;
 
     test_david = 2;
     // 1: check which function is called and which step is executed
diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp
index 72f690d804..349ed6da80 100644
--- a/source/module_hsolver/test/test_hsolver_sdft.cpp
+++ b/source/module_hsolver/test/test_hsolver_sdft.cpp
@@ -78,8 +78,7 @@ Stochastic_Iter<T, Device>::Stochastic_Iter()
 }
 
 template <typename T, typename Device>
-Stochastic_Iter<T, Device>::~Stochastic_Iter(){};
-template class Stochastic_Iter<std::complex<double>, base_device::DEVICE_CPU>;
+Stochastic_Iter<T, Device>::~Stochastic_Iter() {};
 
 template <typename T, typename Device>
 void Stochastic_Iter<T, Device>::init(K_Vectors* pkv_in,
@@ -151,6 +150,8 @@ void Stochastic_Iter<T, Device>::sum_stoband(Stochastic_WF<T, Device>& stowf,
     return;
 }
 
+template class Stochastic_Iter<std::complex<double>, base_device::DEVICE_CPU>;
+
 Charge::Charge(){};
 Charge::~Charge(){};
 
diff --git a/source/module_lr/hsolver_lrtd.hpp b/source/module_lr/hsolver_lrtd.hpp
index 6c8be619eb..9d9802f54d 100644
--- a/source/module_lr/hsolver_lrtd.hpp
+++ b/source/module_lr/hsolver_lrtd.hpp
@@ -76,19 +76,23 @@ namespace LR
 
                 if (method == "dav")
                 {
+                    auto pre_func = [&precondition](T* const dst, const T* const src, const int& ld, const int& nvec)->void
+                        {  hsolver::fvec::div_prevec(dst, src, ld, nvec, precondition.data()); };
                     // Allow 5 tries at most. If ntry > ntry_max = 5, exit diag loop.
                     const int ntry_max = 5;
                     // In non-self consistent calculation, do until totally converged. Else allow 5 eigenvecs to be NOT
                     // converged.
                     const int notconv_max = ("nscf" == PARAM.inp.calculation) ? 0 : 5;
                     // do diag and add davidson iteration counts up to avg_iter
-                    hsolver::DiagoDavid<T> david(precondition.data(), nband, dim, PARAM.inp.pw_diag_ndim, PARAM.inp.use_paw, comm_info);
+                    hsolver::DiagoDavid<T> david(pre_func, nband, dim, PARAM.inp.pw_diag_ndim, PARAM.inp.use_paw, comm_info);
                     hsolver::DiagoIterAssist<T>::avg_iter += static_cast<double>(david.diag(hpsi_func, spsi_func,
                         dim, psi, eigenvalue.data(), diag_ethr, maxiter, ntry_max, 0));
                 }
                 else if (method == "dav_subspace") //need refactor
                 {
-                    hsolver::Diago_DavSubspace<T> dav_subspace(precondition,
+                    auto pre_func = [&precondition](T* const dst, const T* const src, const Real<T>* eig, const int& ld, const int& nvec)->void
+                        {  hsolver::fvec::div_trans_prevec_minus_eigen(dst, src, eig, ld, nvec, precondition.data()); };
+                    hsolver::Diago_DavSubspace<T> dav_subspace(pre_func,
                         nband,
                         dim,
                         PARAM.inp.pw_diag_ndim,
diff --git a/tests/integrate/291_NO_KP_LR/result.ref b/tests/integrate/291_NO_KP_LR/result.ref
index 722b3fff61..1fd774b2fd 100644
--- a/tests/integrate/291_NO_KP_LR/result.ref
+++ b/tests/integrate/291_NO_KP_LR/result.ref
@@ -1 +1 @@
-totexcitationenergyref 0.784274
\ No newline at end of file
+totexcitationenergyref 0.786881
\ No newline at end of file