Feature: Add a new Davidson iteration method called subspace davidson for pw basis (#3903)

* add new_dav method which is same to davidson method for pw basis

* add new_dav method for pw basis which is more efficient than origin dav method

* fix ilaenv interface bug

* update new_dav method 3.5

* implement new_dav method for pw basis (cpu version, one core)

* fix bug of dav method for pw basis

* debug for new davidson method

* opt some value setting for new_dav files

* format and reorganize the code

* fix CUDA compile bug

* format diago_newdav.cpp

* Implement multi-core parallelism of the new davidson method

* fix build bug for without mpi

* replace new-dav of subspace-dav

* change file name from diago_newdav to diago_subspacedav

* fix build bug for tests

* change the name of subspacedav to dav_subspace

* fix build bug in Integration Test

---------

Co-authored-by: Mohan Chen <[email protected]>

Author: haozhihan

source/Makefile.Objects

@@ -286,6 +286,7 @@ OBJS_HCONTAINER=base_matrix.o\
 
 OBJS_HSOLVER=diago_cg.o\
     diago_david.o\
+    diago_dav_subspace.o\
     diagh_consts.o\
     diago_bpcg.o\
     hsolver_pw.o\

source/module_base/lapack_connector.h

@@ -22,6 +22,9 @@
 extern "C"
 {
 
+    int ilaenv_(int* ispec,const char* name,const char* opts,
+                const int* n1,const int* n2,const int* n3,const int* n4);
+
     // solve the generalized eigenproblem Ax=eBx, where A is Hermitian and complex couble
     // zhegv_ & zhegvd_ returns all eigenvalues while zhegvx_ returns selected ones
     void dsygvd_(const int* itype, const char* jobz, const char* uplo, const int* n,
@@ -60,9 +63,12 @@ extern "C"
                  const int* m, double* w, std::complex<double> *z, const int *ldz, 
                  std::complex<double> *work, const int* lwork, double* rwork, int* iwork, int* ifail, int* info);
 
-    void zhegv_(const int* itype,const char* jobz,const char* uplo,const int* n,
-                std::complex<double>* a,const int* lda,std::complex<double>* b,const int* ldb,
-                double* w,std::complex<double>* work,int* lwork,double* rwork,int* info);
+
+    void dsygvx_(const int* itype, const char* jobz, const char* range, const char* uplo,
+                const int* n, double* A, const int* lda, double* B, const int* ldb,
+                const double* vl, const double* vu, const int* il, const int* iu,
+                const double* abstol, const int* m, double* w, double* Z, const int* ldz,
+                double* work, int* lwork, int*iwork, int* ifail, int* info);
 
     void chegvx_(const int* itype,const char* jobz,const char* range,const char* uplo,
              const int* n,std::complex<float> *a,const int* lda,std::complex<float> *b,
@@ -78,6 +84,16 @@ extern "C"
                  std::complex<double> *z,const int *ldz,std::complex<double> *work,const int* lwork,
                  double* rwork,int* iwork,int* ifail,int* info);
 
+    void zhegv_(const int* itype,const char* jobz,const char* uplo,const int* n,
+                std::complex<double>* a,const int* lda,std::complex<double>* b,const int* ldb,
+                double* w,std::complex<double>* work,int* lwork,double* rwork,int* info);
+    void chegv_(const int* itype,const char* jobz,const char* uplo,const int* n,
+                std::complex<float>* a,const int* lda,std::complex<float>* b,const int* ldb,
+                float* w,std::complex<float>* work,int* lwork,float* rwork,int* info);
+	void dsygv_(const int* itype, const char* jobz,const char* uplo, const int* n,
+				double* a,const int* lda,double* b,const int* ldb,
+	 			double* w,double* work,int* lwork,int* info);
+
     // solve the eigenproblem Ax=ex, where A is Hermitian and complex couble
     // zheev_ returns all eigenvalues while zheevx_ returns selected ones
     void zheev_(const char* jobz,const char* uplo,const int* n,std::complex<double> *a,
@@ -86,18 +102,6 @@ extern "C"
     void cheev_(const char* jobz,const char* uplo,const int* n,std::complex<float> *a,
                 const int* lda,float* w,std::complex<float >* work,const int* lwork,
                 float* rwork,int* info);
-
-	// solve the generalized eigenproblem Ax=eBx, where A is Symmetric and real couble
-    // dsygv_ returns all eigenvalues while dsygvx_ returns selected ones
-	void dsygv_(const int* itype, const char* jobz,const char* uplo, const int* n,
-				double* a,const int* lda,double* b,const int* ldb,
-	 			double* w,double* work,int* lwork,int* info);
-	void dsygvx_(const int* itype, const char* jobz, const char* range, const char* uplo,
-                const int* n, double* A, const int* lda, double* B, const int* ldb,
-                const double* vl, const double* vu, const int* il, const int* iu,
-                const double* abstol, int* m, double* w, double* Z, const int* ldz,
-                double* work, int* lwork, int*iwork, int* ifail, int* info);
-    // solve the eigenproblem Ax=ex, where A is Symmetric and real double
 	void dsyev_(const char* jobz,const char* uplo,const int* n,double *a,
                 const int* lda,double* w,double* work,const int* lwork, int* info);
 
@@ -314,23 +318,19 @@ class LapackConnector
 	}
 
 public:
-    // wrap function of fortran lapack routine zhegvd.
     static inline
-    void zhegvd(const int itype, const char jobz, const char uplo, const int n, 
-                std::complex<double>* a, const int lda, 
-                const std::complex<double>* b, const int ldb, double* w, 
-                std::complex<double>* work, int lwork, double* rwork, int lrwork,
-                int* iwork, int liwork, int& info)
+    int ilaenv( int ispec, const char *name,const char *opts,const int n1,const int n2,
+                const int n3,const int n4)
     {
-        zhegvd_(&itype, &jobz, &uplo, &n, 
-                a, &lda, b, &ldb, w,
-                work, &lwork, rwork, &lrwork,
-                iwork, &liwork, &info);
+        const int nb = ilaenv_(&ispec, name, opts, &n1, &n2, &n3, &n4);
+        return nb;
     }
 
+
+
     // wrap function of fortran lapack routine zhegvd. (pointer version)
     static inline
-        void xhegvd(const int itype, const char jobz, const char uplo, const int n,
+    void xhegvd(const int itype, const char jobz, const char uplo, const int n,
             double* a, const int lda,
             const double* b, const int ldb, double* w,
             double* work, int lwork, double* rwork, int lrwork,
@@ -373,23 +373,9 @@ class LapackConnector
                 iwork, &liwork, &info);
     }
 
-    // wrap function of fortran lapack routine zheevx.
-    static inline
-    void zheevx( const int itype, const char jobz, const char range, const char uplo, const int n, 
-                 std::complex<double>* a, const int lda, 
-                 const double vl, const double vu, const int il, const int iu, const double abstol, 
-                 const int m, double* w, std::complex<double>* z, const int ldz,
-                 std::complex<double>* work, const int lwork, double* rwork, int* iwork, int* ifail, int& info)
-    {
-        zheevx_(&jobz, &range, &uplo, &n, 
-                a, &lda, &vl, &vu, &il, &iu, 
-                &abstol, &m, w, z, &ldz, 
-                work, &lwork, rwork, iwork, ifail, &info);
-    }
-
     // wrap function of fortran lapack routine dsyevx.
     static inline
-        void xheevx(const int itype, const char jobz, const char range, const char uplo, const int n,
+    void xheevx(const int itype, const char jobz, const char range, const char uplo, const int n,
             double* a, const int lda,
             const double vl, const double vu, const int il, const int iu, const double abstol,
             const int m, double* w, double* z, const int ldz,
@@ -428,6 +414,98 @@ class LapackConnector
                 &abstol, &m, w, z, &ldz,
                 work, &lwork, rwork, iwork, ifail, &info);
     }
+    
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegvx( const int itype, const char jobz, const char range, const char uplo,
+                 const int n, std::complex<float>* a, const int lda, std::complex<float>* b,
+                 const int ldb, const float vl, const float vu, const int il, const int iu,
+                 const float abstol, const int m, float* w, std::complex<float>* z, const int ldz,
+                 std::complex<float>* work, const int lwork, float* rwork, int* iwork,
+                 int* ifail, int& info)
+    {
+        chegvx_(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl,
+                &vu, &il,&iu, &abstol, &m, w, z, &ldz, work, &lwork, rwork, iwork, ifail, &info);
+    }
+
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegvx( const int itype, const char jobz, const char range, const char uplo,
+                 const int n, std::complex<double>* a, const int lda, std::complex<double>* b,
+                 const int ldb, const double vl, const double vu, const int il, const int iu,
+                 const double abstol, const int m, double* w, std::complex<double>* z, const int ldz,
+                 std::complex<double>* work, const int lwork, double* rwork, int* iwork,
+                 int* ifail, int& info)
+    {
+        zhegvx_(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl,
+                &vu, &il,&iu, &abstol, &m, w, z, &ldz, work, &lwork, rwork, iwork, ifail, &info);
+    }
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegvx( const int itype, const char jobz, const char range, const char uplo,
+                 const int n, double* a, const int lda, double* b,
+                 const int ldb, const double vl, const double vu, const int il, const int iu,
+                 const double abstol, const int m, double* w, double* z, const int ldz,
+                 double* work, const int lwork, double* rwork, int* iwork,
+                 int* ifail, int& info)
+    {
+        // dsygvx_(&itype, &jobz, &range, &uplo, &n, a, &lda, b, &ldb, &vl,
+        //         &vu, &il,&iu, &abstol, &m, w, z, &ldz, work, &lwork, rwork, iwork, ifail, &info);
+    }
+
+
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegv( const int itype, const char jobz, const char uplo,
+                const int n, 
+                double* a, const int lda, 
+                double* b, const int ldb, 
+                double* w, 
+                double* work, int lwork, 
+                double* rwork, int& info)
+    {
+        // TODO
+    }
+
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegv( const int itype, const char jobz, const char uplo,
+                const int n, 
+                std::complex<float>* a, const int lda, 
+                std::complex<float>* b, const int ldb, 
+                float* w, 
+                std::complex<float>* work, int lwork, 
+                float* rwork, int& info)
+    {
+        // TODO
+    }
+    // wrap function of fortran lapack routine xhegvx ( pointer version ).
+    static inline
+    void xhegv( const int itype, const char jobz, const char uplo,
+                const int n, 
+                std::complex<double>* a, const int lda, 
+                std::complex<double>* b, const int ldb, 
+                double* w, 
+                std::complex<double>* work, int lwork, 
+                double* rwork, int& info)
+    {
+        zhegv_(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, w, work, &lwork, rwork, &info);
+    }
+
+
+    // wrap function of fortran lapack routine zhegvd.
+    static inline
+    void zhegvd(const int itype, const char jobz, const char uplo, const int n, 
+                std::complex<double>* a, const int lda, 
+                const std::complex<double>* b, const int ldb, double* w, 
+                std::complex<double>* work, int lwork, double* rwork, int lrwork,
+                int* iwork, int liwork, int& info)
+    {
+        zhegvd_(&itype, &jobz, &uplo, &n, 
+                a, &lda, b, &ldb, w,
+                work, &lwork, rwork, &lrwork,
+                iwork, &liwork, &info);
+    }
 
     // wrap function of fortran lapack routine zhegv ( ModuleBase::ComplexMatrix version ).
     static inline
@@ -543,30 +621,6 @@ class LapackConnector
         delete[] zux;
     }
 
-    // wrap function of fortran lapack routine xhegvx ( pointer version ).
-    static inline
-    void xhegvx( const int itype, const char jobz, const char range, const char uplo,
-                 const int n, const std::complex<float>* a, const int lda, const std::complex<float>* b,
-                 const int ldb, const float vl, const float vu, const int il, const int iu,
-                 const float abstol, const int m, float* w, std::complex<float>* z, const int ldz,
-                 std::complex<float>* work, const int lwork, float* rwork, int* iwork,
-                 int* ifail, int& info, int nbase_x)
-    {
-        chegvx(itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info, nbase_x);
-    }
-
-    // wrap function of fortran lapack routine xhegvx ( pointer version ).
-    static inline
-    void xhegvx( const int itype, const char jobz, const char range, const char uplo,
-                 const int n, const std::complex<double>* a, const int lda, const std::complex<double>* b,
-                 const int ldb, const double vl, const double vu, const int il, const int iu,
-                 const double abstol, const int m, double* w, std::complex<double>* z, const int ldz,
-                 std::complex<double>* work, const int lwork, double* rwork, int* iwork,
-                 int* ifail, int& info, int nbase_x)
-    {
-        zhegvx(itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, rwork, iwork, ifail, info, nbase_x);
-    }
-
 	// calculate the eigenvalues and eigenfunctions of a real symmetric matrix.
     static inline
     void dsygv(	const int itype,const char jobz,const char uplo,const int n,ModuleBase::matrix& a,

source/module_basis/module_ao/ORB_control.cpp

@@ -203,7 +203,7 @@ void ORB_control::setup_2d_division(std::ofstream& ofs_running,
 
     // determine whether 2d-division or not according to ks_solver
     bool div_2d;
-    if (ks_solver == "lapack" || ks_solver == "cg" || ks_solver == "dav") div_2d = false;
+    if (ks_solver == "lapack" || ks_solver == "cg" || ks_solver == "dav" || ks_solver == "dav_subspace") div_2d = false;
 #ifdef __MPI
     else if (ks_solver == "genelpa" || ks_solver == "scalapack_gvx" || ks_solver == "cusolver" || ks_solver == "cg_in_lcao") div_2d = true;
 #endif

source/module_elecstate/elecstate_print.cpp

@@ -289,6 +289,10 @@ void ElecState::print_etot(const bool converged,
     {
         label = "DA";
     }
+    else if (ks_solver_type == "dav_subspace")
+    {
+        label = "DS";
+    }
     else if (ks_solver_type == "scalapack_gvx")
     {
         label = "GV";

source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp

@@ -574,7 +574,7 @@ void diago_PAO_in_pw_k2(const psi::DEVICE_GPU *ctx,
             //GlobalC::hm.diagH_subspace(ik ,starting_nw, nbands, wfcatom, wfcatom, etatom.data());
         }
     }
-    else if(GlobalV::KS_SOLVER=="dav")
+    else if (GlobalV::KS_SOLVER == "dav" || GlobalV::KS_SOLVER == "dav_subspace")
     {
         assert(nbands <= wfcatom.nr);
         // replace by haozhihan 2022-11-23
@@ -685,8 +685,8 @@ void diago_PAO_in_pw_k2(const psi::DEVICE_GPU *ctx,
 			//GlobalC::hm.diagH_subspace(ik ,starting_nw, nbands, wfcatom, wfcatom, etatom.data());
 		}
 	}
-	else if(GlobalV::KS_SOLVER=="dav")
-	{
+    else if (GlobalV::KS_SOLVER == "dav" || GlobalV::KS_SOLVER == "dav_subspace")
+    {
 		assert(nbands <= wfcatom.nr);
 		// replace by haozhihan 2022-11-23
         hsolver::matrixSetToAnother<std::complex<double>, psi::DEVICE_GPU>()(

source/module_hsolver/CMakeLists.txt

@@ -2,6 +2,7 @@ list(APPEND objects
     diagh_consts.cpp
     diago_cg.cpp
     diago_david.cpp
+    diago_dav_subspace.cpp
     diago_bpcg.cpp
     hsolver_pw.cpp
     hsolver_pw_sdft.cpp