separate sorce and destination pointer

maki49 · maki49 · commit 06876435b10b · 2024-11-15T15:15:07.000+08:00
diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp
@@ -319,7 +319,8 @@ void Diago_DavSubspace<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
                          this->dim);
 
     // "precondition!!!"
-    this->precondition(psi_iter + nbase * this->dim, eigenvalue_iter->data(), this->dim, notconv);
+    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_david.cpp b/source/module_hsolver/diago_david.cpp
@@ -489,7 +489,8 @@ void DiagoDavid<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
     //          T is a diagonal stored in array `precondition`
     // to do preconditioning, divide each column of basis by the corresponding element of precondition
     //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-    this->precondition(basis + dim * nbase, dim, notconv);
+    auto* start_ptr = basis + dim * nbase;
+    this->precondition(start_ptr, start_ptr, dim, notconv);
     //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
     // there is a nbase to nbase + notconv band orthogonalise
diff --git a/source/module_hsolver/precondition_funcs.h b/source/module_hsolver/precondition_funcs.h
@@ -26,19 +26,19 @@ namespace hsolver
         /// type 1: directly divide each vector by the precondition vector
         ///---------------------------------------------------------------------------------------------
         template <typename T, typename Device = base_device::DEVICE_CPU>
-        void div_prevec(T* ptr, const size_t& dim, const size_t& nvec,
+        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 (int m = 0; m < nvec; m++)
+            for (size_t m = 0; m < nvec; m++)
             {
-                T* const ptr_m = ptr + m * dim;
-                vector_div_vector_op<T, Device>()(ctx, dim, ptr_m, ptr_m, pre);
+                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 size_t&, const size_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>)>;
@@ -48,7 +48,7 @@ namespace hsolver
         ///$X \to (A-\lambda I)^{-1} X$
         ///---------------------------------------------------------------------------------------------
         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,
+        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>;
@@ -57,27 +57,27 @@ namespace hsolver
             Device* ctx = {};
             base_device::DEVICE_CPU* cpu_ctx = {};
 
-            for (int m = 0; m < nvec; m++)
+            for (size_t m = 0; m < nvec; m++)
             {
-                T* const ptr_m = ptr + m * dim;
+                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, ptr_m, ptr_m, d_pre);
+                    vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, d_pre);
                 }
                 else
 #endif
                 {
-                    vector_div_vector_op<T, Device>()(ctx, dim, ptr_m, ptr_m, pre_trans.data());
+                    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 Real<T>*, const size_t&, const size_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>)>;
@@ -104,7 +104,9 @@ namespace hsolver
             dev_(base_device::get_device_type<Device>({}))
         {
 #if defined(__CUDA) || defined(__ROCM)
-            if (this->dev_ == base_device::GpuDevice) { resmem_real_op()({}, this->d_prevec_, dim_); }
+            if (this->dev_ == base_device::GpuDevice) {
+                resmem_real_op()({}, this->d_prevec_, dim_);
+            }
 #endif
         }
         PreOP(const PreOP&) = delete;
@@ -119,13 +121,13 @@ namespace hsolver
         fvec::Div<T> get() const
         {
             return std::bind(PreOP<T, Device, Kernel_t>::transvec_,
-                std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, this->prevec_);
+                std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, 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::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5,
                 this->prevec_, this->d_prevec_, this->transval_);
         }
 
diff --git a/source/module_lr/hsolver_lrtd.hpp b/source/module_lr/hsolver_lrtd.hpp
@@ -76,8 +76,8 @@ namespace LR
 
                 if (method == "dav")
                 {
-                    auto pre_func = [&precondition](T* ptr, const int& ld, const int& nvec)->void
-                        {  hsolver::fvec::div_prevec(ptr, ld, nvec, precondition.data()); };
+                    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
@@ -90,8 +90,8 @@ namespace LR
                 }
                 else if (method == "dav_subspace") //need refactor
                 {
-                    auto pre_func = [&precondition](T* ptr, const Real<T>* eig, const int& ld, const int& nvec)->void
-                        {  hsolver::fvec::div_trans_prevec_minus_eigen(ptr, eig, ld, nvec, precondition.data()); };
+                    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,

Original file line number	Diff line number	Diff line change
`@@ -26,19 +26,19 @@ namespace hsolver`
`26`	`26`	`/// type 1: directly divide each vector by the precondition vector`
`27`	`27`	`///---------------------------------------------------------------------------------------------`
`28`	`28`	`template <typename T, typename Device = base_device::DEVICE_CPU>`
`29`		`- void div_prevec(T* ptr, const size_t& dim, const size_t& nvec,`
	`29`	`+ void div_prevec(T* const dst, const T* const src, const size_t& dim, const size_t& nvec,`
`30`	`30`	`const Real<T>* const pre)`
`31`	`31`	`{`
`32`	`32`	`Device* ctx = {};`
`33`		`- for (int m = 0; m < nvec; m++)`
	`33`	`+ for (size_t m = 0; m < nvec; m++)`
`34`	`34`	`{`
`35`		`- T* const ptr_m = ptr + m * dim;`
`36`		`- vector_div_vector_op<T, Device>()(ctx, dim, ptr_m, ptr_m, pre);`
	`35`	`+ const size_t offset = m * dim;`
	`36`	`+ vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, pre);`
`37`	`37`	`}`
`38`	`38`	`}`
`39`	`39`	`/// Intereface to be called in the eigensolver`
`40`	`40`	`template <typename T>`
`41`		`- using Div = std::function<void(T*, const size_t&, const size_t&)>;`
	`41`	`+ using Div = std::function<void(T* const, const T* const, const size_t&, const size_t&)>;`
`42`	`42`	`// Kernel function full of dependence`
`43`	`43`	`template <typename T, typename Device = base_device::DEVICE_CPU>`
`44`	`44`	`using DivKernel = std::function<decltype(div_prevec<T, Device>)>;`
`@@ -48,7 +48,7 @@ namespace hsolver`
`48`	`48`	`///$X \to (A-\lambda I)^{-1} X$`
`49`	`49`	`///---------------------------------------------------------------------------------------------`
`50`	`50`	`template <typename T, typename Device = base_device::DEVICE_CPU>`
`51`		`- void div_trans_prevec_minus_eigen(T* ptr, const Real<T>* eig, const size_t& dim, const size_t& nvec,`
	`51`	`+ 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,`
`52`	`52`	`const Real<T>* const pre, Real<T>* const d_pre = nullptr, const std::function<Real<T>(const Real<T>&)>& transval = fval::none<Real<T>>)`
`53`	`53`	`{`
`54`	`54`	`using syncmem_var_h2d_op = base_device::memory::synchronize_memory_op<Real<T>, Device, base_device::DEVICE_CPU>;`
`@@ -57,27 +57,27 @@ namespace hsolver`
`57`	`57`	`Device* ctx = {};`
`58`	`58`	`base_device::DEVICE_CPU* cpu_ctx = {};`
`59`	`59`
`60`		`- for (int m = 0; m < nvec; m++)`
	`60`	`+ for (size_t m = 0; m < nvec; m++)`
`61`	`61`	`{`
`62`		`- T* const ptr_m = ptr + m * dim;`
	`62`	`+ const size_t offset = m * dim;`
`63`	`63`	`for (size_t i = 0; i < dim; i++) { pre_trans[i] = transval(pre[i] - eig[m]); }`
`64`	`64`	`#if defined(__CUDA) \|\| defined(__ROCM)`
`65`	`65`	`if (device == base_device::GpuDevice)`
`66`	`66`	`{`
`67`	`67`	`assert(d_pre);`
`68`	`68`	`syncmem_var_h2d_op()(ctx, cpu_ctx, d_pre, pre_trans.data(), dim);`
`69`		`- vector_div_vector_op<T, Device>()(ctx, dim, ptr_m, ptr_m, d_pre);`
	`69`	`+ vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, d_pre);`
`70`	`70`	`}`
`71`	`71`	`else`
`72`	`72`	`#endif`
`73`	`73`	`{`
`74`		`- vector_div_vector_op<T, Device>()(ctx, dim, ptr_m, ptr_m, pre_trans.data());`
	`74`	`+ vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, pre_trans.data());`
`75`	`75`	`}`
`76`	`76`	`}`
`77`	`77`	`}`
`78`	`78`	`/// Intereface to be called in the eigensolver`
`79`	`79`	`template <typename T>`
`80`		`- using DivTransMinusEig = std::function<void(T, const Real<T>, const size_t&, const size_t&)>;`
	`80`	`+ using DivTransMinusEig = std::function<void(T* const, const T* const, const Real<T>*, const size_t&, const size_t&)>;`
`81`	`81`	`/// Kernel function full of dependence`
`82`	`82`	`template <typename T, typename Device = base_device::DEVICE_CPU>`
`83`	`83`	`using DivTransMinusEigKernel = std::function<decltype(div_trans_prevec_minus_eigen<T, Device>)>;`
`@@ -104,7 +104,9 @@ namespace hsolver`
`104`	`104`	`dev_(base_device::get_device_type<Device>({}))`
`105`	`105`	`{`
`106`	`106`	`#if defined(__CUDA) \|\| defined(__ROCM)`
`107`		`- if (this->dev_ == base_device::GpuDevice) { resmem_real_op()({}, this->d_prevec_, dim_); }`
	`107`	`+ if (this->dev_ == base_device::GpuDevice) {`
	`108`	`+ resmem_real_op()({}, this->d_prevec_, dim_);`
	`109`	`+ }`
`108`	`110`	`#endif`
`109`	`111`	`}`
`110`	`112`	`PreOP(const PreOP&) = delete;`
`@@ -119,13 +121,13 @@ namespace hsolver`
`119`	`121`	`fvec::Div<T> get() const`
`120`	`122`	`{`
`121`	`123`	`return std::bind(PreOP<T, Device, Kernel_t>::transvec_,`
`122`		`- std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, this->prevec_);`
	`124`	`+ std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, this->prevec_);`
`123`	`125`	`}`
`124`	`126`	`template<typename U = Kernel_t, typename std::enable_if<std::is_same<U, fvec::DivTransMinusEigKernel<T, Device>>::value, bool>::type = 0>`
`125`	`127`	`fvec::DivTransMinusEig<T> get() const`
`126`	`128`	`{`
`127`	`129`	`return std::bind(PreOP<T, Device, Kernel_t>::transvec_,`
`128`		`- std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4,`
	`130`	`+ std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5,`
`129`	`131`	`this->prevec_, this->d_prevec_, this->transval_);`
`130`	`132`	`}`
`131`	`133`