Refactor&Feature: templatize ESolver_KS_LCAO , unify ops/opsd, psi/psid; real-valued Davidson

source/Makefile.Objects

@@ -207,6 +207,7 @@ OBJS_ESOLVER=esolver.o\
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_elec.o\
       esolver_ks_lcao_tddft.o\
+      esolver_ks_lcao_tmpfunc.o\
 
 OBJS_GINT=gint.o\
       gint_gamma.o\

source/module_base/blas_connector.h

@@ -193,9 +193,16 @@ class BlasConnector
 		zgemm_(&transb, &transa, &n, &m, &k,
 			&alpha, b, &ldb, a, &lda,
 			&beta, c, &ldc);
-	}
+    }
     static inline
     void gemv(const char trans, const int m, const int n,
+        const double alpha, const double* A, const int lda, const double* X, const int incx,
+        const double beta, double* Y, const int incy)
+    {
+        dgemv_(&trans, &m, &n, &alpha, A, &lda, X, &incx, &beta, Y, &incy);
+    }
+    static inline
+        void gemv(const char trans, const int m, const int n,
               const std::complex<float> alpha, const std::complex<float> *A, const int lda, const std::complex<float> *X, const int incx,
               const std::complex<float> beta, std::complex<float> *Y, const int incy)
     {

source/module_base/global_file.cpp

@@ -71,7 +71,7 @@ void ModuleBase::Global_File::make_dir_out(
 			}
         }
 #ifdef __MPI
-		Parallel_Reduce::reduce_int_all(make_dir);
+        Parallel_Reduce::reduce_all(make_dir);
 #endif
 		if(make_dir>0)break;
 		++times;
@@ -108,7 +108,7 @@ void ModuleBase::Global_File::make_dir_out(
                 }
             }
 #ifdef __MPI
-            Parallel_Reduce::reduce_int_all(make_dir_stru);
+            Parallel_Reduce::reduce_all(make_dir_stru);
 #endif
             if(make_dir_stru>0) break;
             ++times;
@@ -147,7 +147,7 @@ void ModuleBase::Global_File::make_dir_out(
                 }
             }
 #ifdef __MPI
-            Parallel_Reduce::reduce_int_all(make_dir_matrix);
+            Parallel_Reduce::reduce_all(make_dir_matrix);
 #endif
             if(make_dir_matrix>0) break;
             ++times;

source/module_base/global_function_ddotreal.cpp

@@ -22,7 +22,7 @@ namespace GlobalFunc
         pL=(double *)psi_L;
         pR=(double *)psi_R;
         double result=BlasConnector::dot(dim2,pL,1,pR,1);
-        if(reduce)  Parallel_Reduce::reduce_double_pool( result );
+        if (reduce)  Parallel_Reduce::reduce_pool(result);
         return result;
         //======================================================================
         /*std::complex<double> result(0,0);

source/module_base/lapack_connector.h

@@ -24,6 +24,11 @@ extern "C"
 
     // 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,
+        double* a, const int* lda,
+        const double* b, const int* ldb, double* w,
+        double* work, int* lwork,
+        int* iwork, int* liwork, int* info);
 
     void chegvd_(const int* itype, const char* jobz, const char* uplo, const int* n,
              std::complex<float>* a, const int* lda,
@@ -37,6 +42,12 @@ extern "C"
                  std::complex<double>* work, int* lwork, double* rwork, int* lrwork,
                  int* iwork, int* liwork, int* info);
 
+    void dsyevx_(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,
+        double* work, const int* lwork, double* rwork, int* iwork, int* ifail, int* info);
+
     void cheevx_(const char* jobz, const char* range, const char* uplo, const int* n,
              std::complex<float> *a, const int* lda,
              const float* vl, const float* vu, const int* il, const int* iu, const float* abstol,
@@ -313,6 +324,20 @@ class LapackConnector
                 iwork, &liwork, &info);
     }
 
+    // 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,
+            double* a, const int lda,
+            const double* b, const int ldb, double* w,
+            double* work, int lwork, double* rwork, int lrwork,
+            int* iwork, int liwork, int& info)
+    {
+        // call the fortran routine
+        dsygvd_(&itype, &jobz, &uplo, &n,
+            a, &lda, b, &ldb, w,
+            work, &lwork,
+            iwork, &liwork, &info);
+    }
 
     // wrap function of fortran lapack routine zhegvd. (pointer version)
     static inline
@@ -358,7 +383,21 @@ class LapackConnector
                 work, &lwork, rwork, iwork, ifail, &info);
     }
 
-    // wrap function of fortran lapack routine zheevx.
+    // 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,
+            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,
+            double* work, const int lwork, double* rwork, int* iwork, int* ifail, int& info)
+    {
+        dsyevx_(&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 cheevx.
     static inline
     void xheevx( const int itype, const char jobz, const char range, const char uplo, const int n,
                  std::complex<float>* a, const int lda,

source/module_base/memory.cpp

@@ -235,7 +235,7 @@ void Memory::print_all(std::ofstream &ofs)
 
 	const double small = 1.0; 
 #ifdef __MPI
-		Parallel_Reduce::reduce_double_all(Memory::total);
+    Parallel_Reduce::reduce_all(Memory::total);
 #endif
     ofs <<"\n NAME---------------|MEMORY(MB)--------" << std::endl;
 //	std::cout<<"\n"<<std::setw(41)<< " " <<std::setprecision(4)<<total;
@@ -262,7 +262,7 @@ void Memory::print_all(std::ofstream &ofs)
 		}
 		print_flag[k] = true;
 #ifdef __MPI
-		Parallel_Reduce::reduce_double_all(consume[k]);
+        Parallel_Reduce::reduce_all(consume[k]);
 #endif
 	    if ( consume[k] < small ) 
         {

source/module_base/opt_CG.cpp

@@ -134,7 +134,7 @@ double Opt_CG::step_length(
 )
 {
     double dAd = this->inner_product(pdirect, pAd, this->nx);
-    Parallel_Reduce::reduce_double_all(dAd);
+    Parallel_Reduce::reduce_all(dAd);
     ifPD = 0;
     // check for positive-definiteness, very important for convergence
     if (dAd == 0)
@@ -183,7 +183,7 @@ void Opt_CG::stantard_CGdirect(
             temp_gradient[i] = this->pgradient_old[i] + this->alpha * pAd[i];
         }
         this->beta = this->inner_product(temp_gradient, temp_gradient, this->nx) / this->gg;
-        Parallel_Reduce::reduce_double_all(this->beta);
+        Parallel_Reduce::reduce_all(this->beta);
         for (int i = 0; i < this->nx; ++i)
         {
             this->pgradient_old[i] = temp_gradient[i];
@@ -193,7 +193,7 @@ void Opt_CG::stantard_CGdirect(
         delete[] temp_gradient;
     }
     this->gg = this->inner_product(this->pgradient_old, this->pgradient_old, this->nx);
-    Parallel_Reduce::reduce_double_all(this->gg);
+    Parallel_Reduce::reduce_all(this->gg);
     this->iter++;
 }
 
@@ -209,9 +209,9 @@ void Opt_CG::PR_beta(
     double temp_beta = 0.;
     temp_beta = this->inner_product(pgradient, pgradient, this->nx);
     temp_beta -= this->inner_product(pgradient, this->pgradient_old, this->nx);
-    Parallel_Reduce::reduce_double_all(temp_beta);
+    Parallel_Reduce::reduce_all(temp_beta);
     double gg_old = this->inner_product(this->pgradient_old, this->pgradient_old, this->nx);
-    Parallel_Reduce::reduce_double_all(gg_old);
+    Parallel_Reduce::reduce_all(gg_old);
     // temp_beta /= this->inner_product(this->pgradient_old, this->pgradient_old, this->nx);
     temp_beta /= gg_old;
     this->beta = std::max(0., temp_beta);
@@ -230,19 +230,19 @@ void Opt_CG::HZ_beta(
     for (int i = 0; i < this->nx; ++i) y[i] = pgradient[i] - this->pgradient_old[i];
 
     double py = this->inner_product(this->pdirect_old, y, this->nx);
-    Parallel_Reduce::reduce_double_all(py);
+    Parallel_Reduce::reduce_all(py);
     double yy = this->inner_product(y, y, this->nx);
-    Parallel_Reduce::reduce_double_all(yy);
+    Parallel_Reduce::reduce_all(yy);
     double pg = this->inner_product(this->pdirect_old, pgradient, this->nx);
-    Parallel_Reduce::reduce_double_all(pg);
+    Parallel_Reduce::reduce_all(pg);
     double yg = this->inner_product(y, pgradient, this->nx);
-    Parallel_Reduce::reduce_double_all(yg);
+    Parallel_Reduce::reduce_all(yg);
     double temp_beta = (yg - 2 * pg * yy / py) /py;
 
     double pp = this->inner_product(this->pdirect_old, this->pdirect_old, this->nx);
-    Parallel_Reduce::reduce_double_all(pp);
+    Parallel_Reduce::reduce_all(pp);
     double gg = this->inner_product(this->pgradient_old, this->pgradient_old, this->nx);
-    Parallel_Reduce::reduce_double_all(gg);
+    Parallel_Reduce::reduce_all(gg);
     double temp_eta = -1 / (sqrt(pp) * std::min(this->eta, sqrt(gg)));
 
     this->beta = std::max(temp_beta, temp_eta);

source/module_base/opt_TN.hpp

@@ -100,9 +100,9 @@ class Opt_TN
     {
         double epsilon = 0.;
         double xx = this->inner_product(px, px, this->nx);
-        Parallel_Reduce::reduce_double_all(xx);
+        Parallel_Reduce::reduce_all(xx);
         double dd = this->inner_product(pcgDirect, pcgDirect, this->nx);
-        Parallel_Reduce::reduce_double_all(dd);
+        Parallel_Reduce::reduce_all(dd);
         epsilon = 2 * sqrt(this->machPrec) * (1 + sqrt(xx)) / sqrt(dd);
         // epsilon = 2 * sqrt(this->machPrec) * (1 + sqrt(this->inner_product(px, px, this->nx))) 
         //         / sqrt(this->inner_product(pcgDirect, pcgDirect, this->nx));

source/module_base/parallel_reduce.cpp

@@ -4,7 +4,8 @@
 #include "module_base/global_variable.h"
 #include "parallel_global.h"
 
-void Parallel_Reduce::reduce_int_all(int &object)
+template<>
+void Parallel_Reduce::reduce_all<int>(int& object)
 {
 #ifdef __MPI
     MPI_Allreduce(MPI_IN_PLACE , &object , 1, MPI_INT , MPI_SUM , MPI_COMM_WORLD);
@@ -21,16 +22,17 @@ void Parallel_Reduce::reduce_int_diag(int &object)
 }
 
 
-
-void Parallel_Reduce::reduce_double_all(double &object)
+template<>
+void Parallel_Reduce::reduce_all<double>(double& object)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE , &object , 1, MPI_DOUBLE , MPI_SUM , MPI_COMM_WORLD);
 #endif
 	return;
 }
 
-void Parallel_Reduce::reduce_int_all(int *object, const int n)
+template<>
+void Parallel_Reduce::reduce_all<int>(int* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
@@ -48,8 +50,8 @@ void Parallel_Reduce::reduce_int_grid(int *object, const int n)
 }
 
 
-
-void Parallel_Reduce::reduce_double_all(double *object, const int n)
+template<>
+void Parallel_Reduce::reduce_all<double>(double* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
@@ -76,24 +78,26 @@ void Parallel_Reduce::reduce_double_diag(double *object, const int n)
 	return;
 }
 
-
-void Parallel_Reduce::reduce_double_pool(float &object)
+template<>
+void Parallel_Reduce::reduce_pool<float>(float& object)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE , &object , 1, MPI_FLOAT , MPI_SUM , POOL_WORLD);
 #endif
     return;
 }
 
-void Parallel_Reduce::reduce_double_pool(double &object)
+template<>
+void Parallel_Reduce::reduce_pool<double>(double& object)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE , &object , 1, MPI_DOUBLE , MPI_SUM , POOL_WORLD);
 #endif
 	return;
 }
 
-void Parallel_Reduce::reduce_double_pool(double *object, const int n)
+template<>
+void Parallel_Reduce::reduce_pool<double>(double* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE, MPI_SUM, POOL_WORLD);
@@ -130,8 +134,8 @@ void Parallel_Reduce::reduce_double_allpool(double *object, const int n)
 }
 
 
-
-void Parallel_Reduce::reduce_complex_double_all(std::complex<double> &object)
+template<>
+void Parallel_Reduce::reduce_all<std::complex<double>>(std::complex<double>& object)
 {
 #ifdef __MPI
     MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
@@ -140,31 +144,35 @@ void Parallel_Reduce::reduce_complex_double_all(std::complex<double> &object)
 }
 
 //LiuXh add 2019-07-16
-void Parallel_Reduce::reduce_complex_double_all(std::complex <double> *object, const int n)
+template<>
+void Parallel_Reduce::reduce_all<std::complex<double>>(std::complex<double>* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
 #endif
 	return;
 }
 
-void Parallel_Reduce::reduce_complex_double_pool(std::complex<double> &object)
+template<>
+void Parallel_Reduce::reduce_pool<std::complex<double>>(std::complex<double>& object)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD);
 #endif
 	return;
 }
 
-void Parallel_Reduce::reduce_complex_double_pool(std::complex <float> *object, const int n)
+template<>
+void Parallel_Reduce::reduce_pool<std::complex<float>>(std::complex <float>* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_C_FLOAT_COMPLEX, MPI_SUM, POOL_WORLD);
 #endif
     return;
 }
 
-void Parallel_Reduce::reduce_complex_double_pool(std::complex <double> *object, const int n)
+template<>
+void Parallel_Reduce::reduce_pool<std::complex<double>>(std::complex <double>* object, const int n)
 {
 #ifdef __MPI
 	MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD);
@@ -267,4 +275,4 @@ void Parallel_Reduce::gather_min_double_all(double &v)
 	}
 	delete[] value;
 #endif
-}
+}