Merge branch 'develop' of https://gitee.com/deepmodeling/abacus-develop into HSolver

Author: dyzheng

source/Makefile

@@ -14,6 +14,7 @@ VPATH=./src_global\
 :./module_cell\
 :./module_base\
 :./module_md\
+:./module_surchem\
 :./module_deepks\
 :./module_xc\
 :./module_esolver\
@@ -68,6 +69,7 @@ FP_OBJS_0=$(OBJS_FIRST_PRINCIPLES)\
 $(OBJS_PARALLEL)\
 $(OBJS_TOOLS) \
 $(OBJS_COMMON)\
+$(OBJS_SURCHEM)\
 main.o\
 
 FP_OBJS=$(patsubst %.o, ${OBJ_DIR}/%.o, ${FP_OBJS_0})

source/Makefile.Objects

@@ -367,6 +367,15 @@ pzt2s.o\
 pdtrsm.o\
 pzhtrsm.o\
 
+OBJS_SURCHEM=cal_epsilon.o\
+cal_pseudo.o\
+cal_totn.o\
+cal_vcav.o\
+cal_vel.o\
+corrected_energy.o\
+H_correction_pw.o\
+minimize_cg.o\
+
 PDIAG_MR_0=dcopy.o\
 dlae2.o\
 dlaebz.o\

source/module_hsolver/diago_cg.cpp

@@ -37,6 +37,7 @@ void DiagoCG::diag_mock(psi::Psi<std::complex<double>> &phi, double *eigenvalue_
     this->dmx = phi.get_nbasis();
     this->n_band = phi.get_nbands();
     this->eigenvalue = eigenvalue_in;
+    ModuleBase::GlobalFunc::ZEROS(this->eigenvalue, this->n_band);
 
     /// record for how many loops in cg convergence
     double avg = 0.0;

source/module_hsolver/test/CMakeLists.txt

@@ -1,15 +1,31 @@
-remove_definitions(-D__MPI)
 remove_definitions(-D__CUDA)
 remove_definitions(-D__ROCM)
 AddTest(
   TARGET HSolver_cg
   LIBS ${math_libs} base 
-  SOURCES diago_cg_test.cpp ../diago_cg.cpp  ../diago_iter_assist.cpp ../../module_hamilt/hamilt_pw.cpp ../../module_psi/psi.cpp  ../../src_parallel/parallel_reduce.cpp
+  SOURCES diago_cg_test.cpp ../diago_cg.cpp  ../diago_iter_assist.cpp ../../module_hamilt/hamilt_pw.cpp 
+          ../../module_psi/psi.cpp  ../../src_parallel/parallel_reduce.cpp
+          ../../src_parallel/parallel_global.cpp ../../module_pw/unittest/test_tool.cpp
 )
 AddTest(
   TARGET HSolver_dav
   LIBS ${math_libs} base 
-  SOURCES diago_david_test.cpp ../diago_david.cpp  ../diago_iter_assist.cpp ../../module_hamilt/hamilt_pw.cpp ../../module_psi/psi.cpp  ../../src_parallel/parallel_reduce.cpp
+  SOURCES diago_david_test.cpp ../diago_david.cpp  ../diago_iter_assist.cpp ../../module_hamilt/hamilt_pw.cpp 
+          ../../module_psi/psi.cpp  ../../src_parallel/parallel_reduce.cpp ../../src_parallel/parallel_global.cpp 
+          ../../module_pw/unittest/test_tool.cpp
 )
 
-install(FILES data-H DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+install(FILES H-KPoints-Si64.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+install(FILES diago_cg_parallel_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+install(FILES diago_david_parallel_test.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+
+
+find_program(BASH bash)
+add_test(NAME HSolver_cg_parallel
+      COMMAND ${BASH} diago_cg_parallel_test.sh
+      WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+    )
+add_test(NAME HSolver_dav_parallel
+      COMMAND ${BASH} diago_david_parallel_test.sh
+      WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+    )  

source/module_hsolver/test/data-H renamed to source/module_hsolver/test/H-KPoints-Si64.dat

@@ -0,0 +1,13 @@
+#!/bin/bash
+
+np=`cat /proc/cpuinfo | grep "cpu cores" | uniq| awk '{print $NF}'`
+echo "nprocs in this machine is $np"
+
+for i in 6 3 2;do
+    if [[ $i -gt $np ]];then
+        continue
+    fi
+    echo "TEST DIAGO CG in parallel, nprocs=$i"
+    mpirun -np $i ./HSolver_cg
+    break    
+done

source/module_hsolver/test/diago_cg_parallel_test.sh

@@ -8,6 +8,7 @@
 #include "../diago_iter_assist.h"
 #include "diago_mock.h"
 #include "mpi.h"
+#include "module_pw/unittest/test_tool.h"
 #include <complex>
 
 #include "gtest/gtest.h"
@@ -57,21 +58,26 @@ class DiagoCGPrepare
         : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter),
           threshold(threshold)
     {
+#ifdef __MPI	
+		MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+        MPI_Comm_rank(MPI_COMM_WORLD, &mypnum);
+#endif	
     }
 
     int nband, npw, sparsity, maxiter, notconv;
     // eps is the convergence threshold within cg_diago
     double eps, avg_iter;
     bool reorder;
     double threshold;
+    int nprocs=1, mypnum=0;
     // threshold is the comparison standard between cg and lapack
 
     void CompareEigen(double *precondition)
     {
         // calculate eigenvalues by LAPACK;
         double *e_lapack = new double[npw];
         ModuleBase::ComplexMatrix ev = DIAGOTEST::hmatrix;
-        lapackEigen(npw, ev, e_lapack, false);
+        if(mypnum == 0)  lapackEigen(npw, ev, e_lapack, false);
         // initial guess of psi by perturbing lapack psi
         ModuleBase::ComplexMatrix psiguess(nband, npw);
         std::default_random_engine p(1);
@@ -80,7 +86,7 @@ class DiagoCGPrepare
         {
             for (int j = 0; j < npw; j++)
             {
-		double rand = static_cast<double>(u(p))/10.;
+		        double rand = static_cast<double>(u(p))/10.;
                 // psiguess(i,j) = ev(j,i)*(1+rand);
                 psiguess(i, j) = ev(j, i) * rand;
             }
@@ -89,24 +95,47 @@ class DiagoCGPrepare
 	//======================================================================
         double *en = new double[npw];
         int ik = 1;
-	hamilt::Hamilt* ha;
-	ha =new hamilt::HamiltPW;
-	Hamilt_PW* hpw;
-	int* ngk = new int [1];
-	//psi::Psi<std::complex<double>> psi(ngk,ik,nband,npw);
-	psi::Psi<std::complex<double>> psi;
-	psi.resize(ik,nband,npw);
-	//psi.fix_k(0);
+	    hamilt::Hamilt* ha;
+	    ha =new hamilt::HamiltPW;
+	    Hamilt_PW* hpw;
+	    int* ngk = new int [1];
+	    //psi::Psi<std::complex<double>> psi(ngk,ik,nband,npw);
+	    psi::Psi<std::complex<double>> psi;
+	    psi.resize(ik,nband,npw);
+	    //psi.fix_k(0);
         for (int i = 0; i < nband; i++)
         {
             for (int j = 0; j < npw; j++)
             {
-	        psi(i,j)=psiguess(i,j);
-	    }
-	}
-	hsolver::DiagoCG cg(hpw,precondition);
-        cg.diag(ha,psi,en); 
-	//======================================================================
+	            psi(i,j)=psiguess(i,j);
+	        }
+	    }	
+
+        psi::Psi<std::complex<double>> psi_local;
+        double* precondition_local;
+        DIAGOTEST::npw_local = new int[nprocs];
+#ifdef __MPI				
+	    DIAGOTEST::cal_division(DIAGOTEST::npw);
+	    DIAGOTEST::divide_hpsi(psi,psi_local); //will distribute psi and Hmatrix to each process
+	    precondition_local = new double[DIAGOTEST::npw_local[mypnum]];
+	    DIAGOTEST::divide_psi<double>(precondition,precondition_local);	
+#else
+	    DIAGOTEST::hmatrix_local = DIAGOTEST::hmatrix;
+	    DIAGOTEST::npw_local[0] = DIAGOTEST::npw;
+	    psi_local = psi;
+	    precondition_local = new double[DIAGOTEST::npw];
+	    for(int i=0;i<DIAGOTEST::npw;i++) precondition_local[i] = precondition[i];
+#endif
+        hsolver::DiagoCG cg(hpw,precondition_local);
+        psi_local.fix_k(0);
+        double start, end;
+        start = MPI_Wtime();
+        cg.diag(ha,psi_local,en); 
+        end = MPI_Wtime();
+        //if(mypnum == 0) printf("diago time:%7.3f\n",end-start);
+        delete [] DIAGOTEST::npw_local;
+	    delete [] precondition_local;
+	    //======================================================================
         for (int i = 0; i < nband; i++)
         {
             EXPECT_NEAR(en[i], e_lapack[i], threshold);
@@ -142,11 +171,11 @@ INSTANTIATE_TEST_SUITE_P(VerifyCG,
                          DiagoCGTest,
                          ::testing::Values(
                              // nband, npw, sparsity, reorder, eps, maxiter, threshold
-                             DiagoCGPrepare(10, 500, 0, true, 1e-5, 100, 1e-3),
+                             DiagoCGPrepare(10, 500, 0, true, 1e-5, 300, 1e-3),
                              DiagoCGPrepare(20, 500, 6, true, 1e-5, 300, 1e-3),
                              DiagoCGPrepare(20, 1000, 8, true, 1e-5, 300, 1e-3),
                              DiagoCGPrepare(40, 1000, 8, true, 1e-6, 300, 1e-3))); 
-                            // DiagoCGPrepare(40, 2000, 8, true, 1e-5, 500, 1e-2))); 
+                            //DiagoCGPrepare(40, 2000, 8, true, 1e-5, 500, 1e-2))); 
 			    // the last one is passed but time-consumming.
 
 // check that the mock class HPsi work well
@@ -187,6 +216,8 @@ TEST(DiagoCGTest, ZHEEV)
 }
 
 // cg for a 2x2 matrix
+#ifdef __MPI
+#else
 TEST(DiagoCGTest, TwoByTwo)
 {
     int dim = 2;
@@ -206,19 +237,20 @@ TEST(DiagoCGTest, TwoByTwo)
     DIAGOTEST::npw = dim;
     dcp.CompareEigen(hpsi.precond());
 }
+#endif
 
 TEST(DiagoCGTest, readH)
 {
     // read Hamilt matrix from file data-H
     ModuleBase::ComplexMatrix hm;
     std::ifstream ifs;
-    ifs.open("data-H");
+    ifs.open("H-KPoints-Si64.dat");
     DIAGOTEST::readh(ifs, hm);
     ifs.close();
     int dim = hm.nr;
     int nband = 10; // not nband < dim, here dim = 26 in data-H
     // nband, npw, sub, sparsity, reorder, eps, maxiter, threshold
-    DiagoCGPrepare dcp(nband, dim, 0, true, 1e-4, 300, 1e-3);
+    DiagoCGPrepare dcp(nband, dim, 0, true, 1e-5, 500, 1e-3);
     hsolver::DiagoIterAssist::PW_DIAG_NMAX = dcp.maxiter;
     hsolver::DiagoIterAssist::PW_DIAG_THR = dcp.eps;
     HPsi hpsi;
@@ -230,13 +262,28 @@ TEST(DiagoCGTest, readH)
 
 int main(int argc, char **argv)
 {
+	int nproc = 1, myrank = 0;
 
-    MPI_Init(&argc, &argv);
+#ifdef __MPI
+	int nproc_in_pool, kpar=1, mypool, rank_in_pool;
+    setupmpi(argc,argv,nproc, myrank);
+    divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool);
+    GlobalV::NPROC_IN_POOL = nproc;
+#else
+	MPI_Init(&argc, &argv);	
+#endif
 
     testing::InitGoogleTest(&argc, argv);
+    ::testing::TestEventListeners &listeners = ::testing::UnitTest::GetInstance()->listeners();
+    if (myrank != 0) delete listeners.Release(listeners.default_result_printer());
+
     int result = RUN_ALL_TESTS();
+    if (myrank == 0 && result != 0)
+    {
+        std::cout << "ERROR:some tests are not passed" << std::endl;
+        return result;
+	}
 
     MPI_Finalize();
-
-    return result;
+	return 0;
 }

source/module_hsolver/test/diago_cg_test.cpp

@@ -0,0 +1,13 @@
+#!/bin/bash
+
+np=`cat /proc/cpuinfo | grep "cpu cores" | uniq| awk '{print $NF}'`
+echo "nprocs in this machine is $np"
+
+for i in 6 3 2;do
+    if [[ $i -gt $np ]];then
+        continue
+    fi
+    echo "TEST DIAGO davidson in parallel, nprocs=$i"
+    mpirun -np $i ./HSolver_dav
+    break    
+done