Remove Base class DiagH in LCAO code (#5239)

* Remove Base DiagH class of DiagoLapack

* Remove DiagH Base class of DiagoScalapack and lcao test

* Remove DiagH Base class of DiagoCusolver and lcao_cusolver_test

* Update docs in hsolver

* Remove Base DiagH pointer in diago_lapack_test

* Remove Base DiagH of DiagoElpa

* Remove Base DiagH of DiagoElpaNative

* Remove Base DiagH of DiagoCusolverMP

---------

Co-authored-by: Haozhi Han <[email protected]>

Author: Cstandardlib

source/module_hsolver/diago_cusolver.h

@@ -9,9 +9,9 @@
 namespace hsolver
 {
 
-// DiagoCusolver class, derived from DiagH, for diagonalization using CUSOLVER
+// DiagoCusolver class for diagonalization using CUSOLVER
 template <typename T>
-class DiagoCusolver : public DiagH<T>
+class DiagoCusolver
 {
   private:
     // Real is the real part of the complex type T
@@ -24,7 +24,7 @@ class DiagoCusolver : public DiagH<T>
     ~DiagoCusolver();
 
     // Override the diag function for CUSOLVER diagonalization
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 
     // Static variable to keep track of the decomposition state
     static int DecomposedState;

source/module_hsolver/diago_cusolvermp.h

@@ -10,7 +10,7 @@ namespace hsolver
 {
 // DiagoCusolverMP class, derived from DiagH, for diagonalization using CUSOLVERMP
 template <typename T>
-class DiagoCusolverMP : public DiagH<T>
+class DiagoCusolverMP
 {
   private:
     using Real = typename GetTypeReal<T>::type;
@@ -19,8 +19,8 @@ class DiagoCusolverMP : public DiagH<T>
     DiagoCusolverMP()
     {
     }
-    // Override the diag function for CUSOLVERMP diagonalization
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    // the diag function for CUSOLVERMP diagonalization
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 };
 } // namespace hsolver
 #endif // __CUSOLVERMP

source/module_hsolver/diago_elpa.h

@@ -8,16 +8,16 @@ namespace hsolver
 {
 
 template <typename T>
-class DiagoElpa : public DiagH<T>
+class DiagoElpa
 {
   private:
     using Real = typename GetTypeReal<T>::type;
 
   public:
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
-    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm) override;
+    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
     MPI_Comm setmpicomm(); // set mpi comm;
     static int elpa_num_thread;  // need to set mpi_comm or not,-1 not,else the number of mpi needed
 #endif

source/module_hsolver/diago_elpa_native.h

@@ -8,16 +8,16 @@ namespace hsolver
 {
 
 template <typename T>
-class DiagoElpaNative : public DiagH<T>
+class DiagoElpaNative
 {
   private:
     using Real = typename GetTypeReal<T>::type;
 
   public:
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
-    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm) override;
+    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
     MPI_Comm setmpicomm(); // set mpi comm;
     static int elpa_num_thread;  // need to set mpi_comm or not,-1 not,else the number of mpi needed
     static int lastmpinum; // last using mpi;

source/module_hsolver/diago_lapack.h

@@ -20,13 +20,13 @@
 namespace hsolver
 {
 template <typename T>
-class DiagoLapack : public DiagH<T>
+class DiagoLapack
 {
   private:
     using Real = typename GetTypeReal<T>::type;
 
   public:
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 
     void dsygvx_diag(const int ncol,
                      const int nrow,

source/module_hsolver/diago_scalapack.h

@@ -20,15 +20,15 @@
 namespace hsolver
 {
     template<typename T>
-    class DiagoScalapack : public DiagH<T>
+    class DiagoScalapack
 {
 private:
     using Real = typename GetTypeReal<T>::type;
   public:
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in) override;
+    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
-    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm) override;
+    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
 #endif
 
   private:

source/module_hsolver/hsolver_lcao.cpp

@@ -169,7 +169,7 @@ void HSolverLCAO<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>&
         DiagoLapack<T> la;
         la.diag(hm, psi, eigenvalue);
 #else
-        ModuleBase::WARNING_QUIT("HSolverLCAO::solve", "This method of DiagH is not supported!");
+        ModuleBase::WARNING_QUIT("HSolverLCAO::solve", "This type of eigensolver is not supported!");
 #endif
     }
     else
@@ -379,7 +379,7 @@ void HSolverLCAO<T, Device>::parakSolve(hamilt::Hamilt<T>* pHamilt,
             else
             {
                 ModuleBase::WARNING_QUIT("HSolverLCAO::solve",
-                                         "This method of DiagH for k-parallelism diagnolization is not supported!");
+                                         "This type of eigensolver for k-parallelism diagnolization is not supported!");
             }
         }
         MPI_Barrier(MPI_COMM_WORLD);

source/module_hsolver/hsolver_pw_sdft.cpp

@@ -33,7 +33,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt<std::complex<double>>* pHamilt,
     const std::initializer_list<std::string> _methods = {"cg", "dav", "dav_subspace", "bpcg"};
     if (std::find(std::begin(_methods), std::end(_methods), this->method) == std::end(_methods))
     {
-        ModuleBase::WARNING_QUIT("HSolverPW::solve", "This method of DiagH is not supported!");
+        ModuleBase::WARNING_QUIT("HSolverPW::solve", "This type of eigensolver is not supported!");
     }
 
     // part of KSDFT to get KS orbitals

source/module_hsolver/test/diago_lapack_test.cpp

@@ -119,15 +119,15 @@ class DiagoLapackPrepare
         : nlocal(nlocal), nbands(nbands), nb2d(nb2d), sparsity(sparsity), hfname(hfname),
           sfname(sfname), solutionfname(solutionfname)
     {
-        dh = new hsolver::DiagoLapack<T>;
+        // dh = new hsolver::DiagoLapack<T>;
     }
 
     int nlocal, nbands, nb2d, sparsity;
     std::string sfname, hfname, solutionfname;
     std::vector<T> h;
     std::vector<T> s;
     HamiltTEST<T> hmtest;
-    hsolver::DiagH<T>* dh = nullptr;
+    // hsolver::DiagH<T>* dh = nullptr;
     psi::Psi<T> psi;
     std::vector<double> e_solver;
     std::vector<double> e_lapack;
@@ -200,9 +200,11 @@ class DiagoLapackPrepare
 
         for (int i = 0; i < REPEATRUN; i++)
         {
-            dh->diag(&hmtest, psi, e_solver.data());
+            hsolver::DiagoLapack<T> dh;
+            dh.diag(&hmtest, psi, e_solver.data());
+            // dh->diag(&hmtest, psi, e_solver.data());
         }
-        delete dh;
+        // delete dh;
     }
 
     void read_SOLUTION()

source/module_hsolver/test/diago_lcao_cusolver_test.cpp

@@ -74,10 +74,10 @@ class DiagoPrepare
         MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
 
         if (ks_solver == "scalapack_gvx")
-            dh = new hsolver::DiagoScalapack<T>;
+            ;//dh = new hsolver::DiagoScalapack<T>;
 #ifdef __CUDA
         else if (ks_solver == "cusolver")
-            dh = new hsolver::DiagoCusolver<T>;
+            ;//dh = new hsolver::DiagoCusolver<T>;
 #endif
         else
         {
@@ -96,7 +96,7 @@ class DiagoPrepare
     std::vector<T> s;
     std::vector<T> h_local;
     std::vector<T> s_local;
-    hsolver::DiagH<T>* dh = 0;
+    // hsolver::DiagH<T>* dh = 0;
     psi::Psi<T> psi;
     std::vector<double> e_solver;
     std::vector<double> e_lapack;
@@ -222,11 +222,23 @@ class DiagoPrepare
         {
             hmtest.h_local = this->h_local;
             hmtest.s_local = this->s_local;
-            dh->diag(&hmtest, psi, e_solver.data());
+            if (ks_solver == "scalapack_gvx")
+            {
+                hsolver::DiagoScalapack<T> dh;
+                dh.diag(&hmtest, psi, e_solver.data());
+            }
+    #ifdef __CUDA
+            else if (ks_solver == "cusolver")
+                {
+                    hsolver::DiagoCusolver<T> dh;
+                    dh.diag(&hmtest, psi, e_solver.data());
+                }
+    #endif
+            // dh->diag(&hmtest, psi, e_solver.data());
         }
         endtime = MPI_Wtime();
         hsolver_time = (endtime - starttime) / REPEATRUN;
-        delete dh;
+        // delete dh;
     }
 
     void diago_lapack()