Fix: use new genelpa interface (#1213)

* Add new GenELPA package, which support both new and legacy API of
ELPA, to ABACUS. Now ABACUS can use all versions of ELPA packages
instead of whose who are newer than 2017. Also, the program can
choose a better (if not the best) kernel to solve Kohn-Shan equation,
which makes it run (much) more faster than the one with the default
kernel.

More information about GenELPA: <https://github.com/pplab/GenELPA>

* Fix: only support elpa2016+

* Fix: remove legacy support in genelpa.

Co-authored-by: Shen Yu <[email protected]>
Co-authored-by: dyzheng <[email protected]>

Author: 3 people

CMakeLists.txt

@@ -336,6 +336,7 @@ target_link_libraries(${ABACUS_BIN_NAME}
     driver
     xc_
     hsolver
+    genelpa
     elecstate
     hamilt
     psi

modules/FindELPA.cmake

@@ -38,17 +38,4 @@ if(ELPA_FOUND)
 endif()
 
 set(CMAKE_REQUIRED_INCLUDES ${CMAKE_REQUIRED_INCLUDES} ${ELPA_INCLUDE_DIR})
-include(CheckCXXSourceCompiles)
-check_cxx_source_compiles("
-#include <elpa/elpa_version.h>
-#if ELPA_API_VERSION < 20210430
-#error ELPA version is too old.
-#endif
-int main(){}
-"
-ELPA_VERSION_SATISFIES
-)
-if(NOT ELPA_VERSION_SATISFIES)
-    message(FATAL_ERROR "ELPA version is too old. We support version 2017 or higher.")
-endif()
 mark_as_advanced(ELPA_INCLUDE_DIR ELPA_LIBRARY)

source/Makefile

@@ -20,6 +20,7 @@ VPATH=./src_global\
 :./module_xc\
 :./module_esolver\
 :./module_hsolver\
+:./module_hsolver/genelpa\
 :./module_elecstate\
 :./module_psi\
 :./module_hamilt\

source/Makefile.Objects

@@ -277,6 +277,11 @@ hsolver_lcao.o\
 hsolver_pw.o\
 hsolver_pw_sdft.o
 
+OBJ_GENELPA=elpa_new_complex.o\
+elpa_new_real.o\
+elpa_new.o\
+utils.o
+
 OBJ_ELECSTATES=elecstate.o\
 dm2d_to_grid.o\
 elecstate_lcao.o\
@@ -304,6 +309,7 @@ $(OBJ_HSOLVER)\
 $(OBJ_ELECSTATES)\
 $(OBJ_PSI)\
 ${OBJ_OPERATOR}\
+${OBJ_GENELPA}\
 charge.o \
 charge_mixing.o \
 charge_pulay.o \

source/module_deepks/test/CMakeLists.txt

@@ -8,7 +8,7 @@ target_link_libraries(
     test_deepks
     base cell symmetry md surchem xc_
     neighbor orb io relax gint lcao parallel mrrr pdiag pw ri driver esolver hsolver psi elecstate hamilt planewave
-    pthread
+    pthread genelpa
     deepks
     ${ABACUS_LINK_LIBRARIES}
 )

source/module_hamilt/hamilt_lcao.cpp

@@ -134,7 +134,7 @@ template <> void HamiltLCAO<double>::updateHk(const int ik)
         }
         const int inc = 1;
         BlasConnector::copy(this->LM->Sloc.size(), this->LM->Sloc.data(), inc, this->smatrix_k, inc);
-        hsolver::DiagoElpa::is_already_decomposed = false;
+        hsolver::DiagoElpa::DecomposedState = 0;
     }
     ModuleBase::timer::tick("HamiltLCAO", "updateHk");
     return;
@@ -303,4 +303,4 @@ template <> void HamiltLCAO<std::complex<double>>::constructHamilt()
 #endif
 }
 
-} // namespace hamilt
+} // namespace hamilt

source/module_hsolver/CMakeLists.txt

@@ -12,6 +12,8 @@ add_library(
     diago_lapack.cpp
 )
 
+add_subdirectory(genelpa)
+
 IF (BUILD_TESTING)
   add_subdirectory(test)
 endif()

source/module_hsolver/diago_elpa.cpp

@@ -7,52 +7,16 @@
 extern "C"
 {
 #include "module_base/blacs_connector.h"
-#include "my_elpa.h"
 #include "module_base/scalapack_connector.h"
 }
+#include "genelpa/elpa_solver.h"
 
 typedef hamilt::MatrixBlock<double> matd;
 typedef hamilt::MatrixBlock<std::complex<double>> matcd;
 
 namespace hsolver
 {
-bool DiagoElpa::is_already_decomposed = false;
-#ifdef __MPI
-inline int set_elpahandle(elpa_t &handle,
-                          const int *desc,
-                          const int local_nrows,
-                          const int local_ncols,
-                          const int nbands)
-{
-    int error;
-    int nprows, npcols, myprow, mypcol;
-    Cblacs_gridinfo(desc[1], &nprows, &npcols, &myprow, &mypcol);
-    elpa_init(20210430);
-    handle = elpa_allocate(&error);
-    elpa_set_integer(handle, "na", desc[2], &error);
-    elpa_set_integer(handle, "nev", nbands, &error);
-
-    elpa_set_integer(handle, "local_nrows", local_nrows, &error);
-
-    elpa_set_integer(handle, "local_ncols", local_ncols, &error);
-
-    elpa_set_integer(handle, "nblk", desc[4], &error);
-
-    elpa_set_integer(handle, "mpi_comm_parent", MPI_Comm_c2f(MPI_COMM_WORLD), &error);
-
-    elpa_set_integer(handle, "process_row", myprow, &error);
-
-    elpa_set_integer(handle, "process_col", mypcol, &error);
-
-    elpa_set_integer(handle, "blacs_context", desc[1], &error);
-
-    elpa_set_integer(handle, "cannon_for_generalized", 0, &error);
-    /* Setup */
-    elpa_setup(handle); /* Set tunables */
-    return 0;
-}
-#endif
-
+int DiagoElpa::DecomposedState = 0;
 void DiagoElpa::diag(hamilt::Hamilt *phm_in, psi::Psi<std::complex<double>> &psi, double *eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoElpa", "diag");
@@ -62,31 +26,15 @@ void DiagoElpa::diag(hamilt::Hamilt *phm_in, psi::Psi<std::complex<double>> &psi
 
     std::vector<double> eigen(GlobalV::NLOCAL, 0.0);
 
-    static elpa_t handle;
-    static bool has_set_elpa_handle = false;
-    if (!has_set_elpa_handle)
-    {
-        set_elpahandle(handle, h_mat.desc, h_mat.row, h_mat.col, GlobalV::NBANDS);
-        has_set_elpa_handle = true;
-    }
-
-    // compare to old code from pplab, there is no need to copy Sloc2 to another memory,
-    // just change Sloc2, which is a temporary matrix
-    // size_t nloc = h_mat.col * h_mat.row,
-    // BlasConnector::copy(nloc, s_mat, inc, Stmp, inc);
-
+    bool isReal=false;
+    const MPI_Comm COMM_DIAG=MPI_COMM_WORLD; // use all processes
+    ELPA_Solver es((const bool)isReal, COMM_DIAG, (const int)GlobalV::NBANDS, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc);
+    this->DecomposedState=0; // for k pointer, the decomposed s_mat can not be reused
     ModuleBase::timer::tick("DiagoElpa", "elpa_solve");
-    int elpa_derror;
-    elpa_generalized_eigenvectors_dc(handle,
-                                     reinterpret_cast<double _Complex *>(h_mat.p),
-                                     reinterpret_cast<double _Complex *>(s_mat.p),
-                                     eigen.data(),
-                                     reinterpret_cast<double _Complex *>(psi.get_pointer()),
-                                     0,
-                                     &elpa_derror);
+    es.generalized_eigenvector(h_mat.p, s_mat.p, this->DecomposedState, eigen.data(), psi.get_pointer());
     ModuleBase::timer::tick("DiagoElpa", "elpa_solve");
+    es.exit();
 
-    // the eigenvalues.
     const int inc = 1;
     BlasConnector::copy(GlobalV::NBANDS, eigen.data(), inc, eigenvalue_in, inc);
 #else
@@ -103,43 +51,14 @@ void DiagoElpa::diag(hamilt::Hamilt *phm_in, psi::Psi<double> &psi, double *eige
 
     std::vector<double> eigen(GlobalV::NLOCAL, 0.0);
 
-    static elpa_t handle;
-    static bool has_set_elpa_handle = false;
-    if (!has_set_elpa_handle)
-    {
-        set_elpahandle(handle, h_mat.desc, h_mat.row, h_mat.col, GlobalV::NBANDS);
-        has_set_elpa_handle = true;
-    }
-
-    // compare to old code from pplab, there is no need to copy Sloc2 to another memory,
-    // just change Sloc2, which is a temporary matrix
-    // change this judgement to HamiltLCAO
-    /*int is_already_decomposed;
-    if(ifElpaHandle(GlobalC::CHR.get_new_e_iteration(), (GlobalV::CALCULATION=="nscf")))
-    {
-        ModuleBase::timer::tick("DiagoElpa","decompose_S");
-        BlasConnector::copy(pv->nloc, s_mat, inc, Stmp, inc);
-        is_already_decomposed=0;
-        ModuleBase::timer::tick("DiagoElpa","decompose_S");
-    }
-    else
-    {
-        is_already_decomposed=1;
-    }*/
-
+    bool isReal=true;
+    MPI_Comm COMM_DIAG=MPI_COMM_WORLD; // use all processes
+    //ELPA_Solver es(isReal, COMM_DIAG, GlobalV::NBANDS, h_mat.row, h_mat.col, h_mat.desc);
+    ELPA_Solver es((const bool)isReal, COMM_DIAG, (const int)GlobalV::NBANDS, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc);
     ModuleBase::timer::tick("DiagoElpa", "elpa_solve");
-    int elpa_error;
-    elpa_generalized_eigenvectors_d(handle,
-                                    h_mat.p,
-                                    s_mat.p,
-                                    eigen.data(),
-                                    psi.get_pointer(),
-                                    DiagoElpa::is_already_decomposed,
-                                    &elpa_error);
+    es.generalized_eigenvector(h_mat.p, s_mat.p, this->DecomposedState, eigen.data(), psi.get_pointer());
     ModuleBase::timer::tick("DiagoElpa", "elpa_solve");
-
-    //S matrix has been decomposed
-    DiagoElpa::is_already_decomposed = true;
+    es.exit();
 
     const int inc = 1;
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "K-S equation was solved by genelpa2");
@@ -162,4 +81,4 @@ bool DiagoElpa::ifElpaHandle(const bool &newIteration, const bool &ifNSCF)
 }
 #endif
 
-} // namespace hsolver
+} // namespace hsolver

source/module_hsolver/diago_elpa.h

@@ -14,9 +14,9 @@ class DiagoElpa : public DiagH
     void diag(hamilt::Hamilt* phm_in, psi::Psi<double>& psi, double* eigenvalue_in) override;
 
     void diag(hamilt::Hamilt* phm_in, psi::Psi<std::complex<double>>& psi, double* eigenvalue_in) override;
-
-    static bool is_already_decomposed;
 
+    static int DecomposedState;
+
   private:
 #ifdef __MPI
     bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF);
@@ -25,4 +25,4 @@ class DiagoElpa : public DiagH
 
 } // namespace hsolver
 
-#endif
+#endif

source/module_hsolver/genelpa/CMakeLists.txt

@@ -0,0 +1 @@
+add_library(genelpa OBJECT elpa_new.cpp elpa_new_real.cpp elpa_new_complex.cpp utils.cpp)