diff --git a/source/source_hsolver/module_genelpa/cblacs.h b/source/source_hsolver/module_genelpa/cblacs.h index 3b671ef1f9f..db282187727 100644 --- a/source/source_hsolver/module_genelpa/cblacs.h +++ b/source/source_hsolver/module_genelpa/cblacs.h @@ -2,28 +2,29 @@ #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_CBLACS_H // blacs - // Initialization +// Initialization #include "mpi.h" + #include int Csys2blacs_handle(MPI_Comm SysCtxt); -void Cblacs_pinfo(int *myid, int *nprocs); -void Cblacs_get(int icontxt, int what, int *val); -void Cblacs_gridinit(int* icontxt, char *layout, int nprow, int npcol); -void Cblacs_gridmap(int* icontxt, int *usermap, int ldumap, int nprow, int npcol); - // Destruction +void Cblacs_pinfo(int* myid, int* nprocs); +void Cblacs_get(int icontxt, int what, int* val); +void Cblacs_gridinit(int* icontxt, char* layout, int nprow, int npcol); +void Cblacs_gridmap(int* icontxt, int* usermap, int ldumap, int nprow, int npcol); +// Destruction void Cblacs_gridexit(int icontxt); - // Informational and Miscellaneous -void Cblacs_gridinfo(int icontxt, int* nprow, int *npcol, int *myprow, int *mypcol); +// Informational and Miscellaneous +void Cblacs_gridinfo(int icontxt, int* nprow, int* npcol, int* myprow, int* mypcol); int Cblacs_pnum(int icontxt, int prow, int pcol); -void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol); -void Cblacs_barrier(int icontxt, char *scope); - // Point to Point -void Cdgesd2d(int icontxt, int m, int n, double *a, int lda, int rdest, int cdest); -void Cdgerv2d(int icontxt, int m, int n, double *a, int lda, int rsrc, int csrc); -void Czgesd2d(int icontxt, int m, int n, std::complex *a, int lda, int rdest, int cdest); -void Czgerv2d(int icontxt, int m, int n, std::complex *a, int lda, int rsrc, int csrc); - // Combine -//void Cdgamx2d(int icontxt, int scope, int top, int m, int n, -// double *a, int lda, int *ra, int *ca, int rcflag, int rdest, int cdest); +void Cblacs_pcoord(int icontxt, int pnum, int* prow, int* pcol); +void Cblacs_barrier(int icontxt, char* scope); +// Point to Point +void Cdgesd2d(int icontxt, int m, int n, double* a, int lda, int rdest, int cdest); +void Cdgerv2d(int icontxt, int m, int n, double* a, int lda, int rsrc, int csrc); +void Czgesd2d(int icontxt, int m, int n, std::complex* a, int lda, int rdest, int cdest); +void Czgerv2d(int icontxt, int m, int n, std::complex* a, int lda, int rsrc, int csrc); +// Combine +// void Cdgamx2d(int icontxt, int scope, int top, int m, int n, +// double *a, int lda, int *ra, int *ca, int rcflag, int rdest, int cdest); #endif // ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_CBLACS_H diff --git a/source/source_hsolver/module_genelpa/elpa_generic.h b/source/source_hsolver/module_genelpa/elpa_generic.h index 1638dee751b..6945c50fe3b 100644 --- a/source/source_hsolver/module_genelpa/elpa_generic.h +++ b/source/source_hsolver/module_genelpa/elpa_generic.h @@ -2,6 +2,7 @@ #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_ELPA_GENERIC_H #include "elpa_new.h" + #include /*! \brief generic C method for elpa_set * @@ -12,11 +13,11 @@ * \param error on return the error code, which can be queried with elpa_strerr() * \result void */ -inline void elpa_set(elpa_t handle, const char *name, int value, int *error) +inline void elpa_set(elpa_t handle, const char* name, int value, int* error) { elpa_set_integer(handle, name, value, error); } -inline void elpa_set(elpa_t handle, const char *name, double value, int *error) +inline void elpa_set(elpa_t handle, const char* name, double value, int* error) { elpa_set_double(handle, name, value, error); } @@ -30,11 +31,11 @@ inline void elpa_set(elpa_t handle, const char *name, double value, int *error) * \param error on return the error code, which can be queried with elpa_strerr() * \result void */ -inline void elpa_get(elpa_t handle, const char *name, int *value, int *error) +inline void elpa_get(elpa_t handle, const char* name, int* value, int* error) { elpa_get_integer(handle, name, value, error); } -inline void elpa_get(elpa_t handle, const char *name, double *value, int *error) +inline void elpa_get(elpa_t handle, const char* name, double* value, int* error) { elpa_get_double(handle, name, value, error); } @@ -50,68 +51,108 @@ inline void elpa_get(elpa_t handle, const char *name, double *value, int *error) * \result void */ #if ELPA_API_VERSION <= 20210502 // ELPA 2021.05.002 and earlier versions -inline void elpa_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_eigenvectors_d(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_eigenvectors_f(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, double *ev, std::complex *q, int *error) -{ - elpa_eigenvectors_dc(handle, reinterpret_cast(a), ev, reinterpret_cast(q), error); +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + double* ev, + std::complex* q, + int* error) +{ + elpa_eigenvectors_dc(handle, + reinterpret_cast(a), + ev, + reinterpret_cast(q), + error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, float *ev, std::complex *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + float* ev, + std::complex* q, + int* error) { elpa_eigenvectors_fc(handle, reinterpret_cast(a), ev, reinterpret_cast(q), error); } -#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 -inline void elpa_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 +inline void elpa_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_eigenvectors_all_host_arrays_d(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_eigenvectors_all_host_arrays_f(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, double *ev, std::complex *q, int *error) -{ - elpa_eigenvectors_all_host_arrays_dc(handle, reinterpret_cast(a), - ev, reinterpret_cast(q), error); +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + double* ev, + std::complex* q, + int* error) +{ + elpa_eigenvectors_all_host_arrays_dc(handle, + reinterpret_cast(a), + ev, + reinterpret_cast(q), + error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, float *ev, std::complex *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + float* ev, + std::complex* q, + int* error) { - elpa_eigenvectors_all_host_arrays_fc(handle, reinterpret_cast(a), - ev, reinterpret_cast(q), error); + elpa_eigenvectors_all_host_arrays_fc(handle, + reinterpret_cast(a), + ev, + reinterpret_cast(q), + error); } -#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 -inline void elpa_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 +inline void elpa_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_eigenvectors_a_h_a_d(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_eigenvectors_a_h_a_f(handle, a, ev, q, error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, double *ev, std::complex *q, int *error) -{ - elpa_eigenvectors_a_h_a_dc(handle, reinterpret_cast(a), - ev, reinterpret_cast(q), error); +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + double* ev, + std::complex* q, + int* error) +{ + elpa_eigenvectors_a_h_a_dc(handle, + reinterpret_cast(a), + ev, + reinterpret_cast(q), + error); } -inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, float *ev, std::complex *q, int *error) -{ - elpa_eigenvectors_a_h_a_fc(handle, reinterpret_cast(a), - ev, reinterpret_cast(q), error); +inline void elpa_eigenvectors(const elpa_t handle, + std::complex* a, + float* ev, + std::complex* q, + int* error) +{ + elpa_eigenvectors_a_h_a_fc(handle, + reinterpret_cast(a), + ev, + reinterpret_cast(q), + error); } #endif @@ -126,39 +167,37 @@ inline void elpa_eigenvectors(const elpa_t handle, std::complex *a, floa * \result void */ #if ELPA_API_VERSION <= 20210502 // ELPA 2021.05.002 and earlier versions -inline void elpa_skew_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +inline void elpa_skew_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_eigenvectors_d(handle, a, ev, q, error); } -inline void elpa_skew_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_skew_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_eigenvectors_f(handle, a, ev, q, error); } -#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 -inline void elpa_skew_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 +inline void elpa_skew_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_skew_eigenvectors_all_host_arrays_d(handle, a, ev, q, error); } -inline void elpa_skew_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_skew_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_skew_eigenvectors_all_host_arrays_f(handle, a, ev, q, error); } -#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 -inline void elpa_skew_eigenvectors(const elpa_t handle, double *a, double *ev, double *q, int *error) +#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 +inline void elpa_skew_eigenvectors(const elpa_t handle, double* a, double* ev, double* q, int* error) { elpa_skew_eigenvectors_a_h_a_d(handle, a, ev, q, error); } -inline void elpa_skew_eigenvectors(const elpa_t handle, float *a, float *ev, float *q, int *error) +inline void elpa_skew_eigenvectors(const elpa_t handle, float* a, float* ev, float* q, int* error) { elpa_skew_eigenvectors_a_h_a_f(handle, a, ev, q, error); } #endif - - /*! \brief generic C method for elpa_generalized_eigenvectors * * \details @@ -171,26 +210,60 @@ inline void elpa_skew_eigenvectors(const elpa_t handle, float *a, float *ev, f * \param error on return the error code, which can be queried with elpa_strerr() * \result void */ -inline void elpa_generalized_eigenvectors(elpa_t handle, double *a, double *b, double *ev, double *q, int is_already_decomposed, int *error) +inline void elpa_generalized_eigenvectors(elpa_t handle, + double* a, + double* b, + double* ev, + double* q, + int is_already_decomposed, + int* error) { elpa_generalized_eigenvectors_d(handle, a, b, ev, q, is_already_decomposed, error); } -inline void elpa_generalized_eigenvectors(elpa_t handle, float *a, float *b, float *ev, float *q, int is_already_decomposed, int *error) +inline void elpa_generalized_eigenvectors(elpa_t handle, + float* a, + float* b, + float* ev, + float* q, + int is_already_decomposed, + int* error) { elpa_generalized_eigenvectors_f(handle, a, b, ev, q, is_already_decomposed, error); } -inline void elpa_generalized_eigenvectors(elpa_t handle, std::complex *a, std::complex *b, double *ev, std::complex *q, int is_already_decomposed, int *error) -{ - elpa_generalized_eigenvectors_dc(handle, reinterpret_cast(a), reinterpret_cast(b), - ev, reinterpret_cast(q), is_already_decomposed, error); +inline void elpa_generalized_eigenvectors(elpa_t handle, + std::complex* a, + std::complex* b, + double* ev, + std::complex* q, + int is_already_decomposed, + int* error) +{ + elpa_generalized_eigenvectors_dc(handle, + reinterpret_cast(a), + reinterpret_cast(b), + ev, + reinterpret_cast(q), + is_already_decomposed, + error); } -inline void elpa_generalized_eigenvectors(elpa_t handle, std::complex *a, std::complex *b, float *ev, std::complex *q, int is_already_decomposed, int *error) -{ - elpa_generalized_eigenvectors_fc(handle, reinterpret_cast(a), reinterpret_cast(b), - ev, reinterpret_cast(q), is_already_decomposed, error); +inline void elpa_generalized_eigenvectors(elpa_t handle, + std::complex* a, + std::complex* b, + float* ev, + std::complex* q, + int is_already_decomposed, + int* error) +{ + elpa_generalized_eigenvectors_fc(handle, + reinterpret_cast(a), + reinterpret_cast(b), + ev, + reinterpret_cast(q), + is_already_decomposed, + error); } /*! \brief generic C method for elpa_eigenvalues @@ -203,53 +276,53 @@ inline void elpa_generalized_eigenvectors(elpa_t handle, std::complex *a * \result void */ #if ELPA_API_VERSION <= 20210502 // ELPA 2021.05.002 and earlier versions -inline void elpa_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_d(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_f(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, double *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, double* ev, int* error) { elpa_eigenvalues_dc(handle, reinterpret_cast(a), ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, float* ev, int* error) { - elpa_eigenvalues_fc (handle, reinterpret_cast(a), ev, error); + elpa_eigenvalues_fc(handle, reinterpret_cast(a), ev, error); } -#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 -inline void elpa_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 +inline void elpa_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_all_host_arrays_d(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_all_host_arrays_f(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, double *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, double* ev, int* error) { elpa_eigenvalues_all_host_arrays_dc(handle, reinterpret_cast(a), ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, float* ev, int* error) { elpa_eigenvalues_all_host_arrays_fc(handle, reinterpret_cast(a), ev, error); } -#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 -inline void elpa_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 +inline void elpa_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_a_h_a_d(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_a_h_a_f(handle, a, ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, double *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, double* ev, int* error) { elpa_eigenvalues_a_h_a_dc(handle, reinterpret_cast(a), ev, error); } -inline void elpa_eigenvalues(elpa_t handle, std::complex *a, float *ev, int *error) +inline void elpa_eigenvalues(elpa_t handle, std::complex* a, float* ev, int* error) { elpa_eigenvalues_a_h_a_fc(handle, reinterpret_cast(a), ev, error); } @@ -265,29 +338,29 @@ inline void elpa_eigenvalues(elpa_t handle, std::complex *a, float *ev, * \result void */ #if ELPA_API_VERSION <= 20210502 // ELPA 2021.05.002 and earlier versions -inline void elpa_skew_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +inline void elpa_skew_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_d(handle, a, ev, error); } -inline void elpa_skew_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_skew_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_f(handle, a, ev, error); } -#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 -inline void elpa_skew_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +#elif ELPA_API_VERSION < 20220501 // ELPA version between 2021.11.001 and 2022.05.001 +inline void elpa_skew_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_all_host_arrays_d(handle, a, ev, error); } -inline void elpa_skew_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_skew_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_all_host_arrays_f(handle, a, ev, error); } -#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 -inline void elpa_skew_eigenvalues(elpa_t handle, double *a, double *ev, int *error) +#else // ELPA version 2022.05.001, ELPA has its own c++ interface from version 2022.11.001 +inline void elpa_skew_eigenvalues(elpa_t handle, double* a, double* ev, int* error) { elpa_eigenvalues_a_h_a_d(handle, a, ev, error); } -inline void elpa_skew_eigenvalues(elpa_t handle, float *a, float *ev, int *error) +inline void elpa_skew_eigenvalues(elpa_t handle, float* a, float* ev, int* error) { elpa_eigenvalues_a_h_a_f(handle, a, ev, error); } @@ -303,37 +376,37 @@ inline void elpa_skew_eigenvalues(elpa_t handle, float *a, float *ev, int *err * \result void */ -#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 -inline void elpa_cholesky(elpa_t handle, double *a, int *error) +#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 +inline void elpa_cholesky(elpa_t handle, double* a, int* error) { elpa_cholesky_d(handle, a, error); } -inline void elpa_cholesky(elpa_t handle, float *a, int *error) +inline void elpa_cholesky(elpa_t handle, float* a, int* error) { elpa_cholesky_f(handle, a, error); } -inline void elpa_cholesky(elpa_t handle, std::complex *a, int *error) +inline void elpa_cholesky(elpa_t handle, std::complex* a, int* error) { elpa_cholesky_dc(handle, reinterpret_cast(a), error); } -inline void elpa_cholesky(elpa_t handle, std::complex *a, int *error) +inline void elpa_cholesky(elpa_t handle, std::complex* a, int* error) { elpa_cholesky_fc(handle, reinterpret_cast(a), error); } #else -inline void elpa_cholesky(elpa_t handle, double *a, int *error) +inline void elpa_cholesky(elpa_t handle, double* a, int* error) { elpa_cholesky_a_h_a_d(handle, a, error); } -inline void elpa_cholesky(elpa_t handle, float *a, int *error) +inline void elpa_cholesky(elpa_t handle, float* a, int* error) { elpa_cholesky_a_h_a_f(handle, a, error); } -inline void elpa_cholesky(elpa_t handle, std::complex *a, int *error) +inline void elpa_cholesky(elpa_t handle, std::complex* a, int* error) { elpa_cholesky_a_h_a_dc(handle, reinterpret_cast(a), error); } -inline void elpa_cholesky(elpa_t handle, std::complex *a, int *error) +inline void elpa_cholesky(elpa_t handle, std::complex* a, int* error) { elpa_cholesky_a_h_a_fc(handle, reinterpret_cast(a), error); } @@ -356,47 +429,171 @@ inline void elpa_cholesky(elpa_t handle, std::complex *a, int *error) * \param error on return the error code, which can be queried with elpa_strerr() * \result void */ -#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, double *a, double *b, int nrows_b, int ncols_b, double *c, int nrows_c, int ncols_c, int *error) +#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + double* a, + double* b, + int nrows_b, + int ncols_b, + double* c, + int nrows_c, + int ncols_c, + int* error) { elpa_hermitian_multiply_d(handle, uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, c, nrows_c, ncols_c, error); } -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, float *a, float *b, int nrows_b, int ncols_b, float *c, int nrows_c, int ncols_c, int *error) +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + float* a, + float* b, + int nrows_b, + int ncols_b, + float* c, + int nrows_c, + int ncols_c, + int* error) { elpa_hermitian_multiply_df(handle, uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, c, nrows_c, ncols_c, error); } -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, std::complex *a, std::complex *b, int nrows_b, int ncols_b, std::complex *c, int nrows_c, int ncols_c, int *error) -{ - elpa_hermitian_multiply_dc(handle, uplo_a, uplo_c, ncb, reinterpret_cast(a), - reinterpret_cast(b), nrows_b, ncols_b, - reinterpret_cast(c), nrows_c, ncols_c, error); -} -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, std::complex *a, std::complex *b, int nrows_b, int ncols_b, std::complex *c, int nrows_c, int ncols_c, int *error) -{ - elpa_hermitian_multiply_fc(handle, uplo_a, uplo_c, ncb, reinterpret_cast(a), - reinterpret_cast(b), nrows_b, ncols_b, - reinterpret_cast(c), nrows_c, ncols_c, error); +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + std::complex* a, + std::complex* b, + int nrows_b, + int ncols_b, + std::complex* c, + int nrows_c, + int ncols_c, + int* error) +{ + elpa_hermitian_multiply_dc(handle, + uplo_a, + uplo_c, + ncb, + reinterpret_cast(a), + reinterpret_cast(b), + nrows_b, + ncols_b, + reinterpret_cast(c), + nrows_c, + ncols_c, + error); +} +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + std::complex* a, + std::complex* b, + int nrows_b, + int ncols_b, + std::complex* c, + int nrows_c, + int ncols_c, + int* error) +{ + elpa_hermitian_multiply_fc(handle, + uplo_a, + uplo_c, + ncb, + reinterpret_cast(a), + reinterpret_cast(b), + nrows_b, + ncols_b, + reinterpret_cast(c), + nrows_c, + ncols_c, + error); } #else -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, double *a, double *b, int nrows_b, int ncols_b, double *c, int nrows_c, int ncols_c, int *error) +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + double* a, + double* b, + int nrows_b, + int ncols_b, + double* c, + int nrows_c, + int ncols_c, + int* error) { elpa_hermitian_multiply_a_h_a_d(handle, uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, c, nrows_c, ncols_c, error); } -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, float *a, float *b, int nrows_b, int ncols_b, float *c, int nrows_c, int ncols_c, int *error) +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + float* a, + float* b, + int nrows_b, + int ncols_b, + float* c, + int nrows_c, + int ncols_c, + int* error) { elpa_hermitian_multiply_a_h_a_f(handle, uplo_a, uplo_c, ncb, a, b, nrows_b, ncols_b, c, nrows_c, ncols_c, error); } -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, std::complex *a, std::complex *b, int nrows_b, int ncols_b, std::complex *c, int nrows_c, int ncols_c, int *error) -{ - elpa_hermitian_multiply_a_h_a_dc(handle, uplo_a, uplo_c, ncb, reinterpret_cast(a), - reinterpret_cast(b), nrows_b, ncols_b, - reinterpret_cast(c), nrows_c, ncols_c, error); -} -inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int ncb, std::complex *a, std::complex *b, int nrows_b, int ncols_b, std::complex *c, int nrows_c, int ncols_c, int *error) -{ - elpa_hermitian_multiply_a_h_a_fc(handle, uplo_a, uplo_c, ncb, reinterpret_cast(a), - reinterpret_cast(b), nrows_b, ncols_b, - reinterpret_cast(c), nrows_c, ncols_c, error); +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + std::complex* a, + std::complex* b, + int nrows_b, + int ncols_b, + std::complex* c, + int nrows_c, + int ncols_c, + int* error) +{ + elpa_hermitian_multiply_a_h_a_dc(handle, + uplo_a, + uplo_c, + ncb, + reinterpret_cast(a), + reinterpret_cast(b), + nrows_b, + ncols_b, + reinterpret_cast(c), + nrows_c, + ncols_c, + error); +} +inline void elpa_hermitian_multiply(elpa_t handle, + char uplo_a, + char uplo_c, + int ncb, + std::complex* a, + std::complex* b, + int nrows_b, + int ncols_b, + std::complex* c, + int nrows_c, + int ncols_c, + int* error) +{ + elpa_hermitian_multiply_a_h_a_fc(handle, + uplo_a, + uplo_c, + ncb, + reinterpret_cast(a), + reinterpret_cast(b), + nrows_b, + ncols_b, + reinterpret_cast(c), + nrows_c, + ncols_c, + error); } #endif @@ -409,37 +606,37 @@ inline void elpa_hermitian_multiply(elpa_t handle, char uplo_a, char uplo_c, int * \param error on return the error code, which can be queried with elpa_strerr() * \result void */ -#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 -inline void elpa_invert_triangular(elpa_t handle, double *a, int *error) +#if ELPA_API_VERSION < 20220501 // ELPA version before 2022.05.001 +inline void elpa_invert_triangular(elpa_t handle, double* a, int* error) { elpa_invert_trm_d(handle, a, error); } -inline void elpa_invert_triangular(elpa_t handle, float *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, float* a, int* error) { elpa_invert_trm_f(handle, a, error); } -inline void elpa_invert_triangular(elpa_t handle, std::complex *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, std::complex* a, int* error) { elpa_invert_trm_dc(handle, reinterpret_cast(a), error); } -inline void elpa_invert_triangular(elpa_t handle, std::complex *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, std::complex* a, int* error) { elpa_invert_trm_fc(handle, reinterpret_cast(a), error); } #else -inline void elpa_invert_triangular(elpa_t handle, double *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, double* a, int* error) { elpa_invert_trm_a_h_a_d(handle, a, error); } -inline void elpa_invert_triangular(elpa_t handle, float *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, float* a, int* error) { elpa_invert_trm_a_h_a_f(handle, a, error); } -inline void elpa_invert_triangular(elpa_t handle, std::complex *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, std::complex* a, int* error) { elpa_invert_trm_a_h_a_dc(handle, reinterpret_cast(a), error); } -inline void elpa_invert_triangular(elpa_t handle, std::complex *a, int *error) +inline void elpa_invert_triangular(elpa_t handle, std::complex* a, int* error) { elpa_invert_trm_a_h_a_fc(handle, reinterpret_cast(a), error); } diff --git a/source/source_hsolver/module_genelpa/elpa_new.cpp b/source/source_hsolver/module_genelpa/elpa_new.cpp index fa7bb917e4e..9d6cc5f6140 100644 --- a/source/source_hsolver/module_genelpa/elpa_new.cpp +++ b/source/source_hsolver/module_genelpa/elpa_new.cpp @@ -1,13 +1,13 @@ #include "elpa_new.h" #include "elpa_runtime_check.h" - #include "elpa_solver.h" extern "C" { #include "cblacs.h" } #include "utils.h" + #include #include #include @@ -73,8 +73,8 @@ ELPA_Solver::ELPA_Solver(const bool isReal, handle_id = ++total_handle; - //delete the old elpa_handle and reuse the handle_id=0 - if(reuse_handle_0 && total_handle>0) + // delete the old elpa_handle and reuse the handle_id=0 + if (reuse_handle_0 && total_handle > 0) { NEW_ELPA_HANDLE_POOL.erase(0); handle_id = 0; @@ -463,7 +463,7 @@ int ELPA_Solver::read_complex_kernel() int ELPA_Solver::allocate_work() { unsigned long nloc = static_cast(narows) * nacols; // local size - unsigned long maxloc = 0; // maximum local size + unsigned long maxloc = 0; // maximum local size MPI_Allreduce(&nloc, &maxloc, 1, MPI_UNSIGNED_LONG, MPI_MAX, comm); maxloc = nloc; @@ -480,13 +480,13 @@ void ELPA_Solver::timer(int myid, const char function[], const char step[], doub if (t0 < 0) // t0 < 0 means this is the init call before the function { t0 = MPI_Wtime(); - t0 = (double)clock()/CLOCKS_PER_SEC; + t0 = (double)clock() / CLOCKS_PER_SEC; logfile << "DEBUG: Process " << myid << " Call " << function << std::endl; } else { t1 = MPI_Wtime(); - t1 = (double)clock()/CLOCKS_PER_SEC; + t1 = (double)clock() / CLOCKS_PER_SEC; logfile << "DEBUG: Process " << myid << " Step " << step << " " << function << " time: " << t1 - t0 << " s" << std::endl; } diff --git a/source/source_hsolver/module_genelpa/elpa_new.h b/source/source_hsolver/module_genelpa/elpa_new.h index df36ad37a52..8f6b761dd56 100644 --- a/source/source_hsolver/module_genelpa/elpa_new.h +++ b/source/source_hsolver/module_genelpa/elpa_new.h @@ -1,11 +1,8 @@ #ifndef ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_ELPA_NEW_H #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_ELPA_NEW_H - #include - - #if ELPA_API_VERSION >= 20221101 #include #else @@ -14,10 +11,10 @@ extern "C" #include struct elpa_struct; - typedef struct elpa_struct *elpa_t; + typedef struct elpa_struct* elpa_t; struct elpa_autotune_struct; - typedef struct elpa_autotune_struct *elpa_autotune_t; + typedef struct elpa_autotune_struct* elpa_autotune_t; #include #include @@ -30,7 +27,7 @@ extern "C" // #include #undef complex - const char *elpa_strerr(int elpa_error); + const char* elpa_strerr(int elpa_error); // This is a header file related to GPU support. // This header file was not available in the early ELPA. diff --git a/source/source_hsolver/module_genelpa/elpa_new_complex.cpp b/source/source_hsolver/module_genelpa/elpa_new_complex.cpp index 9aec04fa72b..4caf4d69a6a 100644 --- a/source/source_hsolver/module_genelpa/elpa_new_complex.cpp +++ b/source/source_hsolver/module_genelpa/elpa_new_complex.cpp @@ -1,19 +1,16 @@ -#include -#include -#include -#include -#include -#include -#include - #include "elpa_new.h" #include "elpa_solver.h" - #include "source_base/module_external/scalapack_connector.h" -#include "utils.h" - #include "source_base/tool_quit.h" +#include "utils.h" +#include +#include +#include +#include +#include +#include +#include extern std::map NEW_ELPA_HANDLE_POOL; @@ -23,13 +20,13 @@ int ELPA_Solver::eigenvector(std::complex* A, double* EigenValue, std::c int allinfo = 0; double t = 0.0; - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { - t=-1; + t = -1; timer(myid, "elpa_eigenvectors_dc", "1", t); } elpa_eigenvectors(NEW_ELPA_HANDLE_POOL[handle_id], A, EigenValue, EigenVector, &info); - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { timer(myid, "elpa_eigenvectors_dc", "1", t); } @@ -38,62 +35,69 @@ int ELPA_Solver::eigenvector(std::complex* A, double* EigenValue, std::c return allinfo; } -int ELPA_Solver::generalized_eigenvector(std::complex* A, std::complex* B, int& DecomposedState, - double* EigenValue, std::complex* EigenVector) +int ELPA_Solver::generalized_eigenvector(std::complex* A, + std::complex* B, + int& DecomposedState, + double* EigenValue, + std::complex* EigenVector) { int info, allinfo; double t = 0.0; - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { - t=-1; + t = -1; timer(myid, "decomposeRightMatrix", "1", t); } - if(DecomposedState==0) // B is not decomposed - allinfo=decomposeRightMatrix(B, EigenValue, EigenVector, DecomposedState); + if (DecomposedState == 0) // B is not decomposed + allinfo = decomposeRightMatrix(B, EigenValue, EigenVector, DecomposedState); else - allinfo=0; + allinfo = 0; - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { timer(myid, "decomposeRightMatrix", "1", t); } - if(allinfo != 0){ + if (allinfo != 0) + { // if allinfo is still not 0 anyway, report error and quit - if(myid == 0){ - ModuleBase::WARNING_QUIT("ELPA_Solver::generalized_eigenvector", "decomposeRightMatrix failed to decompose right matrix!\n info = " + std::to_string(allinfo)); + if (myid == 0) + { + ModuleBase::WARNING_QUIT("ELPA_Solver::generalized_eigenvector", + "decomposeRightMatrix failed to decompose right matrix!\n info = " + + std::to_string(allinfo)); } } // transform A to A~ - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { - t=-1; + t = -1; timer(myid, "transform A to A~", "2", t); } - if(DecomposedState == 1 || DecomposedState == 2) + if (DecomposedState == 1 || DecomposedState == 2) { // Form H * U^{-1} in work. Only the upper triangle of H is // authoritative, so use the Hermitian matrix-matrix product. - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "A*U^-1", "2.1a", t); } ScalapackConnector::hemm('L', 'U', nFull, 1.0, A, B, 0.0, zwork.data(), desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "A*U^-1", "2.1a", t); } // Form U^{-H} * (H * U^{-1}) in A. - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "U^-H*(H*U^-1)", "2.2a", t); } ScalapackConnector::gemm('C', 'N', nFull, nFull, nFull, 1.0, B, zwork.data(), 0.0, A, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "U^-H*(H*U^-1)", "2.2a", t); } @@ -102,42 +106,43 @@ int ELPA_Solver::generalized_eigenvector(std::complex* A, std::complex1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "B*H", "2.1b", t); } ScalapackConnector::hemm('R', 'U', nFull, 1.0, A, B, 0.0, zwork.data(), desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "B*H", "2.1b", t); } // Form B * work^H = B * H * B^H in A. - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "B*H*B^H", "2.2b", t); } ScalapackConnector::gemm('N', 'C', nFull, nFull, nFull, 1.0, B, zwork.data(), 0.0, A, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "B*H*B^H", "2.2b", t); } } - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { timer(myid, "transform A to A~", "2", t); } // calculate the eigenvalue and eigenvector of A~ - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { - t=-1; + t = -1; timer(myid, "elpa_eigenvectors", "3", t); } - if(loglevel>2) saveMatrix("A_tilde.dat", nFull, A, desc, cblacs_ctxt); - info=eigenvector(A, EigenValue, EigenVector); - if((loglevel>0 && myid==0) || loglevel>1) + if (loglevel > 2) + saveMatrix("A_tilde.dat", nFull, A, desc, cblacs_ctxt); + info = eigenvector(A, EigenValue, EigenVector); + if ((loglevel > 0 && myid == 0) || loglevel > 1) { timer(myid, "elpa_eigenvectors", "3", t); } @@ -146,85 +151,90 @@ int ELPA_Solver::generalized_eigenvector(std::complex* A, std::complex2) saveMatrix("EigenVector_tilde.dat", nFull, EigenVector, desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("EigenVector_tilde.dat", nFull, EigenVector, desc, cblacs_ctxt); - if((loglevel>0 && myid==0) || loglevel>1) + if ((loglevel > 0 && myid == 0) || loglevel > 1) { - t=-1; + t = -1; timer(myid, "composeEigenVector", "4", t); } // transform eigenvector c~ to original eigenvector c - allinfo=composeEigenVector(DecomposedState, B, EigenVector); - if((loglevel>0 && myid==0) || loglevel>1) + allinfo = composeEigenVector(DecomposedState, B, EigenVector); + if ((loglevel > 0 && myid == 0) || loglevel > 1) { timer(myid, "composeEigenVector", "4", t); } return allinfo; } -int ELPA_Solver::decomposeRightMatrix(std::complex* B, double* EigenValue, std::complex* EigenVector, int& DecomposedState) +int ELPA_Solver::decomposeRightMatrix(std::complex* B, + double* EigenValue, + std::complex* EigenVector, + int& DecomposedState) { - int info=0; - int allinfo=0; + int info = 0; + int allinfo = 0; double t = 0.0; // first try cholesky decomposing - if(nFull1) + DecomposedState = 1; + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "pzpotrf_", "1", t); } ScalapackConnector::potrf('U', nFull, B, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "pzpotrf_", "1", t); } MPI_Allreduce(&info, &allinfo, 1, MPI_INT, MPI_MAX, comm); allinfo = info; - if(allinfo != 0) //if pdpotrf fail, try elpa_cholesky_real + if (allinfo != 0) // if pdpotrf fail, try elpa_cholesky_real { - DecomposedState=2; - if(loglevel>1) + DecomposedState = 2; + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "elpa_cholesky_dc", "2", t); } elpa_cholesky(NEW_ELPA_HANDLE_POOL[handle_id], B, &info); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "elpa_cholesky_dc", "2", t); } MPI_Allreduce(&info, &allinfo, 1, MPI_INT, MPI_MAX, comm); allinfo = info; } - } else + } + else { - DecomposedState=2; - if(loglevel>1) + DecomposedState = 2; + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "elpa_cholesky_dc", "1", t); } elpa_cholesky(NEW_ELPA_HANDLE_POOL[handle_id], B, &info); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "elpa_cholesky_dc", "1", t); } MPI_Allreduce(&info, &allinfo, 1, MPI_INT, MPI_MAX, comm); allinfo = info; - if(allinfo != 0) + if (allinfo != 0) { - DecomposedState=1; - if(loglevel>1) + DecomposedState = 1; + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "pzpotrf_", "2", t); } ScalapackConnector::potrf('U', nFull, B, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "pzpotrf_", "2", t); } @@ -233,84 +243,90 @@ int ELPA_Solver::decomposeRightMatrix(std::complex* B, double* EigenValu } } - if(allinfo==0) // calculate U^{-1} + if (allinfo == 0) // calculate U^{-1} { - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "clear low triangle", "1", t); } - for(int j=0; jjGlobal) B[i+j*narows]=0; + int iGlobal = globalIndex(i, nblk, nprows, myprow); + if (iGlobal > jGlobal) + B[i + j * narows] = 0; } } - if(loglevel>1) + if (loglevel > 1) { timer(myid, "clear low triangle", "1", t); } - if(loglevel>2) saveMatrix("U.dat", nFull, B, desc, cblacs_ctxt); - if(loglevel>1) + if (loglevel > 2) + saveMatrix("U.dat", nFull, B, desc, cblacs_ctxt); + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "invert U", "1", t); } elpa_invert_triangular(NEW_ELPA_HANDLE_POOL[handle_id], B, &info); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "invert U", "1", t); } - if(loglevel>2) saveMatrix("U_inv.dat", nFull, B, desc, cblacs_ctxt); - } else { + if (loglevel > 2) + saveMatrix("U_inv.dat", nFull, B, desc, cblacs_ctxt); + } + else + { // if cholesky decomposing failed, try diagonalize // Form B^{-1/2} = Q * diag(ev^{-1/2}) * Q^H in B. - DecomposedState=3; - if(loglevel>1) + DecomposedState = 3; + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "calculate eigenvalue and eigenvector of B", "1", t); } - //elpa_eigenvectors_all_host_arrays_dc(NEW_ELPA_HANDLE_POOL[handle_id], b, - // EigenValue, q, &info); - info=eigenvector(B, EigenValue, EigenVector); - if(loglevel>1) + // elpa_eigenvectors_all_host_arrays_dc(NEW_ELPA_HANDLE_POOL[handle_id], b, + // EigenValue, q, &info); + info = eigenvector(B, EigenValue, EigenVector); + if (loglevel > 1) { timer(myid, "calculate eigenvalue and eigenvector of B", "1", t); } MPI_Allreduce(&info, &allinfo, 1, MPI_INT, MPI_MAX, comm); allinfo = info; // calculate q*ev and put to work - for(int i=0; iDBL_MIN?1.0/sqrt(EigenValue[eidx]):0; - for(int j=0; j DBL_MIN ? 1.0 / sqrt(EigenValue[eidx]) : 0; + for (int j = 0; j < narows; ++j) + zwork[i * lda + j] = EigenVector[i * lda + j] * ev_sqrt; } // Form Q * diag(ev^{-1/2}) * Q^H in B. - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "qevq=qev*q^T", "2", t); } ScalapackConnector::gemm('N', 'C', nFull, nFull, nFull, 1.0, zwork.data(), EigenVector, 0.0, B, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "qevq=qev*q^T", "2", t); } } // if allinfo is still not 0 anyway, report error and quit - if(allinfo != 0) + if (allinfo != 0) { - if(myid == 0){ + if (myid == 0) + { ModuleBase::WARNING_QUIT("decomposeRightMatrix", - "Failed to decompose right matrix!\n info = " + std::to_string(allinfo)); + "Failed to decompose right matrix!\n info = " + std::to_string(allinfo)); } } return allinfo; @@ -319,28 +335,30 @@ int ELPA_Solver::decomposeRightMatrix(std::complex* B, double* EigenValu int ELPA_Solver::composeEigenVector(int DecomposedState, std::complex* B, std::complex* EigenVector) { double t = 0.0; - if(DecomposedState==1 || DecomposedState==2) + if (DecomposedState == 1 || DecomposedState == 2) { // transform the eigenvectors to original general equation, let U^-1*q, and put to q - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "Cpztrmm", "1", t); } ScalapackConnector::trmm('L', 'U', 'N', 'N', nFull, nev, 1.0, B, EigenVector, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "Cpztrmm", "1", t); } - } else { + } + else + { // transform the eigenvectors to original general equation, let b^C*q, and put to q - if(loglevel>1) + if (loglevel > 1) { - t=-1; + t = -1; timer(myid, "Cpzgemm", "1", t); } ScalapackConnector::gemm('C', 'N', nFull, nev, nFull, 1.0, B, zwork.data(), 0.0, EigenVector, desc); - if(loglevel>1) + if (loglevel > 1) { timer(myid, "Cpzgemm", "1", t); } @@ -356,53 +374,57 @@ int ELPA_Solver::composeEigenVector(int DecomposedState, std::complex* B // D: Diagonal matrix of eigenvalue // maxError: maximum absolute value of error // meanError: mean absolute value of error -void ELPA_Solver::verify(std::complex* A, double* EigenValue, std::complex* EigenVector, - double &maxError, double &meanError) +void ELPA_Solver::verify(std::complex* A, + double* EigenValue, + std::complex* EigenVector, + double& maxError, + double& meanError) { - std::complex* V=EigenVector; - const int naloc=narows*nacols; - std::complex* D=new std::complex[naloc]; - std::complex* R=zwork.data(); + std::complex* V = EigenVector; + const int naloc = narows * nacols; + std::complex* D = new std::complex[naloc]; + std::complex* R = zwork.data(); - for(int i=0; i2) saveMatrix("VD.dat", nFull, R, desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("VD.dat", nFull, R, desc, cblacs_ctxt); // R=A*V-V*D=A*V-R ScalapackConnector::hemm('L', 'U', nFull, 1.0, A, V, -1.0, R, desc); - if(loglevel>2) saveMatrix("AV-VD.dat", nFull, R, desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("AV-VD.dat", nFull, R, desc, cblacs_ctxt); // calculate the maximum and mean value of sum_i{R(:,i)*R(:,i)} - double sumError=0; - maxError=0; - for(int i=1; i<=nev; ++i) + double sumError = 0; + maxError = 0; + for (int i = 1; i <= nev; ++i) { std::complex E; - ScalapackConnector::dot(nFull, E, R, 1, i, 1, - R, 1, i, 1, desc); - double abs_E=std::abs(E); - sumError+=abs_E; - maxError=std::max(maxError, abs_E); + ScalapackConnector::dot(nFull, E, R, 1, i, 1, R, 1, i, 1, desc); + double abs_E = std::abs(E); + sumError += abs_E; + maxError = std::max(maxError, abs_E); } - meanError=sumError/nFull; + meanError = sumError / nFull; delete[] D; } @@ -414,57 +436,62 @@ void ELPA_Solver::verify(std::complex* A, double* EigenValue, std::compl // D: Diagonal matrix of eigenvalue // maxError: maximum absolute value of error // meanError: mean absolute value of error -void ELPA_Solver::verify(std::complex* A, std::complex* B, - double* EigenValue, std::complex* EigenVector, - double &maxError, double &meanError) +void ELPA_Solver::verify(std::complex* A, + std::complex* B, + double* EigenValue, + std::complex* EigenVector, + double& maxError, + double& meanError) { - std::complex* V=EigenVector; - const int naloc=narows*nacols; - std::complex* D=new std::complex[naloc]; - std::complex* R=new std::complex[naloc]; + std::complex* V = EigenVector; + const int naloc = narows * nacols; + std::complex* D = new std::complex[naloc]; + std::complex* R = new std::complex[naloc]; - for(int i=0; i2) saveMatrix("BV.dat", nFull, zwork.data(), desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("BV.dat", nFull, zwork.data(), desc, cblacs_ctxt); // R=B*V*D=zwork*D ScalapackConnector::hemm('L', 'U', nFull, 1.0, B, V, 0.0, zwork.data(), desc); - if(loglevel>2) saveMatrix("BVD.dat", nFull, R, desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("BVD.dat", nFull, R, desc, cblacs_ctxt); // R=A*V-B*V*D=A*V-R ScalapackConnector::hemm('L', 'U', nFull, 1.0, B, V, 0.0, zwork.data(), desc); - if(loglevel>2) saveMatrix("AV-BVD.dat", nFull, R, desc, cblacs_ctxt); + if (loglevel > 2) + saveMatrix("AV-BVD.dat", nFull, R, desc, cblacs_ctxt); // calculate the maximum and mean value of sum_i{R(:,i)*R(:,i)} - double sumError=0; - maxError=0; - for(int i=1; i<=nev; ++i) + double sumError = 0; + maxError = 0; + for (int i = 1; i <= nev; ++i) { std::complex E; - ScalapackConnector::dot(nFull, E, R, 1, i, 1, - R, 1, i, 1, desc); - double abs_E=std::abs(E); - sumError+=abs_E; - maxError=std::max(maxError, abs_E); + ScalapackConnector::dot(nFull, E, R, 1, i, 1, R, 1, i, 1, desc); + double abs_E = std::abs(E); + sumError += abs_E; + maxError = std::max(maxError, abs_E); } - meanError=sumError/nFull; + meanError = sumError / nFull; delete[] D; delete[] R; diff --git a/source/source_hsolver/module_genelpa/elpa_new_real.cpp b/source/source_hsolver/module_genelpa/elpa_new_real.cpp index cc6d0242bd0..830bdb49757 100644 --- a/source/source_hsolver/module_genelpa/elpa_new_real.cpp +++ b/source/source_hsolver/module_genelpa/elpa_new_real.cpp @@ -1,10 +1,8 @@ #include "elpa_new.h" #include "elpa_solver.h" - #include "source_base/module_external/scalapack_connector.h" -#include "utils.h" - #include "source_base/tool_quit.h" +#include "utils.h" #include #include @@ -56,10 +54,14 @@ int ELPA_Solver::generalized_eigenvector(double* A, { timer(myid, "decomposeRightMatrix", "1", t); } - if (allinfo != 0){ + if (allinfo != 0) + { // if allinfo is still not 0 anyway, report error and quit - if(myid == 0){ - ModuleBase::WARNING_QUIT("ELPA_Solver::generalized_eigenvector", "decomposeRightMatrix failed to decompose right matrix!\n info = " + std::to_string(allinfo)); + if (myid == 0) + { + ModuleBase::WARNING_QUIT("ELPA_Solver::generalized_eigenvector", + "decomposeRightMatrix failed to decompose right matrix!\n info = " + + std::to_string(allinfo)); } } @@ -321,11 +323,12 @@ int ELPA_Solver::decomposeRightMatrix(double* B, double* EigenValue, double* Eig } // if allinfo is still not 0 anyway, report error and quit - if(allinfo != 0) + if (allinfo != 0) { - if(myid == 0){ + if (myid == 0) + { ModuleBase::WARNING_QUIT("decomposeRightMatrix", - "Failed to decompose right matrix!\n info = " + std::to_string(allinfo)); + "Failed to decompose right matrix!\n info = " + std::to_string(allinfo)); } } return allinfo; @@ -343,7 +346,7 @@ int ELPA_Solver::composeEigenVector(int DecomposedState, double* B, double* Eige timer(myid, "Cpdtrmm", "1", t); } ScalapackConnector::trmm('L', 'U', 'N', 'N', nFull, nev, 1.0, B, EigenVector, desc); - //Cpdtrmm('L', 'U', 'N', 'N', nFull, nev, 1.0, B, EigenVector, desc); + // Cpdtrmm('L', 'U', 'N', 'N', nFull, nev, 1.0, B, EigenVector, desc); if (loglevel > 1) { timer(myid, "Cpdtrmm", "1", t); diff --git a/source/source_hsolver/module_genelpa/elpa_runtime_check.cpp b/source/source_hsolver/module_genelpa/elpa_runtime_check.cpp index c2b94211eea..904397d1c08 100644 --- a/source/source_hsolver/module_genelpa/elpa_runtime_check.cpp +++ b/source/source_hsolver/module_genelpa/elpa_runtime_check.cpp @@ -3,8 +3,6 @@ #include "elpa_new.h" #include "source_base/tool_quit.h" -#include - #include #include #include @@ -12,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -57,8 +56,7 @@ std::pair, std::vector> generate_hs() const double h_imag = distribution(generator); h[static_cast(i) * kN + j] = {h_real, h_imag}; - h[static_cast(j) * kN + i] - = i == j ? Complex{h_real, 0.0} : Complex{h_real, -h_imag}; + h[static_cast(j) * kN + i] = i == j ? Complex{h_real, 0.0} : Complex{h_real, -h_imag}; const double s_real = distribution(generator); const double s_imag = distribution(generator); @@ -203,7 +201,7 @@ bool generate_a_tilde(std::vector& result) double frobenius_norm(const std::vector& matrix) { double norm_squared = 0.0; - for (const Complex value : matrix) + for (const Complex value: matrix) { norm_squared += std::norm(value); } @@ -239,8 +237,7 @@ double relative_residual(const std::vector& matrix, const double eigenvector_norm = vector_norm(eigenvector, kN); const double denominator = (matrix_norm + std::abs(eigenvalue)) * eigenvector_norm; - return denominator > 0.0 ? vector_norm(residual.data(), kN) / denominator - : vector_norm(residual.data(), kN); + return denominator > 0.0 ? vector_norm(residual.data(), kN) / denominator : vector_norm(residual.data(), kN); } double orthogonality_error(const std::vector& eigenvectors) @@ -281,8 +278,7 @@ ElpaRuntimeCheckResult failed_result() bool needs_runtime_check(const int kernel) { - return kernel == ELPA_2STAGE_COMPLEX_AVX_BLOCK2 - || kernel == ELPA_2STAGE_COMPLEX_AVX2_BLOCK2 + return kernel == ELPA_2STAGE_COMPLEX_AVX_BLOCK2 || kernel == ELPA_2STAGE_COMPLEX_AVX2_BLOCK2 || kernel == ELPA_2STAGE_COMPLEX_AVX512_BLOCK2; } @@ -322,22 +318,18 @@ ElpaRuntimeCheckResult check_complex_block2_kernel(const int kernel) } #endif - if (!set_integer(handle, "na", kN, status) - || !set_integer(handle, "nev", kNev, status) - || !set_integer(handle, "local_nrows", kN, status) - || !set_integer(handle, "local_ncols", kN, status) + if (!set_integer(handle, "na", kN, status) || !set_integer(handle, "nev", kNev, status) + || !set_integer(handle, "local_nrows", kN, status) || !set_integer(handle, "local_ncols", kN, status) || !set_integer(handle, "nblk", kNblk, status) || !set_integer(handle, "mpi_comm_parent", MPI_Comm_c2f(MPI_COMM_SELF), status) - || !set_integer(handle, "process_row", 0, status) - || !set_integer(handle, "process_col", 0, status)) + || !set_integer(handle, "process_row", 0, status) || !set_integer(handle, "process_col", 0, status)) { deallocate(status); return failed_result(); } status = elpa_setup(handle); - if (status != ELPA_OK - || !set_integer(handle, "solver", ELPA_SOLVER_2STAGE, status) + if (status != ELPA_OK || !set_integer(handle, "solver", ELPA_SOLVER_2STAGE, status) || !set_integer(handle, "complex_kernel", kernel, status)) { deallocate(status); @@ -355,12 +347,9 @@ ElpaRuntimeCheckResult check_complex_block2_kernel(const int kernel) double maximum_residual = 0.0; for (int band = 0; band < kNev; ++band) { - maximum_residual - = std::max(maximum_residual, - relative_residual(original_matrix, - &eigenvectors[index(0, band, kN)], - eigenvalues[band], - matrix_norm)); + maximum_residual = std::max( + maximum_residual, + relative_residual(original_matrix, &eigenvectors[index(0, band, kN)], eigenvalues[band], matrix_norm)); } const double orthogonality = orthogonality_error(eigenvectors); diff --git a/source/source_hsolver/module_genelpa/elpa_runtime_check.h b/source/source_hsolver/module_genelpa/elpa_runtime_check.h index 9ad6434bfac..48052d6d35b 100644 --- a/source/source_hsolver/module_genelpa/elpa_runtime_check.h +++ b/source/source_hsolver/module_genelpa/elpa_runtime_check.h @@ -1,7 +1,6 @@ #ifndef ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_ELPA_RUNTIME_CHECK_H #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_ELPA_RUNTIME_CHECK_H - #include void validate_elpa_complex_kernel(int kernel, int set_status, MPI_Comm comm); diff --git a/source/source_hsolver/module_genelpa/pblas.h b/source/source_hsolver/module_genelpa/pblas.h index 5efca578455..cf4753331b7 100644 --- a/source/source_hsolver/module_genelpa/pblas.h +++ b/source/source_hsolver/module_genelpa/pblas.h @@ -1,44 +1,183 @@ #ifndef ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_PBLAS_H #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_PBLAS_H -void pddot_(int* n, double* dot, double* x, int* ix, int* jx, int* descx, int* incx, - double* y, int* iy, int* jy, int* descy, int* incy); - -void pzdotc_(int* n, double _Complex* dot, double _Complex* x, int* ix, int* jx, int* descx, int* incx, - double _Complex* y, int* iy, int* jy, int* descy, int* incy); -void pdsymv_(char* uplo, int* n, - double* alpha, double* a, int* ia, int* ja, int* desca, - double* x, int* ix, int* jx, int* descx, int* incx, - double* beta, double* y, int* iy, int* jy, int* descy, int* incy); -void pdtran_(int* m , int* n , - double* alpha , double* a , int* ia , int* ja , int* desca , - double* beta , double* c , int* ic , int* jc , int* descc ); +void pddot_(int* n, + double* dot, + double* x, + int* ix, + int* jx, + int* descx, + int* incx, + double* y, + int* iy, + int* jy, + int* descy, + int* incy); -void pdgemm_(char* transa , char* transb , int* m , int* n , int* k , - double* alpha , double* a , int* ia , int* ja , int* desca , - double* b , int* ib , int* jb , int* descb , - double* beta , double* c , int* ic , int* jc , int* descc ); -void pzgemm_(char* transa , char* transb , int* m , int* n , int* k , - double _Complex* alpha , double _Complex* a , int* ia , int* ja , int* desca , - double _Complex* b , int* ib , int* jb , int* descb , - double _Complex* beta , double _Complex* c , int* ic , int* jc , int* descc ); -void pdsymm_(char* side , char* uplo , int* m , int* n , - double* alpha , double* a , int* ia , int* ja , int* desca , - double* b , int* ib , int* jb , int* descb , - double* beta , double* c , int* ic , int* jc , int* descc ); -void pzsymm_(char* side , char* uplo , int* m , int* n , - double _Complex* alpha , double _Complex* a , int* ia , int* ja , int* desca , - double _Complex* b , int* ib , int* jb , int* descb , - double _Complex* beta , double _Complex* c , int* ic , int* jc , int* descc ); -void pzhemm_(char* side , char* uplo , int* m , int* n , - double _Complex* alpha , double _Complex* a , int* ia , int* ja , int* desca , - double _Complex* b , int* ib , int* jb , int* descb , - double _Complex* beta , double _Complex* c , int* ic , int* jc , int* descc ); -void pdtrmm_(char* side , char* uplo , char* transa , char* diag , int* m , int* n , - double* alpha , double* a , int* ia , int* ja , int* desca , - double* b , int* ib , int* jb , int* descb ); -void pztrmm_(char* side , char* uplo , char* transa , char* diag , int* m , int* n , - double _Complex* alpha , double _Complex* a , int* ia , int* ja , int* desca , - double _Complex* b , int* ib , int* jb , int* descb ); +void pzdotc_(int* n, + double _Complex* dot, + double _Complex* x, + int* ix, + int* jx, + int* descx, + int* incx, + double _Complex* y, + int* iy, + int* jy, + int* descy, + int* incy); +void pdsymv_(char* uplo, + int* n, + double* alpha, + double* a, + int* ia, + int* ja, + int* desca, + double* x, + int* ix, + int* jx, + int* descx, + int* incx, + double* beta, + double* y, + int* iy, + int* jy, + int* descy, + int* incy); +void pdtran_(int* m, + int* n, + double* alpha, + double* a, + int* ia, + int* ja, + int* desca, + double* beta, + double* c, + int* ic, + int* jc, + int* descc); + +void pdgemm_(char* transa, + char* transb, + int* m, + int* n, + int* k, + double* alpha, + double* a, + int* ia, + int* ja, + int* desca, + double* b, + int* ib, + int* jb, + int* descb, + double* beta, + double* c, + int* ic, + int* jc, + int* descc); +void pzgemm_(char* transa, + char* transb, + int* m, + int* n, + int* k, + double _Complex* alpha, + double _Complex* a, + int* ia, + int* ja, + int* desca, + double _Complex* b, + int* ib, + int* jb, + int* descb, + double _Complex* beta, + double _Complex* c, + int* ic, + int* jc, + int* descc); +void pdsymm_(char* side, + char* uplo, + int* m, + int* n, + double* alpha, + double* a, + int* ia, + int* ja, + int* desca, + double* b, + int* ib, + int* jb, + int* descb, + double* beta, + double* c, + int* ic, + int* jc, + int* descc); +void pzsymm_(char* side, + char* uplo, + int* m, + int* n, + double _Complex* alpha, + double _Complex* a, + int* ia, + int* ja, + int* desca, + double _Complex* b, + int* ib, + int* jb, + int* descb, + double _Complex* beta, + double _Complex* c, + int* ic, + int* jc, + int* descc); +void pzhemm_(char* side, + char* uplo, + int* m, + int* n, + double _Complex* alpha, + double _Complex* a, + int* ia, + int* ja, + int* desca, + double _Complex* b, + int* ib, + int* jb, + int* descb, + double _Complex* beta, + double _Complex* c, + int* ic, + int* jc, + int* descc); +void pdtrmm_(char* side, + char* uplo, + char* transa, + char* diag, + int* m, + int* n, + double* alpha, + double* a, + int* ia, + int* ja, + int* desca, + double* b, + int* ib, + int* jb, + int* descb); +void pztrmm_(char* side, + char* uplo, + char* transa, + char* diag, + int* m, + int* n, + double _Complex* alpha, + double _Complex* a, + int* ia, + int* ja, + int* desca, + double _Complex* b, + int* ib, + int* jb, + int* descb); #endif // ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_PBLAS_H diff --git a/source/source_hsolver/module_genelpa/scalapack.h b/source/source_hsolver/module_genelpa/scalapack.h index 03a2b1a321b..56a1d4e771f 100644 --- a/source/source_hsolver/module_genelpa/scalapack.h +++ b/source/source_hsolver/module_genelpa/scalapack.h @@ -2,15 +2,61 @@ #define ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_SCALAPACK_H // scalapack -int numroc_(const int *N, const int *NB, const int *IPROC, const int *ISRCPROC, const int *NPROCS); -void descinit_(int *DESC, const int *M, const int *N, const int *MB, const int *NB, const int *IRSRC, const int *ICSRC, const int *ICTXT, const int *LLD, int *INFO); -void pdpotrf_(const char *UPLO, const int *N, double *A, const int *IA, const int *JA, const int *DESCA, int *INFO); -void pzpotrf_(const char *UPLO, const int *N, double _Complex *A, const int *IA, const int *JA, const int *DESCA, int *INFO); -void pdsyev_(const char *JOBZ, const char *UPLO, int *N, double *A, int *IA, int *JA, int *DESCA, - double *W, double *Z, int *IZ, int *JZ, int *DESCZ, double *WORK, int *LWORK, int *INFO); -void pdgemr2d_(int *M, int *N, double *A, int *IA, int *JA, int *DESCA, - double *B, int *IB, int *JB, int *DESCB, int *ICTXT); -void pzgemr2d_(int *M, int *N, double _Complex *A, int *IA, int *JA, int *DESCA, - double _Complex *B, int *IB, int *JB, int *DESCB, int *ICTXT); +int numroc_(const int* N, const int* NB, const int* IPROC, const int* ISRCPROC, const int* NPROCS); +void descinit_(int* DESC, + const int* M, + const int* N, + const int* MB, + const int* NB, + const int* IRSRC, + const int* ICSRC, + const int* ICTXT, + const int* LLD, + int* INFO); +void pdpotrf_(const char* UPLO, const int* N, double* A, const int* IA, const int* JA, const int* DESCA, int* INFO); +void pzpotrf_(const char* UPLO, + const int* N, + double _Complex* A, + const int* IA, + const int* JA, + const int* DESCA, + int* INFO); +void pdsyev_(const char* JOBZ, + const char* UPLO, + int* N, + double* A, + int* IA, + int* JA, + int* DESCA, + double* W, + double* Z, + int* IZ, + int* JZ, + int* DESCZ, + double* WORK, + int* LWORK, + int* INFO); +void pdgemr2d_(int* M, + int* N, + double* A, + int* IA, + int* JA, + int* DESCA, + double* B, + int* IB, + int* JB, + int* DESCB, + int* ICTXT); +void pzgemr2d_(int* M, + int* N, + double _Complex* A, + int* IA, + int* JA, + int* DESCA, + double _Complex* B, + int* IB, + int* JB, + int* DESCB, + int* ICTXT); #endif // ABACUS_SOURCE_HSOLVER_MODULE_GENELPA_SCALAPACK_H diff --git a/source/source_hsolver/module_genelpa/utils.cpp b/source/source_hsolver/module_genelpa/utils.cpp index a334627c3a4..57cfc2a4dd4 100644 --- a/source/source_hsolver/module_genelpa/utils.cpp +++ b/source/source_hsolver/module_genelpa/utils.cpp @@ -1,8 +1,8 @@ #include "utils.h" #include "source_base/module_external/blacs_connector.h" -#include "source_base/tool_quit.h" #include "source_base/module_external/scalapack_connector.h" +#include "source_base/tool_quit.h" #include #include @@ -30,7 +30,7 @@ void loadMatrix(const char FileName[], int nFull, double* a, int* desca, int bla } } - double* b = nullptr; // buffer + double* b = nullptr; // buffer const int MAX_BUFFER_SIZE = 1e9; // max buffer size is 1GB int N = nFull; @@ -125,7 +125,7 @@ void saveMatrix(const char FileName[], int nFull, double* a, int* desca, int bla matrixFile.width(24); } - double* b = nullptr; // buffer + double* b = nullptr; // buffer const int MAX_BUFFER_SIZE = 1e9; // max buffer size is 1GB int N = nFull; @@ -184,7 +184,7 @@ void loadMatrix(const char FileName[], int nFull, std::complex* a, int* } } - std::complex* b; // buffer + std::complex* b; // buffer const int MAX_BUFFER_SIZE = 1e9; // max buffer size is 1GB int N = nFull; @@ -280,7 +280,7 @@ void saveMatrix(const char FileName[], int nFull, std::complex* a, int* matrixFile.width(24); } - std::complex* b; // buffer + std::complex* b; // buffer const int MAX_BUFFER_SIZE = 1e9; // max buffer size is 1GB int N = nFull; diff --git a/source/source_hsolver/module_pexsi/dist_bcd_matrix.cpp b/source/source_hsolver/module_pexsi/dist_bcd_matrix.cpp index 93a6dcb5b0d..3e6fb8c70e2 100644 --- a/source/source_hsolver/module_pexsi/dist_bcd_matrix.cpp +++ b/source/source_hsolver/module_pexsi/dist_bcd_matrix.cpp @@ -1,5 +1,6 @@ #ifdef __PEXSI #include "dist_bcd_matrix.h" + #include #include extern "C" diff --git a/source/source_hsolver/module_pexsi/dist_bcd_matrix.h b/source/source_hsolver/module_pexsi/dist_bcd_matrix.h index 6822efdccd4..3644d1b9eff 100644 --- a/source/source_hsolver/module_pexsi/dist_bcd_matrix.h +++ b/source/source_hsolver/module_pexsi/dist_bcd_matrix.h @@ -1,9 +1,9 @@ #ifndef DISTBCDMATRIX_H #define DISTBCDMATRIX_H -#include - #include "source_hsolver/module_pexsi/dist_mat_tf.h" + +#include // a Block Cyclic Data Distribution matrix // http://www.netlib.org/utk/papers/factor/node3.html // local matrix elements is stored in column major diff --git a/source/source_hsolver/module_pexsi/dist_ccs_matrix.h b/source/source_hsolver/module_pexsi/dist_ccs_matrix.h index b4a25fdd9ea..ecae2824292 100644 --- a/source/source_hsolver/module_pexsi/dist_ccs_matrix.h +++ b/source/source_hsolver/module_pexsi/dist_ccs_matrix.h @@ -66,7 +66,7 @@ class DistCCSMatrix // total number of processes and the processes with data in int nprocs = 0; - int nproc_data =0; + int nproc_data = 0; MPI_Group group_data; MPI_Comm comm_data; @@ -83,7 +83,7 @@ class DistCCSMatrix int numColLocal; // the first column index in current process - int firstCol=0; + int firstCol = 0; // Array stores the indices to the nonzero row indices in rowptrLocal and nzvalLocal int* colptrLocal = nullptr; diff --git a/source/source_hsolver/module_pexsi/dist_mat_tf.cpp b/source/source_hsolver/module_pexsi/dist_mat_tf.cpp index 8af5d43ef1e..97e0300c8b4 100644 --- a/source/source_hsolver/module_pexsi/dist_mat_tf.cpp +++ b/source/source_hsolver/module_pexsi/dist_mat_tf.cpp @@ -1,18 +1,17 @@ #ifdef __PEXSI #include "dist_mat_tf.h" -#include +#include "dist_bcd_matrix.h" +#include "dist_ccs_matrix.h" #include #include #include #include #include +#include #include -#include "dist_bcd_matrix.h" -#include "dist_ccs_matrix.h" - namespace pexsi { // find the minimum index, the return value will be a non-negtive value index value if it is found, otherwise will be a @@ -397,7 +396,7 @@ int DistMatrixTransformer::transformBCDtoCCS(DistBCDMatrix& SRC_Matrix, receiver_size_process, receiver_displacement_process, buffer2ccsIndex); -// Do transformation + // Do transformation std::vector sender_buffer(sender_size); std::vector receiver_buffer(receiver_size); // put H to sender buffer @@ -425,7 +424,7 @@ int DistMatrixTransformer::transformBCDtoCCS(DistBCDMatrix& SRC_Matrix, &receiver_displacement_process[0], MPI_DOUBLE, COMM_TRANS); -// collect H from receiver buffer + // collect H from receiver buffer delete[] H_ccs; H_ccs = new double[receiver_size]; buffer2CCSvalue(receiver_size, &buffer2ccsIndex[0], &receiver_buffer[0], H_ccs); @@ -455,7 +454,7 @@ int DistMatrixTransformer::transformBCDtoCCS(DistBCDMatrix& SRC_Matrix, &receiver_displacement_process[0], MPI_DOUBLE, COMM_TRANS); -// collect S from receiver buffer + // collect S from receiver buffer delete[] S_ccs; S_ccs = new double[receiver_size]; buffer2CCSvalue(receiver_size, &buffer2ccsIndex[0], &receiver_buffer[0], S_ccs); @@ -524,8 +523,10 @@ int DistMatrixTransformer::transformCCStoBCD(DistCCSMatrix& SRC_Matrix, int g_col = SRC_Matrix.globalCol(icol); int recv_pcol_bcd; int recv_col = DST_Matrix.localCol(g_col, recv_pcol_bcd); - - for (int rowidx = SRC_Matrix.get_colptr_local()[icol] - 1; rowidx < SRC_Matrix.get_colptr_local()[icol + 1] - 1; ++rowidx) + + for (int rowidx = SRC_Matrix.get_colptr_local()[icol] - 1; + rowidx < SRC_Matrix.get_colptr_local()[icol + 1] - 1; + ++rowidx) { int g_row = SRC_Matrix.get_rowind_local()[rowidx] - 1; int recv_prow_bcd; @@ -621,7 +622,9 @@ int DistMatrixTransformer::transformCCStoBCD(DistCCSMatrix& SRC_Matrix, int g_col = SRC_Matrix.globalCol(icol); int recv_pcol_bcd; int recv_col = DST_Matrix.localCol(g_col, recv_pcol_bcd); - for (int rowidx = SRC_Matrix.get_colptr_local()[icol] - 1; rowidx < SRC_Matrix.get_colptr_local()[icol + 1] - 1; ++rowidx) + for (int rowidx = SRC_Matrix.get_colptr_local()[icol] - 1; + rowidx < SRC_Matrix.get_colptr_local()[icol + 1] - 1; + ++rowidx) { int g_row = SRC_Matrix.get_rowind_local()[rowidx] - 1; int recv_prow_bcd; @@ -755,7 +758,6 @@ int DistMatrixTransformer::transformCCStoBCD(DistCCSMatrix& SRC_Matrix, delete[] dst_index; delete[] receiver_index; delete[] receiver_buffer; - } deleteGroupCommTrans(GROUP_TRANS, COMM_TRANS); return 0; diff --git a/source/source_hsolver/module_pexsi/dist_mat_tf.h b/source/source_hsolver/module_pexsi/dist_mat_tf.h index 7d705a68540..86e3ca53622 100644 --- a/source/source_hsolver/module_pexsi/dist_mat_tf.h +++ b/source/source_hsolver/module_pexsi/dist_mat_tf.h @@ -1,8 +1,8 @@ #ifndef DIST_MAT_TF_H #define DIST_MAT_TF_H -#include #include +#include #include // transform a sparse matrix from block cyclic distribution (BCD) to Compressed Column Storage (CCS) distribution // they should have same MPI communicator diff --git a/source/source_hsolver/module_pexsi/pexsi_solver.h b/source/source_hsolver/module_pexsi/pexsi_solver.h index 71637e64ef9..e6091289787 100644 --- a/source/source_hsolver/module_pexsi/pexsi_solver.h +++ b/source/source_hsolver/module_pexsi/pexsi_solver.h @@ -27,104 +27,105 @@ class PEXSI_Solver //========================================================== // PEXSI related variables //========================================================== - /** + /** * @brief Number of terms in the pole expansion. - */ + */ static int pexsi_npole; - /** + /** * @brief Whether inertia counting is used at the very beginning. - */ + */ static bool pexsi_inertia; - /** + /** * @brief Maximum number of PEXSI iterations after each inertia counting procedure. - */ + */ static int pexsi_nmax; - /** + /** * @brief Whether to construct PSelInv communication pattern. - */ + */ static bool pexsi_comm; - /** - * @brief Whether to use symmetric storage space used by the Selected Inversion algorithm for symmetric matrices. - */ + /** + * @brief Whether to use symmetric storage space used by the Selected Inversion algorithm for symmetric matrices. + */ static bool pexsi_storage; - /** - * @brief Ordering strategy for factorization and selected inversion. - */ + /** + * @brief Ordering strategy for factorization and selected inversion. + */ static int pexsi_ordering; - /** - * @brief row permutation strategy for factorization and selected inversion. - */ + /** + * @brief row permutation strategy for factorization and selected inversion. + */ static int pexsi_row_ordering; - /** + /** * @brief Number of processors for PARMETIS/PT-SCOTCH. Only used if the ordering == 0. - */ + */ static int pexsi_nproc; - /** + /** * @brief Matrix structure. * - = 0 : Unsymmetric matrix * - = 1 : Symmetric matrix (default). - */ + */ static bool pexsi_symm; - /** + /** * @brief Transpose. * - = 0 : Factor non transposed matrix (default). * - = 1 : Factor transposed matrix. - */ + */ static bool pexsi_trans; - /** + /** * @brief The pole expansion method to be used. * - = 1 : Cauchy Contour Integral method used. * - = 2 : Moussa optimized method. - */ + */ static int pexsi_method; - /** + /** * @brief The point parallelizaion of PEXSI. * - = 2 : Recommend two points parallelization - */ + */ static int pexsi_nproc_pole; - /** - * @brief Temperature, in the same unit as H - */ + /** + * @brief Temperature, in the same unit as H + */ static double pexsi_temp; - /** + /** * @brief Spectral gap. **Note** This can be set to be 0 in most cases. - */ + */ static double pexsi_gap; - /** + /** * @brief An upper bound for the spectral radius of \f$S^{-1} H\f$. - */ + */ static double pexsi_delta_e; - /** + /** * @brief Initial guess of lower bound for mu. - */ + */ static double pexsi_mu_lower; - /** + /** * @brief Initial guess of upper bound for mu. - */ + */ static double pexsi_mu_upper; - /** + /** * @brief Initial guess for mu (for the solver) (AG) - */ + */ static double pexsi_mu; - /** + /** * @brief Stopping criterion in terms of the chemical potential for the inertia counting procedure. - */ + */ static double pexsi_mu_thr; - /** + /** * @brief If the chemical potential is not in the initial interval, the interval is expanded by muInertiaExpansion. - */ + */ static double pexsi_mu_expand; - /** + /** * @brief Safe guard criterion in terms of the chemical potential to reinvoke the inertia counting procedure. - */ + */ static double pexsi_mu_guard; - /** - * @brief Stopping criterion of the %PEXSI iteration in terms of the number of electrons compared to numElectronExact. - */ + /** + * @brief Stopping criterion of the %PEXSI iteration in terms of the number of electrons compared to + * numElectronExact. + */ static double pexsi_elec_thr; - /** + /** * @brief If the absolute value of CCS matrix element is less than this value, it will be considered as zero. - */ + */ static double pexsi_zero_thr; private: diff --git a/source/source_hsolver/module_pexsi/simple_pexsi.cpp b/source/source_hsolver/module_pexsi/simple_pexsi.cpp index d91a6ad8650..a43103b98d7 100644 --- a/source/source_hsolver/module_pexsi/simple_pexsi.cpp +++ b/source/source_hsolver/module_pexsi/simple_pexsi.cpp @@ -3,15 +3,6 @@ // the Density Matrix and Energy Density Matrix calculated by PEXSI are transformed to 2D block cyclic distribution // #include "mpi.h" #ifdef __PEXSI -#include - -#include -#include -#include -#include -#include -#include - #include "c_pexsi_interface.h" #include "dist_bcd_matrix.h" #include "dist_ccs_matrix.h" @@ -21,6 +12,14 @@ #include "source_base/tool_quit.h" #include "source_hsolver/diago_pexsi.h" +#include +#include +#include +#include +#include +#include +#include + namespace pexsi { inline void strtolower(char* sa, char* sb) @@ -246,7 +245,7 @@ int simplePEXSI(MPI_Comm comm_PEXSI, plan = PPEXSIPlanInitialize(comm_PEXSI, pexsi_prow, pexsi_pcol, outputFileIndex, &info); } ModuleBase::timer::end("Diago_LCAO_Matrix", "PEXSIPlanInit"); - + ModuleBase::timer::end("Diago_LCAO_Matrix", "setup_PEXSI_plan"); // create compressed column storage distribution matrix parameter @@ -256,7 +255,6 @@ int simplePEXSI(MPI_Comm comm_PEXSI, // LiuXh modify 2021-03-30, add DONE(ofs_running,"xx") for test // DONE(ofs_running,"create compressed column storage distribution matrix parameter, finish"); - // create block cyclic distribution matrix parameter DistBCDMatrix SRC_Matrix(comm_2D, group_2D, blacs_ctxt, size, nblk, nrow, ncol, layout); // LiuXh modify 2021-03-30, add DONE(ofs_running,"xx") for test @@ -296,16 +294,16 @@ int simplePEXSI(MPI_Comm comm_PEXSI, int numTotalInertiaIter = 0; // Number of total inertia[out] // LiuXh modify 2021-04-29, add DONE(ofs_running,"xx") for test ModuleBase::timer::start("Diago_LCAO_Matrix", "PEXSIDFT"); - PPEXSIDFTDriver2(plan, // PEXSI plan[in] - &options, // PEXSI Options[in] - numElectronExact, // exact electron number[in] - &mu, // chemical potential[out] - &nelec, // number of electrons[out] - // &muMinInertia, // Lower bound for mu after the last inertia[out] - // &muMaxInertia, // Upper bound for mu after the last inertia[out] - &numTotalInertiaIter, // Number of total inertia[out] - // &numTotalPEXSIIter, // number of total pexsi evaluation procedure[out] - &info); // 0: successful; otherwise: unsuccessful + PPEXSIDFTDriver2(plan, // PEXSI plan[in] + &options, // PEXSI Options[in] + numElectronExact, // exact electron number[in] + &mu, // chemical potential[out] + &nelec, // number of electrons[out] + // &muMinInertia, // Lower bound for mu after the last inertia[out] + // &muMaxInertia, // Upper bound for mu after the last inertia[out] + &numTotalInertiaIter, // Number of total inertia[out] + // &numTotalPEXSIIter, // number of total pexsi evaluation procedure[out] + &info); // 0: successful; otherwise: unsuccessful // LiuXh modify 2021-04-29, add DONE(ofs_running,"xx") for test ModuleBase::timer::end("Diago_LCAO_Matrix", "PEXSIDFT"); diff --git a/source/source_hsolver/test/diago_bpcg_test.cpp b/source/source_hsolver/test/diago_bpcg_test.cpp index 34bed8b4735..0cd150985ca 100644 --- a/source/source_hsolver/test/diago_bpcg_test.cpp +++ b/source/source_hsolver/test/diago_bpcg_test.cpp @@ -1,16 +1,17 @@ -#include "source_base/inverse_matrix.h" -#include "source_base/module_external/lapack_connector.h" -#include "source_base/parallel_comm.h" -#include "source_psi/psi.h" -#include "../diago_iter_assist.h" #include "../diago_bpcg.h" + +#include "../diago_iter_assist.h" #include "diago_mock.h" #include "mpi.h" +#include "source_base/inverse_matrix.h" +#include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" +#include "source_psi/psi.h" -#include #include #include +#include #include /************************************************ @@ -34,20 +35,21 @@ */ // call lapack in order to compare to bpcg -void lapackEigen(int &npw, std::vector> &hm, double *e, bool outtime = false) +void lapackEigen(int& npw, std::vector>& hm, double* e, bool outtime = false) { clock_t start, end; start = clock(); int lwork = 2 * npw; - std::complex *work2 = new std::complex[lwork]; - double *rwork = new double[3 * npw - 2]; + std::complex* work2 = new std::complex[lwork]; + double* rwork = new double[3 * npw - 2]; int info = 0; char tmp_c1 = 'V', tmp_c2 = 'U'; zheev_(&tmp_c1, &tmp_c2, &npw, hm.data(), &npw, e, work2, &lwork, rwork, &info); end = clock(); - if (outtime) { + if (outtime) + { std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; -} + } delete[] rwork; delete[] work2; } @@ -56,13 +58,12 @@ class DiagoBPCGPrepare { public: DiagoBPCGPrepare(int nband, int npw, int sparsity, bool reorder, double eps, int maxiter, double threshold) - : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), - threshold(threshold) + : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), threshold(threshold) { -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +#endif } int nband, npw, sparsity, maxiter, notconv; @@ -70,16 +71,18 @@ class DiagoBPCGPrepare double eps, avg_iter; bool reorder; double threshold; - int nprocs=1, mypnum=0; + int nprocs = 1, mypnum = 0; // threshold is the comparison standard between bpcg and lapack - void CompareEigen(double *precondition) + void CompareEigen(double* precondition) { // calculate eigenvalues by LAPACK; - double *e_lapack = new double[npw]; + double* e_lapack = new double[npw]; auto ev = DIAGOTEST::hmatrix; - if(mypnum == 0) { 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); @@ -88,43 +91,47 @@ class DiagoBPCGPrepare { for (int j = 0; j < npw; j++) { - double rand = static_cast(u(p))/10.; + double rand = static_cast(u(p)) / 10.; // psiguess(i,j) = ev(j,i)*(1+rand); psiguess(i, j) = ev[j * DIAGOTEST::h_nc + i] * rand; } } // run bpcg - //====================================================================== - double *en = new double[npw]; + //====================================================================== + double* en = new double[npw]; int ik = 1; - HSOperatorMock> ha; - int* ngk = new int [1]; - //psi::Psi> psi(ngk,ik,nband,npw); - psi::Psi> psi; - psi.resize(ik,nband,npw); - //psi.fix_k(0); + HSOperatorMock> ha; + int* ngk = new int[1]; + // psi::Psi> psi(ngk,ik,nband,npw); + psi::Psi> 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); - } - } + psi(i, j) = psiguess(i, j); + } + } psi::Psi> psi_local; double* precondition_local; DIAGOTEST::npw_local = new int[nprocs]; -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); - DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix, DIAGOTEST::hmatrix_local); //will distribute psi and Hmatrix to each process - precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; - DIAGOTEST::divide_psi(precondition,precondition_local); +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::divide_hpsi(psi, + psi_local, + DIAGOTEST::hmatrix, + DIAGOTEST::hmatrix_local); // will distribute psi and Hmatrix to each process + precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; + DIAGOTEST::divide_psi(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> bpcg(precondition_local); psi_local.fix_k(0); @@ -138,10 +145,10 @@ class DiagoBPCGPrepare bpcg.diag(ha, psi_local.get_pointer(), en, ethr_band); bpcg.diag(ha, psi_local.get_pointer(), en, ethr_band); end = MPI_Wtime(); - //if(mypnum == 0) printf("diago time:%7.3f\n",end-start); - delete [] DIAGOTEST::npw_local; - delete [] precondition_local; - //====================================================================== + // 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); @@ -159,12 +166,12 @@ class DiagoBPCGTest : public ::testing::TestWithParam TEST_P(DiagoBPCGTest, RandomHamilt) { DiagoBPCGPrepare dcp = GetParam(); - //std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" + // std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" // << dcp.sparsity << ", eps=" << dcp.eps << std::endl; hsolver::DiagoIterAssist>::PW_DIAG_NMAX = dcp.maxiter; hsolver::DiagoIterAssist>::PW_DIAG_THR = dcp.eps; - //std::cout<<"maxiter "<>::PW_DIAG_NMAX<>::PW_DIAG_THR<>::PW_DIAG_NMAX<>::PW_DIAG_THR<> hpsi(dcp.nband, dcp.npw, dcp.sparsity); DIAGOTEST::hmatrix = hpsi.hamilt(); @@ -178,12 +185,12 @@ INSTANTIATE_TEST_SUITE_P(VerifyCG, ::testing::Values( // nband, npw, sparsity, reorder, eps, maxiter, threshold DiagoBPCGPrepare(10, 500, 0, true, 1e-5, 300, 5e-2) - // DiagoBPCGPrepare(20, 500, 6, true, 1e-5, 300, 5e-2) - // DiagoBPCGPrepare(20, 1000, 8, true, 1e-5, 300, 5e-2), - // DiagoBPCGPrepare(40, 1000, 8, true, 1e-6, 300, 5e-2) - )); - //DiagoBPCGPrepare(40, 2000, 8, true, 1e-5, 500, 1e-2))); - // the last one is passed but time-consumming. + // DiagoBPCGPrepare(20, 500, 6, true, 1e-5, 300, 5e-2) + // DiagoBPCGPrepare(20, 1000, 8, true, 1e-5, 300, 5e-2), + // DiagoBPCGPrepare(40, 1000, 8, true, 1e-6, 300, 5e-2) + )); +// DiagoBPCGPrepare(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 // in generating a Hermite matrix @@ -222,7 +229,6 @@ TEST(DiagoBPCGTest, Hamilt) } }*/ - TEST(DiagoBPCGTest, readH) { // read Hamilt matrix from file data-H @@ -252,31 +258,33 @@ TEST(DiagoBPCGTest, readH) dcp.CompareEigen(hpsi.precond()); } -int main(int argc, char **argv) +int main(int argc, char** argv) { - int nproc = 1, myrank = 0; + int nproc = 1, myrank = 0; #ifdef __MPI - int nproc_in_pool, kpar=1, mypool, rank_in_pool; - setupmpi(argc,argv,nproc, myrank); + 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); MPI_Comm_split(MPI_COMM_WORLD, myrank, 0, &BP_WORLD); #else - MPI_Init(&argc, &argv); + 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()); -} + ::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 0; + return 0; } diff --git a/source/source_hsolver/test/diago_cg_float_test.cpp b/source/source_hsolver/test/diago_cg_float_test.cpp index 0651daca593..eb96ab39c03 100644 --- a/source/source_hsolver/test/diago_cg_float_test.cpp +++ b/source/source_hsolver/test/diago_cg_float_test.cpp @@ -1,11 +1,11 @@ -#include "../diago_cg.h" #include "../diag_comm_info.h" +#include "../diago_cg.h" #include "../diago_iter_assist.h" #include "diago_mock.h" #include "mpi.h" -#include "source_base/parallel_comm.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" #include "source_psi/psi.h" @@ -35,20 +35,21 @@ */ // call lapack in order to compare to cg -void lapackEigen(int &npw, std::vector> &hm, float *e, bool outtime = false) +void lapackEigen(int& npw, std::vector>& hm, float* e, bool outtime = false) { clock_t start, end; start = clock(); int lwork = 2 * npw; - std::complex *work2 = new std::complex[lwork]; - float *rwork = new float[3 * npw - 2]; + std::complex* work2 = new std::complex[lwork]; + float* rwork = new float[3 * npw - 2]; int info = 0; char tmp_c1 = 'V', tmp_c2 = 'U'; cheev_(&tmp_c1, &tmp_c2, &npw, hm.data(), &npw, e, work2, &lwork, rwork, &info); end = clock(); - if (outtime) { + if (outtime) + { std::cout << "Lapack Run time: " << (float)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; -} + } delete[] rwork; delete[] work2; } @@ -57,13 +58,12 @@ class DiagoCGPrepare { public: DiagoCGPrepare(int nband, int npw, int sparsity, bool reorder, float eps, int maxiter, float threshold) - : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), - threshold(threshold) + : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), threshold(threshold) { -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +#endif } int nband, npw, sparsity, maxiter, notconv; @@ -71,17 +71,19 @@ class DiagoCGPrepare float eps, avg_iter; bool reorder; float threshold; - int nprocs=1, mypnum=0; + int nprocs = 1, mypnum = 0; // threshold is the comparison standard between cg and lapack - void CompareEigen(float *precondition) + void CompareEigen(float* precondition) { // calculate eigenvalues by LAPACK; - float *e_lapack = new float[npw]; + float* e_lapack = new float[npw]; auto ev = DIAGOTEST::hmatrix_f; - if(mypnum == 0) { lapackEigen(npw, ev, e_lapack, false); -} + if (mypnum == 0) + { + lapackEigen(npw, ev, e_lapack, false); + } // initial guess of psi by perturbing lapack psi std::vector> psiguess(nband * npw); @@ -92,46 +94,50 @@ class DiagoCGPrepare { for (int j = 0; j < npw; j++) { - float rand = static_cast(u(p))/10.; + float rand = static_cast(u(p)) / 10.; // psiguess(i,j) = ev(j,i)*(1+rand); psiguess[i * npw + j] = ev[j * DIAGOTEST::h_nc + i] * rand; } } // run cg - //====================================================================== - float *en = new float[npw]; + //====================================================================== + float* en = new float[npw]; int ik = 1; - HSOperatorMock> ha; - psi::Psi> psi; - psi.resize(ik,nband,npw); - //psi.fix_k(0); + HSOperatorMock> ha; + psi::Psi> 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 * npw + j]; - } - } + psi(i, j) = psiguess[i * npw + j]; + } + } psi::Psi> psi_local; float* precondition_local; DIAGOTEST::npw_local = new int[nprocs]; -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); - DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_f, DIAGOTEST::hmatrix_local_f); //will distribute psi and Hmatrix to each process - precondition_local = new float[DIAGOTEST::npw_local[mypnum]]; +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::divide_hpsi(psi, + psi_local, + DIAGOTEST::hmatrix_f, + DIAGOTEST::hmatrix_local_f); // will distribute psi and Hmatrix to each process + precondition_local = new float[DIAGOTEST::npw_local[mypnum]]; DIAGOTEST::divide_psi(precondition, precondition_local); #else - DIAGOTEST::hmatrix_local_f = DIAGOTEST::hmatrix_f; - DIAGOTEST::npw_local[0] = DIAGOTEST::npw; - psi_local = psi; - precondition_local = new float[DIAGOTEST::npw]; - for(int i=0;i> cg(precondition_local); psi_local.fix_k(0); - // cg.diag(ha,psi_local,en); + // cg.diag(ha,psi_local,en); /**************************************************************/ // New interface of cg method /**************************************************************/ @@ -151,9 +157,9 @@ class DiagoCGPrepare float start, end; start = MPI_Wtime(); - std::vector ethr_band(nband, 1e-5); - cg.diag(ha, psi_local.get_nbasis(), + cg.diag(ha, + psi_local.get_nbasis(), psi_local.get_nbands(), psi_local.get_current_ngk(), psi_local.get_pointer(), @@ -163,13 +169,13 @@ class DiagoCGPrepare /**************************************************************/ end = MPI_Wtime(); - //if(mypnum == 0) printf("diago time:%7.3f\n",end-start); - delete [] DIAGOTEST::npw_local; - delete [] precondition_local; - //====================================================================== + // 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)< TEST_P(DiagoCGFloatTest, RandomHamilt) { DiagoCGPrepare dcp = GetParam(); - //std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" + // std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" // << dcp.sparsity << ", eps=" << dcp.eps << std::endl; hsolver::DiagoIterAssist>::PW_DIAG_NMAX = dcp.maxiter; hsolver::DiagoIterAssist>::PW_DIAG_THR = dcp.eps; - //std::cout<<"maxiter "<>::PW_DIAG_NMAX<>::PW_DIAG_THR<>::PW_DIAG_NMAX<>::PW_DIAG_THR<> hpsi(dcp.nband, dcp.npw, dcp.sparsity); DIAGOTEST::hmatrix_f = hpsi.hamilt(); @@ -203,12 +209,12 @@ INSTANTIATE_TEST_SUITE_P(VerifyCG, ::testing::Values( // nband, npw, sparsity, reorder, eps, maxiter, threshold DiagoCGPrepare(10, 200, 0, true, 1e-6, 300, 1e-0) - // DiagoCGPrepare(10, 200, 6, true, 1e-6, 300, 1e-0) - // DiagoCGPrepare(10, 400, 8, true, 1e-6, 300, 1e-0), - // DiagoCGPrepare(10, 600, 8, true, 1e-6, 300, 1e-0) - )); - //DiagoCGPrepare(40, 2000, 8, true, 1e-5, 500, 1e-2))); - // the last one is passed but time-consumming. + // DiagoCGPrepare(10, 200, 6, true, 1e-6, 300, 1e-0) + // DiagoCGPrepare(10, 400, 8, true, 1e-6, 300, 1e-0), + // DiagoCGPrepare(10, 600, 8, true, 1e-6, 300, 1e-0) + )); +// 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 // in generating a Hermite matrix @@ -299,30 +305,32 @@ TEST(DiagoCGFloatTest, readH) dcp.CompareEigen(hpsi.precond()); } -int main(int argc, char **argv) +int main(int argc, char** argv) { - int nproc = 1, myrank = 0; + int nproc = 1, myrank = 0; #ifdef __MPI - int nproc_in_pool, kpar=1, mypool, rank_in_pool; - setupmpi(argc,argv,nproc, myrank); + 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); #else - MPI_Init(&argc, &argv); + 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()); -} + ::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 0; + return 0; } diff --git a/source/source_hsolver/test/diago_cg_real_test.cpp b/source/source_hsolver/test/diago_cg_real_test.cpp index f1afbd99a45..3e66032876b 100644 --- a/source/source_hsolver/test/diago_cg_real_test.cpp +++ b/source/source_hsolver/test/diago_cg_real_test.cpp @@ -1,11 +1,11 @@ -#include "../diago_cg.h" #include "../diag_comm_info.h" +#include "../diago_cg.h" #include "../diago_iter_assist.h" #include "diago_mock.h" #include "mpi.h" -#include "source_base/parallel_comm.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" #include "source_psi/psi.h" @@ -18,23 +18,23 @@ * unit test of functions in Diago_CG ***********************************************/ - /** - * Class Diago_CG is an approach for eigenvalue problems - * This unittest test the function Diago_CG::diag() for FPTYPE=double, Device=cpu - * with different examples. - * - the Hermite matrices (npw=500,1000) produced using random numbers and with sparsity of 0%, 60%, 80% - * - the Hamiltonian matrix read from "data-H", produced by using out_hs in INPUT of a LCAO calculation - * - a 2x2 Hermite matrix for learning and checking - * - * Note: - * The test is passed when the eignvalues are closed to these calculated by LAPACK. - * It is used together with a header file diago_mock.h. - * The default Hermite matrix generated here is real symmetric, one can add an imaginary part - * by changing two commented out lines in diago_mock.h. - * - */ +/** + * Class Diago_CG is an approach for eigenvalue problems + * This unittest test the function Diago_CG::diag() for FPTYPE=double, Device=cpu + * with different examples. + * - the Hermite matrices (npw=500,1000) produced using random numbers and with sparsity of 0%, 60%, 80% + * - the Hamiltonian matrix read from "data-H", produced by using out_hs in INPUT of a LCAO calculation + * - a 2x2 Hermite matrix for learning and checking + * + * Note: + * The test is passed when the eignvalues are closed to these calculated by LAPACK. + * It is used together with a header file diago_mock.h. + * The default Hermite matrix generated here is real symmetric, one can add an imaginary part + * by changing two commented out lines in diago_mock.h. + * + */ - // call lapack in order to compare to cg +// call lapack in order to compare to cg void lapackEigen(int& npw, std::vector& hm, double* e, bool outtime = false) { int info = 0; @@ -45,30 +45,33 @@ void lapackEigen(int& npw, std::vector& hm, double* e, bool outtime = fa double work_tmp; constexpr int minus_one = -1; - dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, &work_tmp, &minus_one, &info); // get best lwork + dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, &work_tmp, &minus_one, &info); // get best lwork const int lwork = work_tmp; double* work2 = new double[lwork]; dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, &info); end = clock(); - if (info) { std::cout << "ERROR: Lapack solver, info=" << info << std::endl; -} - if (outtime) { std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; -} + if (info) + { + std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + } + if (outtime) + { + std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + } delete[] work2; } class DiagoCGPrepare { -public: + public: DiagoCGPrepare(int nband, int npw, int sparsity, bool reorder, double eps, int maxiter, double threshold) - : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), - threshold(threshold) + : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), threshold(threshold) { -#ifdef __MPI +#ifdef __MPI MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +#endif } int nband, npw, sparsity, maxiter, notconv; @@ -84,8 +87,10 @@ class DiagoCGPrepare // calculate eigenvalues by LAPACK; double* e_lapack = new double[npw]; auto ev = DIAGOTEST::hmatrix_d; - if (mypnum == 0) { lapackEigen(npw, ev, e_lapack, false); -} + if (mypnum == 0) + { + lapackEigen(npw, ev, e_lapack, false); + } // initial guess of psi by perturbing lapack psi ModuleBase::matrix psiguess(nband, npw); std::default_random_engine p(1); @@ -100,14 +105,14 @@ class DiagoCGPrepare } } // run cg - //====================================================================== + //====================================================================== double* en = new double[npw]; int ik = 1; HSOperatorMock ha; int* ngk = new int[1]; psi::Psi psi; psi.resize(ik, nband, npw); - //psi.fix_k(0); + // psi.fix_k(0); for (int i = 0; i < nband; i++) { for (int j = 0; j < npw; j++) @@ -119,9 +124,12 @@ class DiagoCGPrepare psi::Psi psi_local; double* precondition_local; DIAGOTEST::npw_local = new int[nprocs]; -#ifdef __MPI +#ifdef __MPI DIAGOTEST::cal_division(DIAGOTEST::npw); - DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_d, DIAGOTEST::hmatrix_local_d); //will distribute psi and Hmatrix to each process + DIAGOTEST::divide_hpsi(psi, + psi_local, + DIAGOTEST::hmatrix_d, + DIAGOTEST::hmatrix_local_d); // will distribute psi and Hmatrix to each process precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; DIAGOTEST::divide_psi(precondition, precondition_local); #else @@ -129,7 +137,8 @@ class DiagoCGPrepare 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]; + for (int i = 0; i < DIAGOTEST::npw; i++) + precondition_local[i] = precondition[i]; #endif // hsolver::DiagoCG cg(precondition_local); psi_local.fix_k(0); @@ -137,7 +146,6 @@ class DiagoCGPrepare // start = MPI_Wtime(); // cg.diag(ha, psi_local, en); - /**************************************************************/ // New interface of cg method /**************************************************************/ @@ -157,9 +165,9 @@ class DiagoCGPrepare double start, end; start = MPI_Wtime(); - std::vector ethr_band(nband, 1e-5); - cg.diag(ha, psi_local.get_nbasis(), + cg.diag(ha, + psi_local.get_nbasis(), psi_local.get_nbands(), psi_local.get_current_ngk(), psi_local.get_pointer(), @@ -169,7 +177,7 @@ class DiagoCGPrepare /**************************************************************/ end = MPI_Wtime(); - //if(mypnum == 0) printf("diago time:%7.3f\n",end-start); + // if(mypnum == 0) printf("diago time:%7.3f\n",end-start); delete[] DIAGOTEST::npw_local; delete[] precondition_local; //====================================================================== @@ -201,16 +209,16 @@ TEST_P(DiagoCGTest, RandomHamilt) } INSTANTIATE_TEST_SUITE_P(VerifyCG, - DiagoCGTest, - ::testing::Values( - // nband, npw, sparsity, reorder, eps, maxiter, threshold - 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))); -// the last one is passed but time-consumming. + DiagoCGTest, + ::testing::Values( + // nband, npw, sparsity, reorder, eps, maxiter, threshold + 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))); +// the last one is passed but time-consumming. // check that the mock class HPsi work well // in generating a Hermite matrix @@ -291,8 +299,10 @@ int main(int argc, char** argv) testing::InitGoogleTest(&argc, argv); ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners(); - if (myrank != 0) { delete listeners.Release(listeners.default_result_printer()); -} + if (myrank != 0) + { + delete listeners.Release(listeners.default_result_printer()); + } int result = RUN_ALL_TESTS(); if (myrank == 0 && result != 0) diff --git a/source/source_hsolver/test/diago_cg_test.cpp b/source/source_hsolver/test/diago_cg_test.cpp index 853b51c9f71..e3dd5adafc4 100644 --- a/source/source_hsolver/test/diago_cg_test.cpp +++ b/source/source_hsolver/test/diago_cg_test.cpp @@ -4,9 +4,9 @@ #include "../diago_iter_assist.h" #include "diago_mock.h" #include "mpi.h" -#include "source_base/parallel_comm.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" #include "source_psi/psi.h" @@ -36,20 +36,21 @@ */ // call lapack in order to compare to cg -void lapackEigen(int &npw, std::vector> &hm, double *e, bool outtime = false) +void lapackEigen(int& npw, std::vector>& hm, double* e, bool outtime = false) { clock_t start, end; start = clock(); int lwork = 2 * npw; - std::complex *work2 = new std::complex[lwork]; - double *rwork = new double[3 * npw - 2]; + std::complex* work2 = new std::complex[lwork]; + double* rwork = new double[3 * npw - 2]; int info = 0; char tmp_c1 = 'V', tmp_c2 = 'U'; zheev_(&tmp_c1, &tmp_c2, &npw, hm.data(), &npw, e, work2, &lwork, rwork, &info); end = clock(); - if (outtime) { + if (outtime) + { std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; -} + } delete[] rwork; delete[] work2; } @@ -58,13 +59,12 @@ class DiagoCGPrepare { public: DiagoCGPrepare(int nband, int npw, int sparsity, bool reorder, double eps, int maxiter, double threshold) - : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), - threshold(threshold) + : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), threshold(threshold) { -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +#endif } int nband, npw, sparsity, maxiter, notconv; @@ -72,16 +72,18 @@ class DiagoCGPrepare double eps, avg_iter; bool reorder; double threshold; - int nprocs=1, mypnum=0; + int nprocs = 1, mypnum = 0; // threshold is the comparison standard between cg and lapack - void CompareEigen(double *precondition) + void CompareEigen(double* precondition) { // calculate eigenvalues by LAPACK; - double *e_lapack = new double[npw]; + double* e_lapack = new double[npw]; auto ev = DIAGOTEST::hmatrix; - if(mypnum == 0) { 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); @@ -90,41 +92,45 @@ class DiagoCGPrepare { for (int j = 0; j < npw; j++) { - double rand = static_cast(u(p))/10.; + double rand = static_cast(u(p)) / 10.; // psiguess(i,j) = ev(j,i)*(1+rand); psiguess(i, j) = ev[j * DIAGOTEST::h_nc + i] * rand; } } // run cg - //====================================================================== - double *en = new double[npw]; + //====================================================================== + double* en = new double[npw]; int ik = 1; - HSOperatorMock> ha; - psi::Psi> psi; - psi.resize(ik,nband,npw); - //psi.fix_k(0); + HSOperatorMock> ha; + psi::Psi> 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); - } - } + psi(i, j) = psiguess(i, j); + } + } psi::Psi> psi_local; double* precondition_local; DIAGOTEST::npw_local = new int[nprocs]; -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); - DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix, DIAGOTEST::hmatrix_local); //will distribute psi and Hmatrix to each process - precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; - DIAGOTEST::divide_psi(precondition,precondition_local); +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::divide_hpsi(psi, + psi_local, + DIAGOTEST::hmatrix, + DIAGOTEST::hmatrix_local); // will distribute psi and Hmatrix to each process + precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; + DIAGOTEST::divide_psi(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 ethr_band(nband, 1e-5); - cg.diag(ha, psi_local.get_nbasis(), + cg.diag(ha, + psi_local.get_nbasis(), psi_local.get_nbands(), psi_local.get_current_ngk(), psi_local.get_pointer(), @@ -157,12 +163,12 @@ class DiagoCGPrepare precondition_local); /**************************************************************/ - // cg.diag(ha,psi_local,en); + // 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; - //====================================================================== + // 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); @@ -180,12 +186,12 @@ class DiagoCGTest : public ::testing::TestWithParam TEST_P(DiagoCGTest, RandomHamilt) { DiagoCGPrepare dcp = GetParam(); - //std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" + // std::cout << "npw=" << dcp.npw << ", nband=" << dcp.nband << ", sparsity=" // << dcp.sparsity << ", eps=" << dcp.eps << std::endl; hsolver::DiagoIterAssist>::PW_DIAG_NMAX = dcp.maxiter; hsolver::DiagoIterAssist>::PW_DIAG_THR = dcp.eps; - //std::cout<<"maxiter "<>::PW_DIAG_NMAX<>::PW_DIAG_THR<>::PW_DIAG_NMAX<>::PW_DIAG_THR<> hpsi(dcp.nband, dcp.npw, dcp.sparsity); DIAGOTEST::hmatrix = hpsi.hamilt(); @@ -199,12 +205,12 @@ INSTANTIATE_TEST_SUITE_P(VerifyCG, ::testing::Values( // nband, npw, sparsity, reorder, eps, maxiter, threshold 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))); - // the last one is passed but time-consumming. + // 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))); +// the last one is passed but time-consumming. // check that the mock class HPsi work well // in generating a Hermite matrix @@ -295,30 +301,32 @@ TEST(DiagoCGTest, readH) dcp.CompareEigen(hpsi.precond()); } -int main(int argc, char **argv) +int main(int argc, char** argv) { - int nproc = 1, myrank = 0; + int nproc = 1, myrank = 0; #ifdef __MPI - int nproc_in_pool, kpar=1, mypool, rank_in_pool; - setupmpi(argc,argv,nproc, myrank); + 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); #else - MPI_Init(&argc, &argv); + 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()); -} + ::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 0; + return 0; } diff --git a/source/source_hsolver/test/diago_david_float_test.cpp b/source/source_hsolver/test/diago_david_float_test.cpp index 930b4920cb6..8120db8fb28 100644 --- a/source/source_hsolver/test/diago_david_float_test.cpp +++ b/source/source_hsolver/test/diago_david_float_test.cpp @@ -1,22 +1,22 @@ -#include"source_hsolver/diago_david.h" -#include "source_hsolver/diag_comm_info.h" -#include"source_hsolver/diago_iter_assist.h" -#include "source_base/parallel_comm.h" -#include"diago_mock.h" -#include "source_psi/psi.h" -#include"gtest/gtest.h" +#include "diago_mock.h" +#include "mpi.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" -#include"mpi.h" +#include "source_hsolver/diag_comm_info.h" +#include "source_hsolver/diago_david.h" +#include "source_hsolver/diago_iter_assist.h" +#include "source_psi/psi.h" + +#include "gtest/gtest.h" #define CONVTHRESHOLD 1e-0 #define DETAILINFO false - /************************************************ -* unit test of class Diago_David -***********************************************/ + * unit test of class Diago_David + ***********************************************/ /** * Class Diago_David is used to solve the eigenvalues @@ -26,226 +26,242 @@ * - the hamilt matrix (npw=100,500,1000) produced by random with sparsity of 50% * - the hamilt matrix (npw=100,500,1000) produced by random with sparsity of 0% * - the hamilt matrix read from "data-H" - * + * * The test is passed when the eignvalues are closed to these calculated by LAPACK. - * + * */ -//use lapack to calcualte eigenvalue of matrix hm -//NOTE: after finish this function, hm stores the eigen vectors. -void lapackEigen(int &npw, std::vector> &hm, float * e, bool outtime=false) +// use lapack to calcualte eigenvalue of matrix hm +// NOTE: after finish this function, hm stores the eigen vectors. +void lapackEigen(int& npw, std::vector>& hm, float* e, bool outtime = false) { - int lwork = 2 * npw; - std::complex *work2= new std::complex[lwork]; - float* rwork = new float[3*npw-2]; - int info = 0; - - auto tmp = hm; - - clock_t start,end; - start = clock(); - char tmp_c1 = 'V', tmp_c2 = 'U'; - cheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); - end = clock(); - if(info) { std::cout << "ERROR: Lapack solver, info=" << info <* work2 = new std::complex[lwork]; + float* rwork = new float[3 * npw - 2]; + int info = 0; + + auto tmp = hm; + + clock_t start, end; + start = clock(); + char tmp_c1 = 'V', tmp_c2 = 'U'; + cheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); + end = clock(); + if (info) + { + std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + } + if (outtime) + { + std::cout << "Lapack Run time: " << (float)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + } - delete [] rwork; - delete [] work2; + delete[] rwork; + delete[] work2; } -class DiagoDavPrepare +class DiagoDavPrepare { -public: - DiagoDavPrepare(int nband, int npw, int sparsity, int order,float eps,int maxiter): - nband(nband),npw(npw),sparsity(sparsity),order(order),eps(eps),maxiter(maxiter) - { -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + public: + DiagoDavPrepare(int nband, int npw, int sparsity, int order, float eps, int maxiter) + : nband(nband), npw(npw), sparsity(sparsity), order(order), eps(eps), maxiter(maxiter) + { +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif - } - - int nband, npw, sparsity, order, maxiter, notconv; - float eps, avg_iter; - int nprocs=1, mypnum=0; - - void CompareEigen(psi::Psi> &phi, float *precondition) - { - //calculate eigenvalues by LAPACK; - float* e_lapack = new float[npw]; - float* ev; - if(mypnum == 0) { lapackEigen(npw, DIAGOTEST::hmatrix_f, e_lapack,DETAILINFO); -} +#endif + } + + int nband, npw, sparsity, order, maxiter, notconv; + float eps, avg_iter; + int nprocs = 1, mypnum = 0; - //do Diago_David::diag() - float* en = new float[npw]; - HSOperatorMock> phm; + void CompareEigen(psi::Psi>& phi, float* precondition) + { + // calculate eigenvalues by LAPACK; + float* e_lapack = new float[npw]; + float* ev; + if (mypnum == 0) + { + lapackEigen(npw, DIAGOTEST::hmatrix_f, e_lapack, DETAILINFO); + } + + // do Diago_David::diag() + float* en = new float[npw]; + HSOperatorMock> phm; -#ifdef __MPI +#ifdef __MPI const hsolver::diag_comm_info comm_info = {POOL_WORLD, mypnum, nprocs}; #else - const hsolver::diag_comm_info comm_info = {mypnum, nprocs}; + const hsolver::diag_comm_info comm_info = {mypnum, nprocs}; #endif - const int dim = phi.get_current_ngk() ; - const int nband = phi.get_nbands(); - const int ld_psi =phi.get_nbasis(); - hsolver::DiagoDavid> dav(precondition, nband, dim, order, comm_info); + const int dim = phi.get_current_ngk(); + const int nband = phi.get_nbands(); + const int ld_psi = phi.get_nbasis(); + hsolver::DiagoDavid> dav(precondition, nband, dim, order, comm_info); - hsolver::DiagoIterAssist>::PW_DIAG_NMAX = maxiter; + hsolver::DiagoIterAssist>::PW_DIAG_NMAX = maxiter; hsolver::DiagoIterAssist>::PW_DIAG_THR = eps; phi.fix_k(0); - float use_time = 0.0; -#ifdef __MPI - float start = 0.0, end = 0.0; - start = MPI_Wtime(); + float use_time = 0.0; +#ifdef __MPI + float start = 0.0, end = 0.0; + start = MPI_Wtime(); #else - clock_t start, end; - start = clock(); -#endif + clock_t start, end; + start = clock(); +#endif - std::vector ethr_band(phi.get_nbands(), eps); - dav.diag(phm, ld_psi, phi.get_pointer(), en, ethr_band, maxiter); - -#ifdef __MPI - end = MPI_Wtime(); - use_time = end - start; -#else - end = clock(); - use_time = (float)(end-start); -#endif - - if(mypnum == 0) - { - if (DETAILINFO) { std::cout<<"diag Run time: "<< use_time << std::endl; -} - for(int i=0;i {}; +class DiagoDavTest : public ::testing::TestWithParam +{ +}; -TEST_P(DiagoDavTest,RandomHamilt) +TEST_P(DiagoDavTest, RandomHamilt) { - DiagoDavPrepare ddp = GetParam(); - if (DETAILINFO&&ddp.mypnum==0) { std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" - << ddp.sparsity << ", eps=" << ddp.eps << std::endl; -} + DiagoDavPrepare ddp = GetParam(); + if (DETAILINFO && ddp.mypnum == 0) + { + std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" << ddp.sparsity + << ", eps=" << ddp.eps << std::endl; + } HPsi> hpsi(ddp.nband, ddp.npw, ddp.sparsity); - DIAGOTEST::hmatrix_f = hpsi.hamilt(); - DIAGOTEST::npw = ddp.npw; - DIAGOTEST::npw_local = new int[ddp.nprocs]; - psi::Psi> psi = hpsi.psi(); - psi::Psi> psi_local; - float* precondition_local; - -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::hmatrix_f = hpsi.hamilt(); + DIAGOTEST::npw = ddp.npw; + DIAGOTEST::npw_local = new int[ddp.nprocs]; + psi::Psi> psi = hpsi.psi(); + psi::Psi> psi_local; + float* precondition_local; + +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_f, DIAGOTEST::hmatrix_local_f); - precondition_local = new float[DIAGOTEST::npw_local[ddp.mypnum]]; + precondition_local = new float[DIAGOTEST::npw_local[ddp.mypnum]]; DIAGOTEST::divide_psi(hpsi.precond(), precondition_local); #else - DIAGOTEST::hmatrix_local_f = DIAGOTEST::hmatrix_f; - DIAGOTEST::npw_local[0] = DIAGOTEST::npw; - psi_local = psi; - precondition_local = new float[DIAGOTEST::npw]; - for(int i=0;i> hmatrix; - std::ifstream ifs; - std::string filename = "H-KPoints-Si2.dat"; - ifs.open(filename); + std::vector> hmatrix; + std::ifstream ifs; + std::string filename = "H-KPoints-Si2.dat"; + ifs.open(filename); // open file and check status if (!ifs.is_open()) { std::cout << "Error opening file " << filename << std::endl; exit(1); } - DIAGOTEST::readh(ifs,hmatrix); - ifs.close(); - DIAGOTEST::hmatrix_f = hmatrix; - int nband = std::max(DIAGOTEST::npw/20,1); + DIAGOTEST::readh(ifs, hmatrix); + ifs.close(); + DIAGOTEST::hmatrix_f = hmatrix; + int nband = std::max(DIAGOTEST::npw / 20, 1); + + DiagoDavPrepare ddp(nband, DIAGOTEST::npw, 0, 2, 1e-5, 500); - DiagoDavPrepare ddp(nband,DIAGOTEST::npw,0,2,1e-5,500); - HPsi> hpsi(nband, DIAGOTEST::npw); - psi::Psi> psi = hpsi.psi(); - DIAGOTEST::npw_local = new int[ddp.nprocs]; - psi::Psi> psi_local; - float* precondition_local; + psi::Psi> psi = hpsi.psi(); + DIAGOTEST::npw_local = new int[ddp.nprocs]; + psi::Psi> psi_local; + float* precondition_local; -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_f, DIAGOTEST::hmatrix_local_f); - precondition_local = new float[DIAGOTEST::npw_local[ddp.mypnum]]; + precondition_local = new float[DIAGOTEST::npw_local[ddp.mypnum]]; DIAGOTEST::divide_psi(hpsi.precond(), precondition_local); #else - DIAGOTEST::hmatrix_local_f = DIAGOTEST::hmatrix_f; - DIAGOTEST::npw_local[0] = DIAGOTEST::npw; - psi_local = psi; - precondition_local = new float[DIAGOTEST::npw]; - for(int i=0;ilisteners(); - if (myrank != 0) { delete listeners.Release(listeners.default_result_printer()); -} + ::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 0; + return 0; } diff --git a/source/source_hsolver/test/diago_david_real_test.cpp b/source/source_hsolver/test/diago_david_real_test.cpp index 3ff4e33bc7b..aabfdac13b7 100644 --- a/source/source_hsolver/test/diago_david_real_test.cpp +++ b/source/source_hsolver/test/diago_david_real_test.cpp @@ -1,22 +1,22 @@ -#include"source_hsolver/diago_david.h" -#include "source_hsolver/diag_comm_info.h" -#include"source_hsolver/diago_iter_assist.h" -#include "source_base/parallel_comm.h" -#include"diago_mock.h" -#include "source_psi/psi.h" -#include"gtest/gtest.h" +#include "diago_mock.h" +#include "mpi.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" -#include"mpi.h" +#include "source_hsolver/diag_comm_info.h" +#include "source_hsolver/diago_david.h" +#include "source_hsolver/diago_iter_assist.h" +#include "source_psi/psi.h" + +#include "gtest/gtest.h" #define CONVTHRESHOLD 1e-0 #define DETAILINFO false - /************************************************ -* unit test of class Diago_David -***********************************************/ + * unit test of class Diago_David + ***********************************************/ /** * Class Diago_David is used to solve the eigenvalues @@ -31,8 +31,8 @@ * */ - //use lapack to calcualte eigenvalue of matrix hm - //NOTE: after finish this function, hm stores the eigen vectors. +// use lapack to calcualte eigenvalue of matrix hm +// NOTE: after finish this function, hm stores the eigen vectors. void lapackEigen(int& npw, std::vector& hm, double* e, bool outtime = false) { int info = 0; @@ -43,28 +43,32 @@ void lapackEigen(int& npw, std::vector& hm, double* e, bool outtime = fa double work_tmp; constexpr int minus_one = -1; - dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, &work_tmp, &minus_one, &info); // get best lwork + dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, &work_tmp, &minus_one, &info); // get best lwork const int lwork = work_tmp; double* work2 = new double[lwork]; dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, &info); end = clock(); - if (info) { std::cout << "ERROR: Lapack solver, info=" << info << std::endl; -} - if (outtime) { std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; -} + if (info) + { + std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + } + if (outtime) + { + std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + } delete[] work2; } class DiagoDavPrepare { -public: - DiagoDavPrepare(int nband, int npw, int sparsity, int order, double eps, int maxiter) : - nband(nband), npw(npw), sparsity(sparsity), order(order), eps(eps), maxiter(maxiter) + public: + DiagoDavPrepare(int nband, int npw, int sparsity, int order, double eps, int maxiter) + : nband(nband), npw(npw), sparsity(sparsity), order(order), eps(eps), maxiter(maxiter) { -#ifdef __MPI +#ifdef __MPI MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +#endif } int nband, npw, sparsity, order, maxiter, notconv; @@ -73,23 +77,25 @@ class DiagoDavPrepare void CompareEigen(psi::Psi& phi, double* precondition) { - //calculate eigenvalues by LAPACK; + // calculate eigenvalues by LAPACK; double* e_lapack = new double[npw]; double* ev; - if (mypnum == 0) { lapackEigen(npw, DIAGOTEST::hmatrix_d, e_lapack, DETAILINFO); -} + if (mypnum == 0) + { + lapackEigen(npw, DIAGOTEST::hmatrix_d, e_lapack, DETAILINFO); + } - //do Diago_David::diag() + // do Diago_David::diag() double* en = new double[npw]; HSOperatorMock phm; -#ifdef __MPI +#ifdef __MPI const hsolver::diag_comm_info comm_info = {POOL_WORLD, mypnum, nprocs}; #else const hsolver::diag_comm_info comm_info = {mypnum, nprocs}; #endif - const int dim = phi.get_current_ngk(); + const int dim = phi.get_current_ngk(); const int nband = phi.get_nbands(); const int ld_psi = phi.get_nbasis(); hsolver::DiagoDavid dav(precondition, nband, dim, order, comm_info); @@ -99,31 +105,32 @@ class DiagoDavPrepare phi.fix_k(0); double use_time = 0.0; -#ifdef __MPI +#ifdef __MPI double start = 0.0, end = 0.0; start = MPI_Wtime(); #else clock_t start, end; start = clock(); -#endif +#endif - std::vector ethr_band(phi.get_nbands(), eps); dav.diag(phm, ld_psi, phi.get_pointer(), en, ethr_band, maxiter); -#ifdef __MPI +#ifdef __MPI end = MPI_Wtime(); use_time = end - start; -#else +#else end = clock(); use_time = (double)(end - start); -#endif +#endif if (mypnum == 0) { - if (DETAILINFO) { std::cout << "diag Run time: " << use_time << std::endl; -} - for (int i = 0;i < nband;i++) + if (DETAILINFO) + { + std::cout << "diag Run time: " << use_time << std::endl; + } + for (int i = 0; i < nband; i++) { EXPECT_NEAR(en[i], e_lapack[i], CONVTHRESHOLD); } @@ -133,14 +140,18 @@ class DiagoDavPrepare } }; -class DiagoDavTest : public ::testing::TestWithParam {}; +class DiagoDavTest : public ::testing::TestWithParam +{ +}; TEST_P(DiagoDavTest, RandomHamilt) { DiagoDavPrepare ddp = GetParam(); - if (DETAILINFO && ddp.mypnum == 0) { std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" - << ddp.sparsity << ", eps=" << ddp.eps << std::endl; -} + if (DETAILINFO && ddp.mypnum == 0) + { + std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" << ddp.sparsity + << ", eps=" << ddp.eps << std::endl; + } HPsi hpsi(ddp.nband, ddp.npw, ddp.sparsity); DIAGOTEST::hmatrix_d = hpsi.hamilt(); @@ -150,7 +161,7 @@ TEST_P(DiagoDavTest, RandomHamilt) psi::Psi psi_local; double* precondition_local; -#ifdef __MPI +#ifdef __MPI DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_d, DIAGOTEST::hmatrix_local_d); precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; @@ -160,7 +171,8 @@ TEST_P(DiagoDavTest, RandomHamilt) 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] = (hpsi.precond())[i]; + for (int i = 0; i < DIAGOTEST::npw; i++) + precondition_local[i] = (hpsi.precond())[i]; #endif ddp.CompareEigen(psi_local, precondition_local); @@ -168,14 +180,15 @@ TEST_P(DiagoDavTest, RandomHamilt) delete[] precondition_local; } - -INSTANTIATE_TEST_SUITE_P(VerifyDiag, DiagoDavTest, ::testing::Values( - //DiagoDavPrepare(int nband, int npw, int sparsity, int order,double eps,int maxiter) - DiagoDavPrepare(10, 100, 0, 4, 1e-5, 500), - DiagoDavPrepare(20, 500, 7, 4, 1e-5, 500) - //DiagoDavPrepare(50,1000,8,4,1e-5,500) - //DiagoDavPrepare(20,2000,8,4,1e-5,500) -)); +INSTANTIATE_TEST_SUITE_P(VerifyDiag, + DiagoDavTest, + ::testing::Values( + // DiagoDavPrepare(int nband, int npw, int sparsity, int order,double eps,int maxiter) + DiagoDavPrepare(10, 100, 0, 4, 1e-5, 500), + DiagoDavPrepare(20, 500, 7, 4, 1e-5, 500) + // DiagoDavPrepare(50,1000,8,4,1e-5,500) + // DiagoDavPrepare(20,2000,8,4,1e-5,500) + )); TEST(DiagoDavRealSystemTest, dataH) { @@ -202,7 +215,7 @@ TEST(DiagoDavRealSystemTest, dataH) psi::Psi psi_local; double* precondition_local; -#ifdef __MPI +#ifdef __MPI DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_d, DIAGOTEST::hmatrix_local_d); precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; @@ -212,7 +225,8 @@ TEST(DiagoDavRealSystemTest, dataH) 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] = (hpsi.precond())[i]; + for (int i = 0; i < DIAGOTEST::npw; i++) + precondition_local[i] = (hpsi.precond())[i]; #endif ddp.CompareEigen(psi_local, precondition_local); @@ -235,8 +249,10 @@ int main(int argc, char** argv) testing::InitGoogleTest(&argc, argv); ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners(); - if (myrank != 0) { delete listeners.Release(listeners.default_result_printer()); -} + if (myrank != 0) + { + delete listeners.Release(listeners.default_result_printer()); + } int result = RUN_ALL_TESTS(); if (myrank == 0 && result != 0) diff --git a/source/source_hsolver/test/diago_david_test.cpp b/source/source_hsolver/test/diago_david_test.cpp index 7d4b8cbfea0..732bc479e70 100644 --- a/source/source_hsolver/test/diago_david_test.cpp +++ b/source/source_hsolver/test/diago_david_test.cpp @@ -1,22 +1,23 @@ -#include"source_hsolver/diago_david.h" -#include "source_hsolver/diag_comm_info.h" -#include"source_hsolver/diago_iter_assist.h" -#include "source_base/parallel_comm.h" -#include"diago_mock.h" -#include "source_psi/psi.h" -#include"gtest/gtest.h" +#include "source_hsolver/diago_david.h" + +#include "diago_mock.h" +#include "mpi.h" #include "source_base/inverse_matrix.h" #include "source_base/module_external/lapack_connector.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" -#include"mpi.h" +#include "source_hsolver/diag_comm_info.h" +#include "source_hsolver/diago_iter_assist.h" +#include "source_psi/psi.h" + +#include "gtest/gtest.h" #define CONVTHRESHOLD 1e-3 #define DETAILINFO false - /************************************************ -* unit test of class Diago_David -***********************************************/ + * unit test of class Diago_David + ***********************************************/ /** * Class Diago_David is used to solve the eigenvalues @@ -26,230 +27,243 @@ * - the hamilt matrix (npw=100,500,1000) produced by random with sparsity of 50% * - the hamilt matrix (npw=100,500,1000) produced by random with sparsity of 0% * - the hamilt matrix read from "data-H" - * + * * The test is passed when the eignvalues are closed to these calculated by LAPACK. - * + * */ -//use lapack to calcualte eigenvalue of matrix hm -//NOTE: after finish this function, hm stores the eigen vectors. - +// use lapack to calcualte eigenvalue of matrix hm +// NOTE: after finish this function, hm stores the eigen vectors. void lapackEigen(int& npw, std::vector>& hm, double* e, bool outtime = false) { - int lwork = 2 * npw; - std::complex *work2= new std::complex[lwork]; - double* rwork = new double[3*npw-2]; - int info = 0; - - auto tmp = hm; - - clock_t start,end; - start = clock(); - char tmp_c1 = 'V', tmp_c2 = 'U'; - zheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); - end = clock(); - if(info) { std::cout << "ERROR: Lapack solver, info=" << info <* work2 = new std::complex[lwork]; + double* rwork = new double[3 * npw - 2]; + int info = 0; + + auto tmp = hm; + + clock_t start, end; + start = clock(); + char tmp_c1 = 'V', tmp_c2 = 'U'; + zheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); + end = clock(); + if (info) + { + std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + } + if (outtime) + { + std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + } - delete [] rwork; - delete [] work2; + delete[] rwork; + delete[] work2; } -class DiagoDavPrepare +class DiagoDavPrepare { -public: - DiagoDavPrepare(int nband, int npw, int sparsity, int order,double eps,int maxiter): - nband(nband),npw(npw),sparsity(sparsity),order(order),eps(eps),maxiter(maxiter) - { -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + public: + DiagoDavPrepare(int nband, int npw, int sparsity, int order, double eps, int maxiter) + : nband(nband), npw(npw), sparsity(sparsity), order(order), eps(eps), maxiter(maxiter) + { +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif - } - - int nband, npw, sparsity, order, maxiter, notconv; - double eps, avg_iter; - int nprocs=1, mypnum=0; - - void CompareEigen(psi::Psi> &phi, double *precondition) - { - //calculate eigenvalues by LAPACK; - double* e_lapack = new double[npw]; - double* ev = nullptr; - if(mypnum == 0) - { - lapackEigen(npw, DIAGOTEST::hmatrix, e_lapack,DETAILINFO); - } - - //do Diago_David::diag() - double* en = new double[npw]; - HSOperatorMock> phm; - -#ifdef __MPI +#endif + } + + int nband, npw, sparsity, order, maxiter, notconv; + double eps, avg_iter; + int nprocs = 1, mypnum = 0; + + void CompareEigen(psi::Psi>& phi, double* precondition) + { + // calculate eigenvalues by LAPACK; + double* e_lapack = new double[npw]; + double* ev = nullptr; + if (mypnum == 0) + { + lapackEigen(npw, DIAGOTEST::hmatrix, e_lapack, DETAILINFO); + } + + // do Diago_David::diag() + double* en = new double[npw]; + HSOperatorMock> phm; + +#ifdef __MPI const hsolver::diag_comm_info comm_info = {POOL_WORLD, mypnum, nprocs}; #else const hsolver::diag_comm_info comm_info = {mypnum, nprocs}; #endif - const int dim = phi.get_current_ngk(); - const int nband = phi.get_nbands(); - const int ld_psi = phi.get_nbasis(); - hsolver::DiagoDavid> dav(precondition, nband, dim, order, comm_info); + const int dim = phi.get_current_ngk(); + const int nband = phi.get_nbands(); + const int ld_psi = phi.get_nbasis(); + hsolver::DiagoDavid> dav(precondition, nband, dim, order, comm_info); - hsolver::DiagoIterAssist>::PW_DIAG_NMAX = maxiter; + hsolver::DiagoIterAssist>::PW_DIAG_NMAX = maxiter; hsolver::DiagoIterAssist>::PW_DIAG_THR = eps; phi.fix_k(0); - double use_time = 0.0; -#ifdef __MPI - double start = 0.0, end = 0.0; - start = MPI_Wtime(); + double use_time = 0.0; +#ifdef __MPI + double start = 0.0, end = 0.0; + start = MPI_Wtime(); #else - clock_t start, end; - start = clock(); -#endif + clock_t start, end; + start = clock(); +#endif - std::vector ethr_band(phi.get_nbands(), eps); - dav.diag(phm, ld_psi, phi.get_pointer(), en, ethr_band, maxiter); - -#ifdef __MPI - end = MPI_Wtime(); - use_time = end - start; -#else - end = clock(); - use_time = (double)(end-start); -#endif - - if(mypnum == 0) - { - if (DETAILINFO) { std::cout<<"diag Run time: "<< use_time << std::endl; -} - for(int i=0;i {}; +class DiagoDavTest : public ::testing::TestWithParam +{ +}; -TEST_P(DiagoDavTest,RandomHamilt) +TEST_P(DiagoDavTest, RandomHamilt) { - DiagoDavPrepare ddp = GetParam(); - if (DETAILINFO&&ddp.mypnum==0) { std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" - << ddp.sparsity << ", eps=" << ddp.eps << std::endl; -} + DiagoDavPrepare ddp = GetParam(); + if (DETAILINFO && ddp.mypnum == 0) + { + std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" << ddp.sparsity + << ", eps=" << ddp.eps << std::endl; + } HPsi> hpsi(ddp.nband, ddp.npw, ddp.sparsity); - DIAGOTEST::hmatrix = hpsi.hamilt(); - DIAGOTEST::npw = ddp.npw; - DIAGOTEST::npw_local = new int[ddp.nprocs]; - psi::Psi> psi = hpsi.psi(); - psi::Psi> psi_local; - double* precondition_local; - -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::hmatrix = hpsi.hamilt(); + DIAGOTEST::npw = ddp.npw; + DIAGOTEST::npw_local = new int[ddp.nprocs]; + psi::Psi> psi = hpsi.psi(); + psi::Psi> psi_local; + double* precondition_local; + +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix, DIAGOTEST::hmatrix_local); - precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; - DIAGOTEST::divide_psi(hpsi.precond(),precondition_local); + precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; + DIAGOTEST::divide_psi(hpsi.precond(), 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> hmatrix; - std::ifstream ifs; - std::string filename = "H-KPoints-Si2.dat"; - ifs.open(filename); + std::vector> hmatrix; + std::ifstream ifs; + std::string filename = "H-KPoints-Si2.dat"; + ifs.open(filename); // open file and check status if (!ifs.is_open()) { std::cout << "Error opening file " << filename << std::endl; exit(1); } - DIAGOTEST::readh(ifs,hmatrix); - ifs.close(); - DIAGOTEST::hmatrix = hmatrix; - int nband = std::max(DIAGOTEST::npw/20,1); + DIAGOTEST::readh(ifs, hmatrix); + ifs.close(); + DIAGOTEST::hmatrix = hmatrix; + int nband = std::max(DIAGOTEST::npw / 20, 1); + + DiagoDavPrepare ddp(nband, DIAGOTEST::npw, 0, 2, 1e-5, 500); - DiagoDavPrepare ddp(nband,DIAGOTEST::npw,0,2,1e-5,500); - HPsi> hpsi(nband, DIAGOTEST::npw); - psi::Psi> psi = hpsi.psi(); - DIAGOTEST::npw_local = new int[ddp.nprocs]; - psi::Psi> psi_local; - double* precondition_local; + psi::Psi> psi = hpsi.psi(); + DIAGOTEST::npw_local = new int[ddp.nprocs]; + psi::Psi> psi_local; + double* precondition_local; -#ifdef __MPI - DIAGOTEST::cal_division(DIAGOTEST::npw); +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix, DIAGOTEST::hmatrix_local); - precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; - DIAGOTEST::divide_psi(hpsi.precond(),precondition_local); + precondition_local = new double[DIAGOTEST::npw_local[ddp.mypnum]]; + DIAGOTEST::divide_psi(hpsi.precond(), 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;ilisteners(); - if (myrank != 0) { delete listeners.Release(listeners.default_result_printer()); -} + ::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 0; + return 0; } diff --git a/source/source_hsolver/test/diago_elpa_utils.h b/source/source_hsolver/test/diago_elpa_utils.h index 843c53b6a51..c1fab37c259 100644 --- a/source/source_hsolver/test/diago_elpa_utils.h +++ b/source/source_hsolver/test/diago_elpa_utils.h @@ -1,6 +1,6 @@ +#include "mpi.h" #include "source_base/module_external/blacs_connector.h" #include "source_base/module_external/lapack_connector.h" -#include "mpi.h" #include #include @@ -17,7 +17,7 @@ namespace LCAO_DIAGO_TEST { -void process_2d(int &nprows, int &npcols, int &myprow, int &mypcol, int &icontxt) // map the processes to 2D matrix +void process_2d(int& nprows, int& npcols, int& myprow, int& mypcol, int& icontxt) // map the processes to 2D matrix { int nprocs, mypnum; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); @@ -30,7 +30,7 @@ void process_2d(int &nprows, int &npcols, int &myprow, int &mypcol, int &icontxt } nprows = nprocs / npcols; - int *usermap = new int[nprocs]; + int* usermap = new int[nprocs]; for (int i = 0; i < nprows; i++) { for (int j = 0; j < npcols; j++) @@ -57,7 +57,7 @@ inline int na_rc(int n, int nblk, int nprows, int myprow) // calculate the row/c // distribute hmatrix to each process // hmatrix is row-first, and new_matrix is column-first template -void distribute_data(T *hmatrix, T *new_matrix, int &nFull, int &nblk, int &na_rows, int &na_cols, int &icontxt) +void distribute_data(T* hmatrix, T* new_matrix, int& nFull, int& nblk, int& na_rows, int& na_cols, int& icontxt) { int mypnum; MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); @@ -65,12 +65,12 @@ void distribute_data(T *hmatrix, T *new_matrix, int &nFull, int &nblk, int &na_r Cblacs_gridinfo(icontxt, &nprows, &npcols, &myprow, &mypcol); int SENDPROW = 0, SENDPCOL = 0, tag = 0; - T *tmp = new T[nblk]; + T* tmp = new T[nblk]; // do iteration for matrix, distribute old_matrix to each process, pass a block each time for (int row = 0; row < nFull; row++) { - int recv_prow = (row / nblk) % nprows; // the row number of recive process + int recv_prow = (row / nblk) % nprows; // the row number of recive process int nm_row = ((row / nblk) / nprows) * nblk + row % nblk; // row number of block in new_matrix for (int col = 0; col < nFull; col += nblk) { @@ -121,7 +121,8 @@ void distribute_data(T *hmatrix, T *new_matrix, int &nFull, int &nblk, int &na_r delete[] tmp; } -template inline void print_matrix(std::ofstream &fp, T *matrix, int &nrow, int &ncol, bool row_first) +template +inline void print_matrix(std::ofstream& fp, T* matrix, int& nrow, int& ncol, bool row_first) { int coef_row = row_first ? ncol : 1; int coef_col = row_first ? 1 : nrow; @@ -135,17 +136,18 @@ template inline void print_matrix(std::ofstream &fp, T *matrix, int &n } } -template bool read_hs(std::string fname, T &matrix) +template +bool read_hs(std::string fname, T& matrix) { int ndim; std::ifstream inf(fname); - if(! inf.is_open()) + if (!inf.is_open()) { std::cout << "Error: open file " << fname << " failed, skip!" << std::endl; return false; } inf >> ndim; - matrix.resize(ndim*ndim); + matrix.resize(ndim * ndim); for (int i = 0; i < ndim; i++) { for (int j = i; j < ndim; j++) @@ -161,16 +163,16 @@ template bool read_hs(std::string fname, T &matrix) return true; } -void lapack_diago(double *hmatrix, double *smatrix, double *e, int &nFull) +void lapack_diago(double* hmatrix, double* smatrix, double* e, int& nFull) { - const int itype = 1; // solve A*X=(lambda)*B*X + const int itype = 1; // solve A*X=(lambda)*B*X const char jobz = 'V'; // 'N':only calc eigenvalue, 'V': eigenvalues and eigenvectors const char uplo = 'U'; // Upper triangles int lwork = (nFull + 2) * nFull, info = 0; - double *ev = new double[lwork]; + double* ev = new double[lwork]; - double *a = new double[nFull * nFull]; - double *b = new double[nFull * nFull]; + double* a = new double[nFull * nFull]; + double* b = new double[nFull * nFull]; for (int i = 0; i < nFull * nFull; i++) { a[i] = hmatrix[i]; @@ -189,17 +191,17 @@ void lapack_diago(double *hmatrix, double *smatrix, double *e, int &nFull) delete[] ev; } -void lapack_diago(std::complex *hmatrix, std::complex *smatrix, double *e, int &nFull) +void lapack_diago(std::complex* hmatrix, std::complex* smatrix, double* e, int& nFull) { - const int itype = 1; // solve A*X=(lambda)*B*X + const int itype = 1; // solve A*X=(lambda)*B*X const char jobz = 'V'; // 'N':only calc eigenvalue, 'V': eigenvalues and eigenvectors const char uplo = 'U'; // Upper triangles int lwork = (nFull + 1) * nFull, info = 0; - double *rwork = new double[3 * nFull - 2]; - std::complex *ev = new std::complex[lwork]; + double* rwork = new double[3 * nFull - 2]; + std::complex* ev = new std::complex[lwork]; - std::complex *a = new std::complex[nFull * nFull]; - std::complex *b = new std::complex[nFull * nFull]; + std::complex* a = new std::complex[nFull * nFull]; + std::complex* b = new std::complex[nFull * nFull]; for (int i = 0; i < nFull * nFull; i++) { a[i] = hmatrix[i]; diff --git a/source/source_hsolver/test/diago_lapack_test.cpp b/source/source_hsolver/test/diago_lapack_test.cpp index f5dead8b518..bc91b4524d9 100644 --- a/source/source_hsolver/test/diago_lapack_test.cpp +++ b/source/source_hsolver/test/diago_lapack_test.cpp @@ -3,13 +3,14 @@ // Remove some useless functions and dependencies. Serialized the full code // and refactored some function. +#include "source_hsolver/diago_lapack.h" + +#include "source_base/matrix_block.h" + #include "gtest/gtest.h" #include #include -#include "source_base/matrix_block.h" -#include "source_hsolver/diago_lapack.h" - #define PASSTHRESHOLD 1e-5 #define DETAILINFO false #define PRINT_HS false @@ -37,17 +38,18 @@ class HamiltTEST // The serialized version of functions from diago_elpa_utils -template bool read_hs(std::string fname, T &matrix) +template +bool read_hs(std::string fname, T& matrix) { int ndim; std::ifstream inf(fname); - if(! inf.is_open()) + if (!inf.is_open()) { std::cout << "Error: open file " << fname << " failed, skip!" << std::endl; return false; } inf >> ndim; - matrix.resize(ndim*ndim); + matrix.resize(ndim * ndim); for (int i = 0; i < ndim; i++) { for (int j = i; j < ndim; j++) @@ -63,11 +65,11 @@ template bool read_hs(std::string fname, T &matrix) return true; } -bool read_solution(std::string fname, std::vector &result) +bool read_solution(std::string fname, std::vector& result) { int ndim; std::ifstream inf(fname); - if(! inf.is_open()) + if (!inf.is_open()) { std::cout << "Error: open file " << fname << " failed, skip!" << std::endl; return false; @@ -81,8 +83,8 @@ bool read_solution(std::string fname, std::vector &result) return true; } - -template inline void print_matrix(std::ofstream &fp, T *matrix, int &nrow, int &ncol, bool row_first) +template +inline void print_matrix(std::ofstream& fp, T* matrix, int& nrow, int& ncol, bool row_first) { int coef_row = row_first ? ncol : 1; int coef_col = row_first ? 1 : nrow; @@ -96,11 +98,10 @@ template inline void print_matrix(std::ofstream &fp, T *matrix, int &n } } - template class DiagoLapackPrepare { - public: + public: DiagoLapackPrepare(int nlocal, int nbands, int nb2d, @@ -108,8 +109,8 @@ class DiagoLapackPrepare std::string hfname, std::string sfname, std::string solutionfname) - : nlocal(nlocal), nbands(nbands), nb2d(nb2d), sparsity(sparsity), hfname(hfname), - sfname(sfname), solutionfname(solutionfname) + : nlocal(nlocal), nbands(nbands), nb2d(nb2d), sparsity(sparsity), hfname(hfname), sfname(sfname), + solutionfname(solutionfname) { // dh = new hsolver::DiagoLapack; } @@ -141,9 +142,10 @@ class DiagoLapackPrepare } nlocal = hdim; nbands = nlocal / 2; - if (readhfile && readsfile) { + if (readhfile && readsfile) + { return true; - } + } return false; } @@ -157,7 +159,8 @@ class DiagoLapackPrepare void print_hs() { - if (!PRINT_HS) { + if (!PRINT_HS) + { return; } std::ofstream fp("hmatrix.dat"); @@ -211,7 +214,8 @@ class DiagoLapackPrepare maxerror = error; iindex = i; } - if (error > PASSTHRESHOLD) { + if (error > PASSTHRESHOLD) + { pass = false; } } @@ -248,13 +252,22 @@ TEST_P(DiagoLapackGammaOnlyTest, LCAO) EXPECT_TRUE(pass) << out_info.str(); } -INSTANTIATE_TEST_SUITE_P( - DiagoLapackTest, - DiagoLapackGammaOnlyTest, - ::testing::Values( - DiagoLapackPrepare(0, 0, 1, 0, "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat", "GammaOnly-Si2-Solution.dat"), - DiagoLapackPrepare(0, 0, 32, 0, "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat", "GammaOnly-Si64-Solution.dat"))); - +INSTANTIATE_TEST_SUITE_P(DiagoLapackTest, + DiagoLapackGammaOnlyTest, + ::testing::Values(DiagoLapackPrepare(0, + 0, + 1, + 0, + "H-GammaOnly-Si2.dat", + "S-GammaOnly-Si2.dat", + "GammaOnly-Si2-Solution.dat"), + DiagoLapackPrepare(0, + 0, + 32, + 0, + "H-GammaOnly-Si64.dat", + "S-GammaOnly-Si64.dat", + "GammaOnly-Si64-Solution.dat"))); class DiagoLapackKPointsTest : public ::testing::TestWithParam>> { @@ -271,12 +284,22 @@ TEST_P(DiagoLapackKPointsTest, LCAO) EXPECT_TRUE(pass) << out_info.str(); } -INSTANTIATE_TEST_SUITE_P( - DiagoLapackTest, - DiagoLapackKPointsTest, - ::testing::Values( - DiagoLapackPrepare>(0, 0, 1, 0, "H-KPoints-Si2.dat", "S-KPoints-Si2.dat", "KPoints-Si2-Solution.dat"), - DiagoLapackPrepare>(0, 0, 32, 0, "H-KPoints-Si64.dat", "S-KPoints-Si64.dat", "KPoints-Si64-Solution.dat"))); +INSTANTIATE_TEST_SUITE_P(DiagoLapackTest, + DiagoLapackKPointsTest, + ::testing::Values(DiagoLapackPrepare>(0, + 0, + 1, + 0, + "H-KPoints-Si2.dat", + "S-KPoints-Si2.dat", + "KPoints-Si2-Solution.dat"), + DiagoLapackPrepare>(0, + 0, + 32, + 0, + "H-KPoints-Si64.dat", + "S-KPoints-Si64.dat", + "KPoints-Si64-Solution.dat"))); int main(int argc, char** argv) { diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp index bab49060153..6edea073b9c 100644 --- a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp +++ b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp @@ -1,7 +1,7 @@ +#include "mpi.h" #include "source_base/matrix_block.h" #include "source_hsolver/diago_scalapack.h" #include "source_hsolver/test/diago_elpa_utils.h" -#include "mpi.h" #include "string.h" #include "gtest/gtest.h" @@ -68,10 +68,10 @@ class DiagoPrepare MPI_Comm_rank(MPI_COMM_WORLD, &myrank); if (ks_solver == "scalapack_gvx") - ;//dh = new hsolver::DiagoScalapack; + ; // dh = new hsolver::DiagoScalapack; #ifdef __CUDA else if (ks_solver == "cusolver") - ;//dh = new hsolver::DiagoCusolver; + ; // dh = new hsolver::DiagoCusolver; #endif else { @@ -239,13 +239,13 @@ class DiagoPrepare hsolver::DiagoScalapack dh(nlocal, nbands); dh.diag(h_mat, s_mat, psi, e_solver.data()); } - #ifdef __CUDA +#ifdef __CUDA else if (ks_solver == "cusolver") - { - hsolver::DiagoCusolver dh(nlocal, nbands); - dh.diag(h_mat, s_mat, psi, e_solver.data()); - } - #endif + { + hsolver::DiagoCusolver dh(nlocal, nbands); + dh.diag(h_mat, s_mat, psi, e_solver.data()); + } +#endif // dh->diag(&hmtest, psi, e_solver.data()); } endtime = MPI_Wtime(); diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp index 18b4a30fb0c..32f1c13792c 100644 --- a/source/source_hsolver/test/diago_lcao_test.cpp +++ b/source/source_hsolver/test/diago_lcao_test.cpp @@ -1,8 +1,8 @@ +#include "mpi.h" #include "source_base/matrix_block.h" -#include "source_hsolver/diago_scalapack.h" #include "source_hsolver/diago_lapack.h" +#include "source_hsolver/diago_scalapack.h" #include "source_hsolver/test/diago_elpa_utils.h" -#include "mpi.h" #include "string.h" #include "gtest/gtest.h" @@ -69,7 +69,7 @@ class DiagoPrepare if (ks_solver == "scalapack_gvx") ; -// dh = new hsolver::DiagoScalapack; + // dh = new hsolver::DiagoScalapack; else if (ks_solver == "lapack") ; #ifdef __ELPA @@ -269,13 +269,13 @@ class DiagoPrepare hsolver::DiagoLapack la(nlocal, nbands); la.diag(h_mat, s_mat, psi, e_solver.data()); } - #ifdef __ELPA +#ifdef __ELPA else if (ks_solver == "genelpa") { hsolver::DiagoElpa dh(nlocal, nbands); dh.diag(h_mat, s_mat, psi, e_solver.data()); } - #endif +#endif // dh.diag(&hmtest, psi, e_solver.data()); } endtime = MPI_Wtime(); @@ -371,7 +371,7 @@ INSTANTIATE_TEST_SUITE_P( // DiagoPrepare(0, 0, 32, 0, "scalapack_gvx", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat"), DiagoPrepare(0, 0, 1, 0, "lapack", "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat") // DiagoPrepare(0, 0, 32, 0, "lapack", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat") - )); + )); class DiagoKPointsTest : public ::testing::TestWithParam>> { @@ -419,8 +419,7 @@ class DiagoElpaNativeUpperTest : public ::testing::TestWithParam TEST_P(DiagoElpaNativeUpperTest, PreservesInputsAndIgnoresLowerTriangle) { - DiagoPrepare> dp(0, 0, GetParam(), 0, "genelpa", - "H-KPoints-Si2.dat", "S-KPoints-Si2.dat"); + DiagoPrepare> dp(0, 0, GetParam(), 0, "genelpa", "H-KPoints-Si2.dat", "S-KPoints-Si2.dat"); ASSERT_TRUE(dp.produce_HS()); dp.diago_native_with_poisoned_lower(); if (dp.myrank == 0) @@ -498,12 +497,7 @@ TEST(DiagoElpaComplexTest, FallbackUsesAuthoritativeUpperTriangle) dp.pb2d(); dp.distribute_data(); int decomposed_state = 3; - ELPA_Solver solver(false, - MPI_COMM_WORLD, - dp.nbands, - dp.hmtest.nrow, - dp.hmtest.ncol, - dp.hmtest.desc); + ELPA_Solver solver(false, MPI_COMM_WORLD, dp.nbands, dp.hmtest.nrow, dp.hmtest.ncol, dp.hmtest.desc); solver.generalized_eigenvector(dp.h_local.data(), dp.s_local.data(), decomposed_state, diff --git a/source/source_hsolver/test/diago_mock.h b/source/source_hsolver/test/diago_mock.h index cc4e0853e7b..4b25df3fe72 100644 --- a/source/source_hsolver/test/diago_mock.h +++ b/source/source_hsolver/test/diago_mock.h @@ -1,4 +1,3 @@ -#include #include "mpi.h" #include "source_base/parallel_reduce.h" #include "source_hsolver/hs_operator.h" @@ -8,243 +7,292 @@ #include #include #include +#include #include namespace DIAGOTEST { - std::vector hmatrix_d; - std::vector hmatrix_local_d; - std::vector> hmatrix; - std::vector> hmatrix_local; - std::vector> hmatrix_f; - std::vector> hmatrix_local_f; - int h_nr; - int h_nc; - int npw; - int* npw_local; //number of plane wave distributed to each process +std::vector hmatrix_d; +std::vector hmatrix_local_d; +std::vector> hmatrix; +std::vector> hmatrix_local; +std::vector> hmatrix_f; +std::vector> hmatrix_local_f; +int h_nr; +int h_nc; +int npw; +int* npw_local; // number of plane wave distributed to each process - void readh(std::ifstream& inf, std::vector& hm) +void readh(std::ifstream& inf, std::vector& hm) +{ + int npw_; + inf >> npw_; + DIAGOTEST::npw = npw_; + hm.resize(npw * npw); + h_nr = npw; + h_nc = npw; + for (int i = 0; i < npw; i++) { - int npw_; - inf >> npw_; - DIAGOTEST::npw = npw_; - hm.resize(npw * npw); - h_nr = npw; - h_nc = npw; - for (int i = 0;i < npw;i++) + for (int j = i; j < npw; j++) { - for (int j = i;j < npw;j++) + inf >> hm[i * npw + j]; + if (i != j) { - inf >> hm[i * npw + j]; - if (i != j) { hm[j * npw + i] = hm[i * npw + j]; } + hm[j * npw + i] = hm[i * npw + j]; } } } +} - void readh(std::ifstream& inf, std::vector>& hm) +void readh(std::ifstream& inf, std::vector>& hm) +{ + int npw_; + inf >> npw_; + DIAGOTEST::npw = npw_; + hm.resize(npw * npw); + h_nr = npw; + h_nc = npw; + for (int i = 0; i < npw; i++) { - int npw_; - inf >> npw_; - DIAGOTEST::npw = npw_; - hm.resize(npw * npw); - h_nr = npw; - h_nc = npw; - for(int i=0;i> hm[i * npw + j]; + if (i != j) { - inf >> hm[i * npw + j]; - if (i != j) {hm[j * npw + i] = conj(hm[i * npw + j]);} + hm[j * npw + i] = conj(hm[i * npw + j]); } - double real = hm[i * npw + i].real(); - hm[i * npw + i] = std::complex {real,0.0}; } + double real = hm[i * npw + i].real(); + hm[i * npw + i] = std::complex{real, 0.0}; } +} - void readh(std::ifstream &inf, std::vector> &hm) +void readh(std::ifstream& inf, std::vector>& hm) +{ + int npw_; + inf >> npw_; + DIAGOTEST::npw = npw_; + hm.resize(npw * npw); + h_nr = npw; + h_nc = npw; + for (int i = 0; i < npw; i++) { - int npw_; - inf >> npw_; - DIAGOTEST::npw = npw_; - hm.resize(npw * npw); - h_nr = npw; - h_nc = npw; - for(int i=0;i> hm[i * npw + j]; + if (i != j) { - inf >> hm[i * npw + j]; - if (i != j) {hm[j * npw + i] = conj(hm[i * npw + j]);} + hm[j * npw + i] = conj(hm[i * npw + j]); } - float real = hm[i * npw + i].real(); - hm[i * npw + i] = std::complex {real,0.0}; } + float real = hm[i * npw + i].real(); + hm[i * npw + i] = std::complex{real, 0.0}; } +} - template - void divide_psi(T *psi, T *psi_local) - { - int nprocs=1, mypnum=0; -#ifdef __MPI - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); -#endif +template +void divide_psi(T* psi, T* psi_local) +{ + int nprocs = 1, mypnum = 0; +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); +#endif - if(mypnum == 0) - { - for(int j=0;j::value) { - if (std::is_same::value) { - MPI_Send(&(psi[start_point]),npw_local[j],MPI_DOUBLE,j,0,MPI_COMM_WORLD); - } else if(std::is_same>::value) { - MPI_Send(&(psi[start_point]),npw_local[j],MPI_DOUBLE_COMPLEX,j,0,MPI_COMM_WORLD); - } else if (std::is_same::value) { - MPI_Send(&(psi[start_point]), npw_local[j], MPI_FLOAT, j, 0, MPI_COMM_WORLD); - } else if (std::is_same>::value) { - MPI_Send(&(psi[start_point]), npw_local[j], MPI_C_FLOAT_COMPLEX, j, 0, MPI_COMM_WORLD); -} - start_point += npw_local[j]; + MPI_Send(&(psi[start_point]), npw_local[j], MPI_DOUBLE, j, 0, MPI_COMM_WORLD); + } + else if (std::is_same>::value) + { + MPI_Send(&(psi[start_point]), npw_local[j], MPI_DOUBLE_COMPLEX, j, 0, MPI_COMM_WORLD); + } + else if (std::is_same::value) + { + MPI_Send(&(psi[start_point]), npw_local[j], MPI_FLOAT, j, 0, MPI_COMM_WORLD); + } + else if (std::is_same>::value) + { + MPI_Send(&(psi[start_point]), npw_local[j], MPI_C_FLOAT_COMPLEX, j, 0, MPI_COMM_WORLD); } + start_point += npw_local[j]; } - else + } + else + { + int recv_len = mypnum < (npw % nprocs) ? npw / nprocs + 1 : npw / nprocs; + if (std::is_same::value) { - int recv_len = mypnum < (npw%nprocs) ? npw/nprocs + 1 : npw/nprocs; - if (std::is_same::value) { - MPI_Recv(psi_local, npw_local[mypnum],MPI_DOUBLE,0,0,MPI_COMM_WORLD, MPI_STATUS_IGNORE); - } else if(std::is_same>::value) { - MPI_Recv(psi_local, npw_local[mypnum], MPI_DOUBLE_COMPLEX, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - } else if (std::is_same::value) { - MPI_Recv(psi_local, npw_local[mypnum], MPI_FLOAT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - } else if (std::is_same>::value) { - MPI_Recv(psi_local, npw_local[mypnum], MPI_C_FLOAT_COMPLEX, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); -} -#endif - } + MPI_Recv(psi_local, npw_local[mypnum], MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } + else if (std::is_same>::value) + { + MPI_Recv(psi_local, npw_local[mypnum], MPI_DOUBLE_COMPLEX, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } + else if (std::is_same::value) + { + MPI_Recv(psi_local, npw_local[mypnum], MPI_FLOAT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } + else if (std::is_same>::value) + { + MPI_Recv(psi_local, npw_local[mypnum], MPI_C_FLOAT_COMPLEX, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + } +#endif } - template void divide_psi(double* psi, double* psi_local); - template void divide_psi>(std::complex* psi, std::complex* psi_local); - template void divide_psi>(std::complex* psi, std::complex* psi_local); +} +template void divide_psi(double* psi, double* psi_local); +template void divide_psi>(std::complex* psi, std::complex* psi_local); +template void divide_psi>(std::complex* psi, std::complex* psi_local); #ifdef __MPI - void cal_division(int &npw) +void cal_division(int& npw) +{ + int nprocs, mypnum; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); + for (int i = 0; i < nprocs; i++) { - int nprocs, mypnum; - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); - for(int i=0;i - void divide_hpsi(psi::Psi& psi, psi::Psi& psi_local, std::vector& hmatrix, std::vector& hmatrix_local) - { - int nprocs, mypnum; - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); +template +void divide_hpsi(psi::Psi& psi, psi::Psi& psi_local, std::vector& hmatrix, std::vector& hmatrix_local) +{ + int nprocs, mypnum; + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); - int npw = psi.get_nbasis(); - int nbands = psi.get_nbands(); - int nk = psi.get_nk(); + int npw = psi.get_nbasis(); + int nbands = psi.get_nbands(); + int nk = psi.get_nk(); - hmatrix_local.resize(npw * npw_local[mypnum]); - h_nc = npw_local[mypnum]; - h_nr = npw; - psi_local.resize(nk,nbands,npw_local[mypnum]); - for(int i=0;i(&(hmatrix[i * npw]), &(hmatrix_local[i * npw_local[mypnum]])); + if (i < nbands) { - divide_psi(&(hmatrix[i * npw]), &(hmatrix_local[i * npw_local[mypnum]])); - if(i(psi.get_pointer(i), psi_local.get_pointer(i)); - } + psi.fix_k(k); + psi_local.fix_k(k); + divide_psi(psi.get_pointer(i), psi_local.get_pointer(i)); } } } -#endif } +#endif +} // namespace DIAGOTEST #include "source_base/macros.h" -template +template class HPsi { /** - * This calss used to produce the Hermite matrix, the initial - * guess wave function, and the precondition by the random + * This calss used to produce the Hermite matrix, the initial + * guess wave function, and the precondition by the random * number. The elements of Hermite matrix and wave function are * between -1.0 to 1.0, and the preconddition is between 1.0 to 2.0. - * + * * The parameters in construct function or function create() * are same: * - int nband/nbd: number of calculated bands * - int npw: number of plane wave - * - int sparsity: the sparsity of Halmit matrix, between 0 and 10. + * - int sparsity: the sparsity of Halmit matrix, between 0 and 10. * (0 means no sparsity, 10 means a diagonal matrix) - * + * * After instantiation a HPsi, one can use below functions: * - hamilt(): return the Hermite matrix (type: std::vector) * - precond(): return the precondition (type: Real Pointer) - * + * */ using Real = typename GetTypeReal::type; - public: - HPsi(int nband,int npw, int sparsity=7): - nband(nband),npw(npw),sparsity(sparsity) {genhmatrix();genpsi();genprecondition();} - HPsi(){}; - ~HPsi() {delete [] precondition;} - void create(int nbd,int npw, int sparsity=7) + public: + HPsi(int nband, int npw, int sparsity = 7) : nband(nband), npw(npw), sparsity(sparsity) + { + genhmatrix(); + genpsi(); + genprecondition(); + } + HPsi() {}; + ~HPsi() + { + delete[] precondition; + } + + void create(int nbd, int npw, int sparsity = 7) { this->nband = nbd; this->npw = npw; this->sparsity = sparsity; - genhmatrix();genpsi();genprecondition(); + genhmatrix(); + genpsi(); + genprecondition(); } - //return the matrix - Real* precond() { return precondition; } - std::vector hamilt() { return hmatrix; } - //ModuleBase::ComplexMatrix psi() {return psimatrix;} + // return the matrix + Real* precond() + { + return precondition; + } + std::vector hamilt() + { + return hmatrix; + } + // ModuleBase::ComplexMatrix psi() {return psimatrix;} psi::Psi psi() { psi::Psi psitmp(1, nband, npw, npw, true); - for(int i=0;i u(min, max); - for (int i = 0;i < npw;i++) + for (int i = 0; i < npw; i++) { precondition[i] = 1.0 + static_cast(u(e)) / max; } } -private: + private: int npw; int nband; std::vector hmatrix; @@ -255,7 +303,7 @@ class HPsi int sparsity; }; -template<> +template <> void HPsi::genhmatrix() { hmatrix.resize(npw * npw); @@ -263,27 +311,36 @@ void HPsi::genhmatrix() DIAGOTEST::h_nc = npw; std::default_random_engine e(100); std::uniform_int_distribution u(min, max); - if (sparsity < 0) { sparsity = 0; -} - if (sparsity > 10) { sparsity = 10; -} - for (int i = 0;i < npw;i++) + if (sparsity < 0) { - for (int j = 0;j <= i;j++) + sparsity = 0; + } + if (sparsity > 10) + { + sparsity = 10; + } + for (int i = 0; i < npw; i++) + { + for (int j = 0; j <= i; j++) { double mincoef = 0.0; double realp = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; - if (int(u(e) % 10) > int(sparsity - 1)) { mincoef = 1.0; -} - if (i == j) { + if (int(u(e) % 10) > int(sparsity - 1)) + { + mincoef = 1.0; + } + if (i == j) + { hmatrix[i * npw + j] = realp; - } else { + } + else + { hmatrix[j * npw + i] = hmatrix[i * npw + j] = mincoef * realp; -} + } } } } -template<> +template <> void HPsi>::genhmatrix() { hmatrix.resize(npw * npw); @@ -291,33 +348,39 @@ void HPsi>::genhmatrix() DIAGOTEST::h_nc = npw; std::default_random_engine e(100); std::uniform_int_distribution u(min, max); - if (sparsity < 0) { sparsity = 0; -} - if (sparsity > 10) { sparsity = 10; -} - for (int i = 0;i < npw;i++) + if (sparsity < 0) + { + sparsity = 0; + } + if (sparsity > 10) + { + sparsity = 10; + } + for (int i = 0; i < npw; i++) { - for (int j = 0;j <= i;j++) + for (int j = 0; j <= i; j++) { double mincoef = 0.0; double realp = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; - //double imagp= pow(-1.0,u(e)%2) * static_cast(u(e))/max; - if (int(u(e) % 10) > int(sparsity - 1)) { mincoef = 1.0; -} + // double imagp= pow(-1.0,u(e)%2) * static_cast(u(e))/max; + if (int(u(e) % 10) > int(sparsity - 1)) + { + mincoef = 1.0; + } if (i == j) { - hmatrix[i * npw + j] = std::complex{ realp,0.0 }; + hmatrix[i * npw + j] = std::complex{realp, 0.0}; } else { - //hmatrix(i,j) = mincoef*std::complex{realp,imagp}; - hmatrix[i * npw + j] = mincoef * std::complex{ realp, 0.0 }; + // hmatrix(i,j) = mincoef*std::complex{realp,imagp}; + hmatrix[i * npw + j] = mincoef * std::complex{realp, 0.0}; hmatrix[j * npw + i] = conj(hmatrix[i * npw + j]); } } } } -template<> +template <> void HPsi>::genhmatrix() { hmatrix.resize(npw * npw); @@ -325,78 +388,86 @@ void HPsi>::genhmatrix() DIAGOTEST::h_nc = npw; std::default_random_engine e(100); std::uniform_int_distribution u(min, max); - if (sparsity < 0) { sparsity = 0; -} - if (sparsity > 10) { sparsity = 10; -} - for (int i = 0;i < npw;i++) + if (sparsity < 0) + { + sparsity = 0; + } + if (sparsity > 10) { - for (int j = 0;j <= i;j++) + sparsity = 10; + } + for (int i = 0; i < npw; i++) + { + for (int j = 0; j <= i; j++) { float mincoef = 0.0; float realp = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; - if (int(u(e) % 10) > int(sparsity - 1)) { mincoef = 1.0; -} + if (int(u(e) % 10) > int(sparsity - 1)) + { + mincoef = 1.0; + } if (i == j) { - hmatrix[i * npw + j] = std::complex{ realp,0.0 }; + hmatrix[i * npw + j] = std::complex{realp, 0.0}; } else { - hmatrix[i * npw + j] = mincoef * std::complex{ realp, 0.0 }; + hmatrix[i * npw + j] = mincoef * std::complex{realp, 0.0}; hmatrix[j * npw + i] = conj(hmatrix[i * npw + j]); } } } } -template<> +template <> void HPsi::genpsi() { { psimatrix.resize(nband * npw); std::default_random_engine e(10); std::uniform_int_distribution u(min, max); - for (int i = 0;i < nband;i++) { - for (int j = 0;j < npw;j++) { + for (int i = 0; i < nband; i++) + { + for (int j = 0; j < npw; j++) + { psimatrix[i * npw + j] = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; -} -} + } + } } } -template<> +template <> void HPsi>::genpsi() { { psimatrix.resize(nband * npw); std::default_random_engine e(10); std::uniform_int_distribution u(min, max); - for (int i = 0;i < nband;i++) + for (int i = 0; i < nband; i++) { - for (int j = 0;j < npw;j++) + for (int j = 0; j < npw; j++) { double realp = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; - //double imagp=pow(-1.0,u(e)%2) * static_cast(u(e))/max; + // double imagp=pow(-1.0,u(e)%2) * static_cast(u(e))/max; double imagp = 0.0; - psimatrix[i * npw + j] = std::complex{ realp,imagp }; + psimatrix[i * npw + j] = std::complex{realp, imagp}; } } } } -template<> +template <> void HPsi>::genpsi() { { psimatrix.resize(nband * npw); std::default_random_engine e(10); std::uniform_int_distribution u(min, max); - for (int i = 0;i < nband;i++) + for (int i = 0; i < nband; i++) { - for (int j = 0;j < npw;j++) + for (int j = 0; j < npw; j++) { float realp = pow(-1.0, u(e) % 2) * static_cast(u(e)) / max; float imagp = 0.0; - psimatrix[i * npw + j] = std::complex{ realp,imagp }; + psimatrix[i * npw + j] = std::complex{realp, imagp}; } } } @@ -408,12 +479,25 @@ template class HPsi>; namespace DIAGOTEST { - /// the process-local slice of the test matrix for each scalar type - template std::vector& hmatrix_local_of(); - template <> std::vector& hmatrix_local_of() { return hmatrix_local_d; } - template <> std::vector>& hmatrix_local_of>() { return hmatrix_local; } - template <> std::vector>& hmatrix_local_of>() { return hmatrix_local_f; } +/// the process-local slice of the test matrix for each scalar type +template +std::vector& hmatrix_local_of(); +template <> +std::vector& hmatrix_local_of() +{ + return hmatrix_local_d; +} +template <> +std::vector>& hmatrix_local_of>() +{ + return hmatrix_local; +} +template <> +std::vector>& hmatrix_local_of>() +{ + return hmatrix_local_f; } +} // namespace DIAGOTEST /** * The test matrix seen through hsolver::HSOperator, which is all the iterative diff --git a/source/source_hsolver/test/diago_pexsi_test.cpp b/source/source_hsolver/test/diago_pexsi_test.cpp index 14b23d05b0d..ffe96a1be0c 100644 --- a/source/source_hsolver/test/diago_pexsi_test.cpp +++ b/source/source_hsolver/test/diago_pexsi_test.cpp @@ -1,7 +1,7 @@ #ifdef __PEXSI -#include "source_base/matrix_block.h" #include "source_hsolver/diago_pexsi.h" +#include "source_base/matrix_block.h" #include "source_base/module_external/scalapack_connector.h" #include "source_base/parallel_global.h" #include "source_basis/module_ao/parallel_orbitals.h" @@ -382,8 +382,7 @@ class PexsiGammaOnlyTest : public ::testing::TestWithParam TEST_P(PexsiGammaOnlyTest, LCAO) { const auto& test_case = GetParam(); - PexsiPrepare dp( - 0, 0, test_case.nb2d, 0, test_case.hfname, test_case.sfname, test_case.dmname); + PexsiPrepare dp(0, 0, test_case.nb2d, 0, test_case.hfname, test_case.sfname, test_case.dmname); std::stringstream out_info; if (DETAILINFO && dp.myrank == 0) @@ -417,14 +416,8 @@ INSTANTIATE_TEST_SUITE_P( DiagoTest, PexsiGammaOnlyTest, ::testing::Values( - PexsiTestCase{2, - "PEXSI-H-GammaOnly-Si2.dat", - "PEXSI-S-GammaOnly-Si2.dat", - "PEXSI-DM-GammaOnly-Si2.dat"}, - PexsiTestCase{1, - "PEXSI-H-GammaOnly-Si2.dat", - "PEXSI-S-GammaOnly-Si2.dat", - "PEXSI-DM-GammaOnly-Si2.dat"})); + PexsiTestCase{2, "PEXSI-H-GammaOnly-Si2.dat", "PEXSI-S-GammaOnly-Si2.dat", "PEXSI-DM-GammaOnly-Si2.dat"}, + PexsiTestCase{1, "PEXSI-H-GammaOnly-Si2.dat", "PEXSI-S-GammaOnly-Si2.dat", "PEXSI-DM-GammaOnly-Si2.dat"})); int main(int argc, char** argv) { diff --git a/source/source_hsolver/test/hsolver_pw_sup.h b/source/source_hsolver/test/hsolver_pw_sup.h index 75ec3a90500..e733ea3dc30 100644 --- a/source/source_hsolver/test/hsolver_pw_sup.h +++ b/source/source_hsolver/test/hsolver_pw_sup.h @@ -1,30 +1,34 @@ #include "source_basis/module_pw/pw_basis_k.h" -namespace ModulePW { - -PW_Basis::PW_Basis(){}; -PW_Basis::~PW_Basis(){}; -void PW_Basis::initgrids( - const double lat0_in, // unit length (unit in bohr) - const ModuleBase::Matrix3 - latvec_in, // Unitcell lattice vectors (unit in lat0) - const double gridecut // unit in Ry, ecut to set up grids -) { +namespace ModulePW +{ + +PW_Basis::PW_Basis() {}; +PW_Basis::~PW_Basis() {}; +void PW_Basis::initgrids(const double lat0_in, // unit length (unit in bohr) + const ModuleBase::Matrix3 latvec_in, // Unitcell lattice vectors (unit in lat0) + const double gridecut // unit in Ry, ecut to set up grids +) +{ return; } -void PW_Basis::initgrids( - const double lat0_in, - const ModuleBase::Matrix3 latvec_in, // Unitcell lattice vectors - const int nx_in, - int ny_in, - int nz_in) { +void PW_Basis::initgrids(const double lat0_in, + const ModuleBase::Matrix3 latvec_in, // Unitcell lattice vectors + const int nx_in, + int ny_in, + int nz_in) +{ return; } -void PW_Basis::distribute_r() { return; } +void PW_Basis::distribute_r() +{ + return; +} -PW_Basis_K::PW_Basis_K() { +PW_Basis_K::PW_Basis_K() +{ this->nks = 1; this->npwk_max = 3; // used for update_precondition @@ -34,30 +38,30 @@ PW_Basis_K::PW_Basis_K() { this->tpiba2 = 1.0; } -PW_Basis_K::~PW_Basis_K() { +PW_Basis_K::~PW_Basis_K() +{ delete[] this->gk2; delete[] this->npwk; } -double& PW_Basis_K::getgk2(const int ik, const int igl) const { +double& PW_Basis_K::getgk2(const int ik, const int igl) const +{ this->gk2[igl] = (ik + igl) * 1.5; return this->gk2[igl]; } - } // namespace ModulePW +#include "source_hsolver/diag_comm_info.h" #include "source_hsolver/diago_cg.h" #include "source_hsolver/diago_david.h" -#include "source_hsolver/diag_comm_info.h" #include "source_hsolver/diago_iter_assist.h" - -namespace hsolver { +namespace hsolver +{ template -DiagoCG::DiagoCG(const std::string& basis_type, - const std::string& calculation) +DiagoCG::DiagoCG(const std::string& basis_type, const std::string& calculation) #ifdef __MPI : diag_comm_(MPI_COMM_SELF, 0, 1) #else @@ -78,7 +82,8 @@ DiagoCG::DiagoCG(const std::string& basis_type, const diag_comm_info& diag_comm, const Real& pw_diag_thr, const int& pw_diag_nmax) - : diag_comm_(diag_comm) { + : diag_comm_(diag_comm) +{ basis_type_ = basis_type; calculation_ = calculation; need_subspace_ = need_subspace; @@ -90,7 +95,8 @@ DiagoCG::DiagoCG(const std::string& basis_type, } template -DiagoCG::~DiagoCG() { +DiagoCG::~DiagoCG() +{ delete this->one_; delete this->zero_; delete this->neg_one_; @@ -104,12 +110,15 @@ double DiagoCG::diag(const HSOperator& op, T* psi_in, Real* eigenvalue_in, const std::vector& ethr_band, - const Real* prec) { + const Real* prec) +{ // do something - for (int ib = 0; ib < nband; ib++) { + for (int ib = 0; ib < nband; ib++) + { eigenvalue_in[ib] = 0.0; T* psi_band = psi_in + static_cast(ib) * static_cast(ld_psi); - for (int ig = 0; ig < ld_psi; ig++) { + for (int ig = 0; ig < ld_psi; ig++) + { psi_band[ig] += T(2.0, 0.0); eigenvalue_in[ib] += psi_band[ig].real(); } @@ -128,7 +137,9 @@ DiagoDavid::DiagoDavid(const Real* precondition_in, const int dim_in, const int david_ndim_in, const diag_comm_info& diag_comm_in) - : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), diag_comm(diag_comm_in) { + : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), + diag_comm(diag_comm_in) +{ this->device = base_device::get_device_type(this->ctx); this->precondition = precondition_in; @@ -140,7 +151,8 @@ DiagoDavid::DiagoDavid(const Real* precondition_in, } template -DiagoDavid::~DiagoDavid() { +DiagoDavid::~DiagoDavid() +{ delmem_complex_op()(this->hpsi); delmem_complex_op()(this->spsi); delmem_complex_op()(this->hcc); @@ -157,7 +169,8 @@ int DiagoDavid::diag(const HSOperator& op, const std::vector& ethr_band, const int david_maxiter, const int ntry_max, - const int notconv_max) { + const int notconv_max) +{ // do nothing, we dont need it // do something // for (int ib = 0; ib < psi.get_nbands(); ib++) { diff --git a/source/source_hsolver/test/hsolver_supplementary_mock.h b/source/source_hsolver/test/hsolver_supplementary_mock.h index 76da6b23eba..55107b8bd51 100644 --- a/source/source_hsolver/test/hsolver_supplementary_mock.h +++ b/source/source_hsolver/test/hsolver_supplementary_mock.h @@ -13,14 +13,11 @@ const double* ElecState::getRho(int spin) const return &(this->charge->rho[spin][0]); } - void ElecState::init_nelec_spin() { return; } - - void ElecState::init_ks(Charge* chg_in, // pointer for class Charge const K_Vectors* klist_in, int nk_in, diff --git a/source/source_hsolver/test/parallel_k2d_test.cpp b/source/source_hsolver/test/parallel_k2d_test.cpp index fa5f6479583..e788b439d94 100644 --- a/source/source_hsolver/test/parallel_k2d_test.cpp +++ b/source/source_hsolver/test/parallel_k2d_test.cpp @@ -14,15 +14,23 @@ * Test fixture for class Parallel_K2D */ -class MPIContext { +class MPIContext +{ public: - MPIContext() { + MPIContext() + { MPI_Comm_rank(MPI_COMM_WORLD, &_rank); MPI_Comm_size(MPI_COMM_WORLD, &_size); } - int GetRank() const { return _rank; } - int GetSize() const { return _size; } + int GetRank() const + { + return _rank; + } + int GetSize() const + { + return _size; + } int KPAR; int NPROC_IN_POOL; @@ -34,33 +42,41 @@ class MPIContext { int _size; }; -class ParaPrepare { +class ParaPrepare +{ public: - ParaPrepare(int KPAR_in, int nkstot_in) - : KPAR_(KPAR_in), nkstot_(nkstot_in) {} + ParaPrepare(int KPAR_in, int nkstot_in) : KPAR_(KPAR_in), nkstot_(nkstot_in) + { + } int KPAR_; int nkstot_; }; -class ParallelK2DTest : public ::testing::TestWithParam { +class ParallelK2DTest : public ::testing::TestWithParam +{ protected: Parallel_K2D k2d; MPIContext mpi; int NPROC; int MY_RANK; - void SetUp() override { + void SetUp() override + { NPROC = mpi.GetSize(); MY_RANK = mpi.GetRank(); } - void TearDown() override {} + void TearDown() override + { + } }; -TEST_P(ParallelK2DTest, DividePools) { +TEST_P(ParallelK2DTest, DividePools) +{ ParaPrepare pp = GetParam(); mpi.KPAR = pp.KPAR_; - //k2d.Pkpoints = new Parallel_Kpoints; - //k2d.get_p2D_pool() = new Parallel_2D; - if (mpi.KPAR > NPROC) { + // k2d.Pkpoints = new Parallel_Kpoints; + // k2d.get_p2D_pool() = new Parallel_2D; + if (mpi.KPAR > NPROC) + { std::string output; testing::internal::CaptureStdout(); EXPECT_EXIT(Parallel_Global::divide_mpi_groups(this->NPROC, @@ -72,10 +88,10 @@ TEST_P(ParallelK2DTest, DividePools) { testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); - EXPECT_THAT( - output, - testing::HasSubstr("must be greater than the number of groups")); - } else { + EXPECT_THAT(output, testing::HasSubstr("must be greater than the number of groups")); + } + else + { k2d.set_kpar(mpi.KPAR); k2d.set_para_env(pp.nkstot_, 10, 1, this->NPROC, this->MY_RANK, 1); /* @@ -111,7 +127,7 @@ TEST_P(ParallelK2DTest, DividePools) { k2d.Pkpoints->get_startpro_pool(ipool) << std::endl; } */ - //k2d.get_p2D_pool()->init(10, 10, 1, POOL_WORLD, 0); + // k2d.get_p2D_pool()->init(10, 10, 1, POOL_WORLD, 0); EXPECT_EQ(this->NPROC, 8); EXPECT_EQ(k2d.get_p2D_pool()->dim0, 2); EXPECT_EQ(k2d.get_p2D_pool()->dim1, 2); @@ -121,8 +137,8 @@ TEST_P(ParallelK2DTest, DividePools) { EXPECT_EQ(k2d.get_kpar(), 10); k2d.unset_para_env(); } - //delete k2d.Pkpoints; - //delete k2d.get_p2D_pool(); + // delete k2d.Pkpoints; + // delete k2d.get_p2D_pool(); } INSTANTIATE_TEST_SUITE_P(TESTPK, @@ -131,7 +147,8 @@ INSTANTIATE_TEST_SUITE_P(TESTPK, // KPAR, nkstot ParaPrepare(2, 16))); -int main(int argc, char** argv) { +int main(int argc, char** argv) +{ MPI_Init(&argc, &argv); testing::InitGoogleTest(&argc, argv); int result = RUN_ALL_TESTS(); diff --git a/source/source_hsolver/test/test_diago_david.cpp b/source/source_hsolver/test/test_diago_david.cpp index b2556513944..3fda440fce8 100644 --- a/source/source_hsolver/test/test_diago_david.cpp +++ b/source/source_hsolver/test/test_diago_david.cpp @@ -1,10 +1,10 @@ -#include "source_hsolver/diago_david.h" -#include "source_hsolver/diago_dav_subspace.h" -#include "source_hsolver/diago_cg.h" #include "source_base/parallel_cell.h" #include "source_base/parallel_comm.h" #include "source_base/parallel_global.h" #include "source_base/parallel_reduce.h" +#include "source_hsolver/diago_cg.h" +#include "source_hsolver/diago_dav_subspace.h" +#include "source_hsolver/diago_david.h" #include #include @@ -20,8 +20,12 @@ template class DiagonalOperator final : public hsolver::HSOperator { public: - DiagonalOperator(const int offset, const int size) : offset_(offset), size_(size) {} - void update_k(const int ik) override {} + DiagonalOperator(const int offset, const int size) : offset_(offset), size_(size) + { + } + void update_k(const int ik) override + { + } void hpsi(const T* x, T* hx, const int ld, const int nvec) const override { EXPECT_EQ(ld, size_); @@ -61,7 +65,9 @@ std::complex initial_value(const int row, const int band) } template -class DavidsonOwnershipTest : public ::testing::Test {}; +class DavidsonOwnershipTest : public ::testing::Test +{ +}; using ScalarTypes = ::testing::Types>; TYPED_TEST_SUITE(DavidsonOwnershipTest, ScalarTypes); diff --git a/source/source_hsolver/test/test_diago_hs_para.cpp b/source/source_hsolver/test/test_diago_hs_para.cpp index a32d329fbdf..69905c05319 100644 --- a/source/source_hsolver/test/test_diago_hs_para.cpp +++ b/source/source_hsolver/test/test_diago_hs_para.cpp @@ -1,24 +1,31 @@ #ifndef TEST_DIAGO_PXXXGVX_H #define TEST_DIAGO_PXXXGVX_H -#include -#include -#include -#include -#include -#include -#include - #include "../diag_hs_para.h" #include "source_hsolver/kernels/hegvd_op.h" +#include +#include +#include +#include +#include +#include +#include + template -typename std::enable_if::value || std::is_same::value>::type -generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution::type>& dis, std::vector& h_mat, std::vector& s_mat) { +typename std::enable_if::value || std::is_same::value>::type generate_random_hs_impl( + int d, + std::mt19937& gen, + std::uniform_real_distribution::type>& dis, + std::vector& h_mat, + std::vector& s_mat) +{ // For S matrix, firstly we generate a random symmetric matrix s_tmp, then we set S = s_tmp * s_tmp^T + n * I - std::vector s_tmp(d*d,0); - for (int i = 0; i < d; ++i) { - for (int j = i; j < d; ++j) { + std::vector s_tmp(d * d, 0); + for (int i = 0; i < d; ++i) + { + for (int j = i; j < d; ++j) + { typename GetTypeReal::type value1 = static_cast::type>(dis(gen)); h_mat[i * d + j] = value1; h_mat[j * d + i] = value1; @@ -30,13 +37,17 @@ generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution } // set S = s_tmp * s_tmp^T + n * I - for (int i = 0; i < d; ++i) { - for (int j = 0; j < d; ++j) { + for (int i = 0; i < d; ++i) + { + for (int j = 0; j < d; ++j) + { s_mat[i * d + j] = 0; - for (int k = 0; k < d; ++k) { + for (int k = 0; k < d; ++k) + { s_mat[i * d + j] += s_tmp[i * d + k] * s_tmp[j * d + k]; } - if (i == j) { + if (i == j) + { s_mat[i * d + j] += 2.0; } } @@ -44,11 +55,19 @@ generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution } template -typename std::enable_if>::value || std::is_same>::value>::type -generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution::type>& dis, std::vector& h_mat, std::vector& s_mat) { - std::vector s_tmp(d*d,0); - for (int i = 0; i < d; ++i) { - for (int j = i; j < d; ++j) { +typename std::enable_if>::value + || std::is_same>::value>::type + generate_random_hs_impl(int d, + std::mt19937& gen, + std::uniform_real_distribution::type>& dis, + std::vector& h_mat, + std::vector& s_mat) +{ + std::vector s_tmp(d * d, 0); + for (int i = 0; i < d; ++i) + { + for (int j = i; j < d; ++j) + { typename GetTypeReal::type value1 = static_cast::type>(dis(gen)); typename GetTypeReal::type value2 = static_cast::type>(dis(gen)); h_mat[i * d + j] = T(value1, value2); @@ -56,7 +75,8 @@ generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution { h_mat[j * d + i] = T(value1, -value2); } - else{ + else + { h_mat[j * d + i] = T(value1, 0); } @@ -68,13 +88,17 @@ generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution } // set S = s_tmp * s_tmp^T + n * I - for (int i = 0; i < d; ++i) { - for (int j = 0; j < d; ++j) { + for (int i = 0; i < d; ++i) + { + for (int j = 0; j < d; ++j) + { s_mat[i * d + j] = T(0, 0); - for (int k = 0; k < d; ++k) { + for (int k = 0; k < d; ++k) + { s_mat[i * d + j] += s_tmp[i * d + k] * std::conj(s_tmp[j * d + k]); } - if (i == j) { + if (i == j) + { s_mat[i * d + j] += T(2.0, 0); } } @@ -82,31 +106,48 @@ generate_random_hs_impl(int d, std::mt19937& gen, std::uniform_real_distribution } template -void generate_random_hs(int d, int random_seed ,std::vector& h_mat, std::vector& s_mat) { +void generate_random_hs(int d, int random_seed, std::vector& h_mat, std::vector& s_mat) +{ std::mt19937 gen(random_seed); - std::uniform_real_distribution::type> dis(-1.0,1.0); + std::uniform_real_distribution::type> dis(-1.0, 1.0); h_mat.resize(d * d); s_mat.resize(d * d); generate_random_hs_impl(d, gen, dis, h_mat, s_mat); } - template -typename std::enable_if::value || std::is_same::value>::type -verify_results(const std::vector& h_psi, const std::vector& s_psi, const std::vector::type>& ekb, int lda, int nbands, double threshold) { - for (int i = 0; i < lda; ++i) { - for (int j = 0; j < nbands; ++j) { +typename std::enable_if::value || std::is_same::value>::type verify_results( + const std::vector& h_psi, + const std::vector& s_psi, + const std::vector::type>& ekb, + int lda, + int nbands, + double threshold) +{ + for (int i = 0; i < lda; ++i) + { + for (int j = 0; j < nbands; ++j) + { ASSERT_NEAR(h_psi[j * lda + i], ekb[j] * s_psi[j * lda + i], threshold); } } } template -typename std::enable_if>::value || std::is_same>::value>::type -verify_results(const std::vector& h_psi, const std::vector& s_psi, const std::vector::type>& ekb, int lda, int nbands, double threshold) { - for (int i = 0; i < lda; ++i) { - for (int j = 0; j < nbands; ++j) { +typename std::enable_if>::value + || std::is_same>::value>::type + verify_results(const std::vector& h_psi, + const std::vector& s_psi, + const std::vector::type>& ekb, + int lda, + int nbands, + double threshold) +{ + for (int i = 0; i < lda; ++i) + { + for (int j = 0; j < nbands; ++j) + { ASSERT_NEAR(h_psi[j * lda + i].real(), ekb[j] * s_psi[j * lda + i].real(), threshold); ASSERT_NEAR(h_psi[j * lda + i].imag(), ekb[j] * s_psi[j * lda + i].imag(), threshold); } @@ -114,35 +155,41 @@ verify_results(const std::vector& h_psi, const std::vector& s_psi, const s } template -void test_diago_hs(int lda, int nb, int random_seed, int nbands, int diag_type, MPI_Comm comm) { +void test_diago_hs(int lda, int nb, int random_seed, int nbands, int diag_type, MPI_Comm comm) +{ // diag_type should be 1 (for elpa) or 2 (for scalapack) int my_rank; MPI_Comm_rank(comm, &my_rank); std::vector h_mat, s_mat, wfc, h_psi, s_psi; std::vector::type> ekb(lda); - if (my_rank==0) + if (my_rank == 0) { h_mat.resize(lda * lda); s_mat.resize(lda * lda); wfc.resize(lda * lda); generate_random_hs(lda, random_seed, h_mat, s_mat); } - hsolver::diago_hs_para(h_mat.data(), s_mat.data(), lda, nbands,ekb.data(), wfc.data(), comm, diag_type, nb); + hsolver::diago_hs_para(h_mat.data(), s_mat.data(), lda, nbands, ekb.data(), wfc.data(), comm, diag_type, nb); // Verify results - if (my_rank == 0){ + if (my_rank == 0) + { double threshold = 1e-6; - if (std::is_same>::value || std::is_same::value) { + if (std::is_same>::value || std::is_same::value) + { threshold = 1e-12; } h_psi.resize(lda * nbands, 0); s_psi.resize(lda * nbands, 0); - for (int i = 0; i < lda; ++i) { - for (int j = 0; j < nbands; ++j) { - for (int k = 0; k < lda; ++k) { + for (int i = 0; i < lda; ++i) + { + for (int j = 0; j < nbands; ++j) + { + for (int k = 0; k < lda; ++k) + { h_psi[j * lda + i] += h_mat[k * lda + i] * wfc[j * lda + k]; s_psi[j * lda + i] += s_mat[k * lda + i] * wfc[j * lda + k]; } @@ -153,7 +200,8 @@ void test_diago_hs(int lda, int nb, int random_seed, int nbands, int diag_type, } template -void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, int loop_numb) { +void test_performance(int lda, int nb, int nbands, MPI_Comm comm, int case_numb, int loop_numb) +{ // generate 10 random H/S, and do the diagonalization 100 times by using elpa/scalapack and lapack. int my_rank, nproc; MPI_Comm_rank(comm, &my_rank); @@ -166,7 +214,7 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, std::vector::type> ekb_scalap(lda); std::vector::type> ekb_lapack(lda); - if (my_rank==0) + if (my_rank == 0) { std::cout << "\nMatrix size: " << lda << " x " << lda << std::endl; std::cout << "Number of bands: " << nbands << std::endl; @@ -184,12 +232,15 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, std::vector time_scalap(case_numb, 0); std::vector time_lapack(case_numb, 0); - if (my_rank == 0) { std::cout << "Random matrix "; -} + if (my_rank == 0) + { + std::cout << "Random matrix "; + } for (int randomi = 0; randomi < case_numb; ++randomi) { - if (my_rank == 0) { + if (my_rank == 0) + { std::cout << randomi << " "; generate_random_hs(lda, randomi, h_mat, s_mat); } @@ -199,9 +250,17 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, // ELPA MPI_Barrier(comm); start = std::chrono::high_resolution_clock::now(); - for (int j=0;j(h_mat.data(), s_mat.data(), lda, nbands,ekb_elpa.data(), wfc.data(), comm, 1, nb); + hsolver::diago_hs_para(h_mat.data(), + s_mat.data(), + lda, + nbands, + ekb_elpa.data(), + wfc.data(), + comm, + 1, + nb); MPI_Barrier(comm); } MPI_Barrier(comm); @@ -211,40 +270,49 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, // scalapack start = std::chrono::high_resolution_clock::now(); - for (int j=0;j(h_mat.data(), s_mat.data(), lda, nbands,ekb_scalap.data(), wfc.data(), comm, 2, nb); + hsolver::diago_hs_para(h_mat.data(), + s_mat.data(), + lda, + nbands, + ekb_scalap.data(), + wfc.data(), + comm, + 2, + nb); MPI_Barrier(comm); } MPI_Barrier(comm); end = std::chrono::high_resolution_clock::now(); time_scalap[randomi] = std::chrono::duration_cast(end - start).count(); - //LApack + // LApack if (my_rank == 0) { std::vector h_tmp, s_tmp; start = std::chrono::high_resolution_clock::now(); base_device::DEVICE_CPU* ctx = {}; - for (int j=0;j()(ctx, - lda, - lda, - h_tmp.data(), - s_tmp.data(), - nbands, - ekb_lapack.data(), - wfc.data()); + hsolver::hegvx_op()(ctx, + lda, + lda, + h_tmp.data(), + s_tmp.data(), + nbands, + ekb_lapack.data(), + wfc.data()); } end = std::chrono::high_resolution_clock::now(); time_lapack[randomi] = std::chrono::duration_cast(end - start).count(); - //COMPARE EKB - for (int i = 0; i < nbands; ++i) { + // COMPARE EKB + for (int i = 0; i < nbands; ++i) + { typename GetTypeReal::type diff_scalap_lapack = std::abs(ekb_scalap[i] - ekb_lapack[i]); #ifdef __ELPA typename GetTypeReal::type diff_elpa_lapack = std::abs(ekb_elpa[i] - ekb_lapack[i]); @@ -262,26 +330,31 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, } } MPI_Barrier(comm); - } if (my_rank == 0) { #ifdef __ELPA std::cout << "\nELPA Time : "; - for (int i=0; i < case_numb;i++) - {std::cout << time_elpa[i] << " ";} + for (int i = 0; i < case_numb; i++) + { + std::cout << time_elpa[i] << " "; + } std::cout << std::endl; #endif std::cout << "scalapack Time: "; - for (int i=0; i < case_numb;i++) - {std::cout << time_scalap[i] << " ";} + for (int i = 0; i < case_numb; i++) + { + std::cout << time_scalap[i] << " "; + } std::cout << std::endl; std::cout << "lapack Time : "; - for (int i=0; i < case_numb;i++) - {std::cout << time_lapack[i] << " ";} + for (int i = 0; i < case_numb; i++) + { + std::cout << time_lapack[i] << " "; + } std::cout << std::endl; // print out the average time and speedup @@ -290,7 +363,7 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, #endif double avg_time_scalap = 0; double avg_time_lapack = 0; - for (int i=0; i < case_numb;i++) + for (int i = 0; i < case_numb; i++) { #ifdef __ELPA avg_time_elpa += time_elpa[i]; @@ -306,48 +379,57 @@ void test_performance(int lda, int nb, int nbands, MPI_Comm comm,int case_numb, avg_time_lapack /= case_numb; std::cout << "Average Lapack Time : " << avg_time_lapack << " ms" << std::endl; #ifdef __ELPA - std::cout << "Average ELPA Time : " << avg_time_elpa << " ms, Speedup: " << avg_time_lapack / avg_time_elpa << std::endl; + std::cout << "Average ELPA Time : " << avg_time_elpa << " ms, Speedup: " << avg_time_lapack / avg_time_elpa + << std::endl; #endif - std::cout << "Average Scalapack Time: " << avg_time_scalap << " ms, Speedup: " << avg_time_lapack / avg_time_scalap << std::endl; + std::cout << "Average Scalapack Time: " << avg_time_scalap + << " ms, Speedup: " << avg_time_lapack / avg_time_scalap << std::endl; } } -//test_diago_hs(int lda, int nb, int random_seed, int nbands, int diag_type, MPI_Comm comm) -TEST(DiagoPxxxgvxElpaTest, Double) { +// test_diago_hs(int lda, int nb, int random_seed, int nbands, int diag_type, MPI_Comm comm) +TEST(DiagoPxxxgvxElpaTest, Double) +{ test_diago_hs(16, 4, 0, 10, 1, MPI_COMM_WORLD); test_diago_hs(20, 6, 0, 18, 1, MPI_COMM_WORLD); } -TEST(DiagoPxxxgvxElpaTest, ComplexDouble) { +TEST(DiagoPxxxgvxElpaTest, ComplexDouble) +{ test_diago_hs>(16, 6, 0, 10, 1, MPI_COMM_WORLD); test_diago_hs>(20, 8, 0, 18, 1, MPI_COMM_WORLD); } -TEST(DiagoPxxxgvxScalapackTest, Double) { - test_diago_hs(16, 4, 0, 10, 2,MPI_COMM_WORLD); - test_diago_hs(20, 6, 0, 18, 2,MPI_COMM_WORLD); +TEST(DiagoPxxxgvxScalapackTest, Double) +{ + test_diago_hs(16, 4, 0, 10, 2, MPI_COMM_WORLD); + test_diago_hs(20, 6, 0, 18, 2, MPI_COMM_WORLD); } -TEST(DiagoPxxxgvxScalapackTest, ComplexDouble) { +TEST(DiagoPxxxgvxScalapackTest, ComplexDouble) +{ test_diago_hs>(16, 4, 0, 10, 2, MPI_COMM_WORLD); } -TEST(DiagoPxxxgvxScalapackTest, Float) { - test_diago_hs(16, 4, 0, 10,2,MPI_COMM_WORLD); +TEST(DiagoPxxxgvxScalapackTest, Float) +{ + test_diago_hs(16, 4, 0, 10, 2, MPI_COMM_WORLD); } -TEST(DiagoPxxxgvxScalapackTest, ComplexFloat) { - test_diago_hs>(16, 4, 0, 10,2,MPI_COMM_WORLD); +TEST(DiagoPxxxgvxScalapackTest, ComplexFloat) +{ + test_diago_hs>(16, 4, 0, 10, 2, MPI_COMM_WORLD); } -//TEST(DiagoPxxxgvxPerformanceTest, Double) { -// int ndim = 200; -// int nband = 180; -// int case_numb = 10; -// int loop_numb = 10; -// test_performance>(ndim, 32, nband, MPI_COMM_WORLD, case_numb, loop_numb); -//} +// TEST(DiagoPxxxgvxPerformanceTest, Double) { +// int ndim = 200; +// int nband = 180; +// int case_numb = 10; +// int loop_numb = 10; +// test_performance>(ndim, 32, nband, MPI_COMM_WORLD, case_numb, loop_numb); +// } -int main(int argc, char** argv) { +int main(int argc, char** argv) +{ MPI_Init(&argc, &argv); int myrank; MPI_Comm_rank(MPI_COMM_WORLD, &myrank); @@ -355,7 +437,8 @@ int main(int argc, char** argv) { testing::InitGoogleTest(&argc, argv); ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners(); - if (myrank != 0) { + if (myrank != 0) + { delete listeners.Release(listeners.default_result_printer()); } @@ -365,11 +448,11 @@ int main(int argc, char** argv) { std::cout << "ERROR:some tests are not passed" << std::endl; MPI_Finalize(); return result; - } + } MPI_Finalize(); - return 0; + return 0; } #endif // TEST_DIAGO_PXXXGVX_H \ No newline at end of file diff --git a/source/source_hsolver/test/test_hsolver.cpp b/source/source_hsolver/test/test_hsolver.cpp index 97629b3a002..65cacab0f53 100644 --- a/source/source_hsolver/test/test_hsolver.cpp +++ b/source/source_hsolver/test/test_hsolver.cpp @@ -1,12 +1,11 @@ -#include +#include "source_hsolver/hsolver.h" #include #include +#include #include #include -#include "source_hsolver/hsolver.h" - /************************************************ * unit test of the free functions in hsolver.h ***********************************************/ diff --git a/source/source_hsolver/test/test_hsolver_linear.cpp b/source/source_hsolver/test/test_hsolver_linear.cpp index 14efc6a96d8..c96665e8b7e 100644 --- a/source/source_hsolver/test/test_hsolver_linear.cpp +++ b/source/source_hsolver/test/test_hsolver_linear.cpp @@ -40,7 +40,11 @@ template class DenseOperator final : public hsolver::LinearOperator { public: - DenseOperator(const std::vector& matrix, const int n, const int start, const int dim, const hsolver::diag_comm_info& comm) + DenseOperator(const std::vector& matrix, + const int n, + const int start, + const int dim, + const hsolver::diag_comm_info& comm) : matrix_(matrix), n_(n), start_(start), dim_(dim), comm_(comm) { } @@ -117,7 +121,8 @@ TYPED_TEST(LinearSolveTest, DenseReferenceAndReusedBatches) const hsolver::diag_comm_info comm = linear_test::world_comm(); const double tolerance = std::is_same::value ? 1e-12 : 2e-6; const double solution_tol = std::is_same::value ? 1e-10 : 1e-4; - for (const hsolver::LinearMethod method: {hsolver::LinearMethod::bicgstab, hsolver::LinearMethod::cgs, hsolver::LinearMethod::gmres}) + for (const hsolver::LinearMethod method: + {hsolver::LinearMethod::bicgstab, hsolver::LinearMethod::cgs, hsolver::LinearMethod::gmres}) { for (const bool preconditioned: {false, true}) { @@ -208,7 +213,8 @@ TYPED_TEST(LinearSolveTest, DenseReferenceAndReusedBatches) for (int i = 0; i < dim; ++i) { full[band * n + start + i] = x[band * ld + i]; - EXPECT_LT(std::abs(full[band * n + start + i] - reference[band * n + start + i]), solution_tol); + EXPECT_LT(std::abs(full[band * n + start + i] - reference[band * n + start + i]), + solution_tol); } for (int i = dim; i < ld; ++i) { @@ -261,7 +267,8 @@ TEST(LinearSolveFailure, ReportsOriginalColumnAndTrueResidual) const int start = n * comm.rank / comm.nproc; const int dim = n * (comm.rank + 1) / comm.nproc - start; const int ld = dim + 1; - for (const hsolver::LinearMethod method: {hsolver::LinearMethod::bicgstab, hsolver::LinearMethod::cgs, hsolver::LinearMethod::gmres}) + for (const hsolver::LinearMethod method: + {hsolver::LinearMethod::bicgstab, hsolver::LinearMethod::cgs, hsolver::LinearMethod::gmres}) { for (const bool singular: {false, true}) { @@ -284,7 +291,8 @@ TEST(LinearSolveFailure, ReportsOriginalColumnAndTrueResidual) b[ld + i] = T(1, 0.1 * (start + i)); } const IdentityOperator identity(dim); - const hsolver::LinearSolveResult result = solver.solve(op, identity, ld, 3, dim, x.data(), b.data(), nullptr, true); + const hsolver::LinearSolveResult result + = solver.solve(op, identity, ld, 3, dim, x.data(), b.data(), nullptr, true); const hsolver::LinearSolveStatus expected = singular ? hsolver::LinearSolveStatus::breakdown : hsolver::LinearSolveStatus::max_iterations; EXPECT_EQ(result.status, expected); @@ -351,7 +359,8 @@ void check_reconstruction(const bool corrupt_initial) const IdentityOperator identity(dim); hsolver::HSolverLinear solver(options, comm); const T* initial = corrupt_initial ? bad_residual.data() : nullptr; - const hsolver::LinearSolveResult result = solver.solve(op, identity, ld, bands, dim, x.data(), b.data(), initial, corrupt_initial); + const hsolver::LinearSolveResult result + = solver.solve(op, identity, ld, bands, dim, x.data(), b.data(), initial, corrupt_initial); ASSERT_EQ(result.status, hsolver::LinearSolveStatus::converged); EXPECT_EQ(result.failed_band, -1); EXPECT_GT(result.restarts, 0); @@ -369,7 +378,8 @@ void check_reconstruction(const bool corrupt_initial) { error += std::norm(linear_test::Complex(ax[band * ld + i]) - linear_test::Complex(b[band * ld + i])); norm += std::norm(linear_test::Complex(b[band * ld + i])); - EXPECT_LT(std::abs(linear_test::Complex(x[band * ld + i]) - reference[band * n + start + i]), 10 * options.tolerance); + EXPECT_LT(std::abs(linear_test::Complex(x[band * ld + i]) - reference[band * n + start + i]), + 10 * options.tolerance); } #ifdef __MPI Parallel_Common::reduce_data(&error, 1, comm.comm); @@ -383,7 +393,8 @@ void check_reconstruction(const bool corrupt_initial) options.max_iterations = 1; hsolver::HSolverLinear limited(options, comm); std::fill(x.begin(), x.end(), T(0)); - const hsolver::LinearSolveResult exhausted = limited.solve(op, identity, ld, bands, dim, x.data(), b.data(), initial, true); + const hsolver::LinearSolveResult exhausted + = limited.solve(op, identity, ld, bands, dim, x.data(), b.data(), initial, true); EXPECT_EQ(exhausted.status, hsolver::LinearSolveStatus::max_iterations); EXPECT_EQ(exhausted.failed_band, 1); EXPECT_EQ(exhausted.iterations, options.max_iterations); diff --git a/source/source_hsolver/test/test_hsolver_pw.cpp b/source/source_hsolver/test/test_hsolver_pw.cpp index 694349948d8..e11574875f4 100644 --- a/source/source_hsolver/test/test_hsolver_pw.cpp +++ b/source/source_hsolver/test/test_hsolver_pw.cpp @@ -1,8 +1,3 @@ -#include -#include -#include -#include - #include "hsolver_pw_sup.h" #include "hsolver_supplementary_mock.h" #include "source_hsolver/diag_comm_info.h" @@ -12,6 +7,10 @@ #include #include +#include +#include +#include +#include /// H = S = identity: the simplest operator the solvers can be handed template @@ -32,66 +31,63 @@ class IdentityHSOperator : public hsolver::HSOperator - void PW_Basis_K::recip_to_real(const Device* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const FPTYPE factor) const - { - // Simple mock implementation that does nothing - // In a real test, you might want to implement behavior that simulates the actual function - } - - // Mock implementation for real_to_recip - template - void PW_Basis_K::real_to_recip(const Device* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const FPTYPE factor) const - { - // Simple mock implementation that does nothing - } +namespace ModulePW +{ +// Mock implementation for recip_to_real +template +void PW_Basis_K::recip_to_real(const Device* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const FPTYPE factor) const +{ + // Simple mock implementation that does nothing + // In a real test, you might want to implement behavior that simulates the actual function +} - // Explicit template instantiations - template void PW_Basis_K::recip_to_real( - const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const float factor) const; - - template void PW_Basis_K::recip_to_real( - const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; - - template void PW_Basis_K::real_to_recip( - const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const float factor) const; - - template void PW_Basis_K::real_to_recip( - const base_device::DEVICE_CPU* ctx, - const std::complex* in, - std::complex* out, - const int ik, - const bool add, - const double factor) const; +// Mock implementation for real_to_recip +template +void PW_Basis_K::real_to_recip(const Device* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const FPTYPE factor) const +{ + // Simple mock implementation that does nothing } +// Explicit template instantiations +template void PW_Basis_K::recip_to_real(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const float factor) const; + +template void PW_Basis_K::recip_to_real(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; + +template void PW_Basis_K::real_to_recip(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const float factor) const; + +template void PW_Basis_K::real_to_recip(const base_device::DEVICE_CPU* ctx, + const std::complex* in, + std::complex* out, + const int ik, + const bool add, + const double factor) const; +} // namespace ModulePW + /************************************************ * unit test of HSolverPW class ***********************************************/ @@ -112,7 +108,8 @@ namespace ModulePW { */ // mock diago_hs_para -namespace hsolver { +namespace hsolver +{ template void diago_hs_para(T* h, T* s, @@ -123,7 +120,8 @@ void diago_hs_para(T* h, const MPI_Comm& comm, const int diag_subspace, const int block_size = 0) -{} +{ +} template void diago_hs_para(double* h, double* s, const int lda, @@ -153,7 +151,7 @@ template void diago_hs_para(float* h, const MPI_Comm& comm, const int diag_subspace, const int block_size); - + template void diago_hs_para>(std::complex* h, std::complex* s, const int lda, @@ -164,9 +162,10 @@ template void diago_hs_para>(std::complex* h, const int diag_subspace, const int block_size); -} +} // namespace hsolver -class TestHSolverPW : public ::testing::Test { +class TestHSolverPW : public ::testing::Test +{ public: // HSolverPW declares this fixture a friend, but a TEST_F body lives in a // class derived from it and friendship is not inherited, so the call into @@ -273,7 +272,7 @@ class TestHSolverPW : public ::testing::Test { // 1, // true); - + // // EXPECT_EQ(this->hs_d.initialed_psi, true); // EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, // 0.0); @@ -365,13 +364,14 @@ class TestHSolverPW : public ::testing::Test { // // EXPECT_EQ(test_diagethr_d, 0.0001); // } -TEST_F(TestHSolverPW, SolveLcaoInPW) { +TEST_F(TestHSolverPW, SolveLcaoInPW) +{ pwbk.nks = 1; // initial memory and data elecstate_test.ekb.create(1, 2); - elecstate_test.wg.create(1,2); - elecstate_test.klist=new K_Vectors; - elecstate_test.skip_weights=true; + elecstate_test.wg.create(1, 2); + elecstate_test.klist = new K_Vectors; + elecstate_test.skip_weights = true; elecstate_test.pot = new elecstate::Potential; // 1 kpt, 2 bands, 3 basis psi_test_cf.resize(1, 2, 3); @@ -387,15 +387,14 @@ TEST_F(TestHSolverPW, SolveLcaoInPW) { // 1, 2, 3 / 4, 5, 6 / 7, 8, 9 would be rank deficient, so the diagonal is // lifted to keep the three subspace vectors linearly independent - for (int iband = 0; iband < transform_test_cd.get_nbands(); iband++) { - for (int ibasis = 0; ibasis < transform_test_cd.get_nbasis(); - ibasis++) { + for (int iband = 0; iband < transform_test_cd.get_nbands(); iband++) + { + for (int ibasis = 0; ibasis < transform_test_cd.get_nbasis(); ibasis++) + { const double value = iband * transform_test_cd.get_nbasis() + ibasis + 1 + (iband == ibasis ? 10.0 : 0.0); - transform_test_cd - .get_pointer()[iband * transform_test_cd.get_nbasis() + ibasis] + transform_test_cd.get_pointer()[iband * transform_test_cd.get_nbasis() + ibasis] = std::complex(value, 0.0); - transform_test_cf - .get_pointer()[iband * transform_test_cf.get_nbasis() + ibasis] + transform_test_cf.get_pointer()[iband * transform_test_cf.get_nbasis() + ibasis] = std::complex(value, 0.0); } } @@ -404,10 +403,13 @@ TEST_F(TestHSolverPW, SolveLcaoInPW) { auto check_orthonormal = [](const auto& p, const double tol) { const int nb = p.get_nbands(); const int nbasis = p.get_nbasis(); - for (int i = 0; i < nb; i++) { - for (int j = 0; j < nb; j++) { + for (int i = 0; i < nb; i++) + { + for (int j = 0; j < nb; j++) + { std::complex dot = 0.0; - for (int ig = 0; ig < nbasis; ig++) { + for (int ig = 0; ig < nbasis; ig++) + { dot += std::conj(std::complex(p.get_pointer()[i * nbasis + ig])) * std::complex(p.get_pointer()[j * nbasis + ig]); } @@ -430,13 +432,7 @@ TEST_F(TestHSolverPW, SolveLcaoInPW) { const hsolver::diag_comm_info diag_comm(0, 1); #endif std::ostringstream log; - hs_f_lip.solve(hamilt_test_f, - psi_test_cf, - &elecstate_test, - transform_test_cf, - diag_comm, - log, - true); + hs_f_lip.solve(hamilt_test_f, psi_test_cf, &elecstate_test, transform_test_cf, diag_comm, log, true); EXPECT_NE(log.str().find("Average iterative diagonalization steps"), std::string::npos); EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, 0.0); check_orthonormal(psi_test_cf, 1e-5); @@ -445,13 +441,7 @@ TEST_F(TestHSolverPW, SolveLcaoInPW) { elecstate_test.ekb.c[0] = 1.0; elecstate_test.ekb.c[1] = 2.0; - hs_d_lip.solve(hamilt_test_d, - psi_test_cd, - &elecstate_test, - transform_test_cd, - diag_comm, - log, - true); + hs_d_lip.solve(hamilt_test_d, psi_test_cd, &elecstate_test, transform_test_cd, diag_comm, log, true); EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, 0.0); check_orthonormal(psi_test_cd, 1e-10); EXPECT_NEAR(elecstate_test.ekb.c[0], 1.0, 1e-10); @@ -467,9 +457,7 @@ TEST_F(TestHSolverPW, NpwxLessThanNbandsDeath) std::vector precond(3, 0.0); std::vector eigenvalues(5, 0.0); // Expect death from WARNING_QUIT due to npwx < nbands - EXPECT_EXIT( - hamiltSolvePsiK(hs_d, hamilt_test_d, psi_test_cd, precond, eigenvalues.data(), 1), - ::testing::ExitedWithCode(1), - ".*" - ); + EXPECT_EXIT(hamiltSolvePsiK(hs_d, hamilt_test_d, psi_test_cd, precond, eigenvalues.data(), 1), + ::testing::ExitedWithCode(1), + ".*"); } diff --git a/source/source_hsolver/test/test_hsolver_pw_tddft.cpp b/source/source_hsolver/test/test_hsolver_pw_tddft.cpp index 8500e142329..5b15d8c5b56 100644 --- a/source/source_hsolver/test/test_hsolver_pw_tddft.cpp +++ b/source/source_hsolver/test/test_hsolver_pw_tddft.cpp @@ -84,7 +84,11 @@ void initialize_basis(const std::vector& sizes, ModulePW::PW_Basis_K* basis } template -std::vector dense_cn(const std::vector& h, const std::vector& previous, int n, int bands, double dt) +std::vector dense_cn(const std::vector& h, + const std::vector& previous, + int n, + int bands, + double dt) { std::vector lhs(n * n); std::vector rhs_matrix(n * n); @@ -102,7 +106,10 @@ std::vector dense_cn(const std::vector& h, const std::v } template -void check_dense_step(const DenseHamiltonian& op, const psi::Psi& previous, const psi::Psi& current, const double dt) +void check_dense_step(const DenseHamiltonian& op, + const psi::Psi& previous, + const psi::Psi& current, + const double dt) { const double tolerance = std::is_same>::value ? 1e-4 : 1e-10; const int bands = current.get_nbands(); diff --git a/source/source_hsolver/test/test_linear_low_rank.cpp b/source/source_hsolver/test/test_linear_low_rank.cpp index 6a4d110fba0..8ad9b8cb390 100644 --- a/source/source_hsolver/test/test_linear_low_rank.cpp +++ b/source/source_hsolver/test/test_linear_low_rank.cpp @@ -73,7 +73,11 @@ TYPED_TEST(LinearLowRankTest, CoarseAndHistoryCorrectionsMatchDenseReference) hsolver::LinearLowRank preconditioner(algebra, local_d.data(), dim); if (history) { - preconditioner.prepare_response(ld, response.rank(), response.directions(), response.images(), &correction_workspace); + preconditioner.prepare_response(ld, + response.rank(), + response.directions(), + response.images(), + &correction_workspace); } else { diff --git a/source/source_hsolver/test/test_orthonormal.cpp b/source/source_hsolver/test/test_orthonormal.cpp index 50ff4231d75..b6624fcac5e 100644 --- a/source/source_hsolver/test/test_orthonormal.cpp +++ b/source/source_hsolver/test/test_orthonormal.cpp @@ -85,7 +85,8 @@ TYPED_TEST(OrthonormalTest, CorrectsComplexOverlapsWithoutFallback) const hsolver::diag_comm_info comm = self_comm(); hsolver::LinearAlgebra algebra(comm); hsolver::Orthonormal orth(comm, algebra); - for (const hsolver::OrthMethod method: {hsolver::OrthMethod::cholesky, hsolver::OrthMethod::lowdin, hsolver::OrthMethod::newton_schulz}) + for (const hsolver::OrthMethod method: + {hsolver::OrthMethod::cholesky, hsolver::OrthMethod::lowdin, hsolver::OrthMethod::newton_schulz}) { SCOPED_TRACE(hsolver::orth_method_name(method)); std::vector psi = orbitals(); @@ -113,7 +114,8 @@ TYPED_TEST(OrthonormalTest, NewtonSchulzFallsBackOutsideItsRange) std::vector psi(ld * bands, T(0)); psi[0] = T(2); psi[ld + 1] = T(2); - const hsolver::OrthResult result = orth.apply(psi.data(), ld, dim, bands, hsolver::OrthMethod::newton_schulz, false); + const hsolver::OrthResult result + = orth.apply(psi.data(), ld, dim, bands, hsolver::OrthMethod::newton_schulz, false); ASSERT_EQ(result.status, hsolver::OrthStatus::accepted); EXPECT_EQ(result.actual, hsolver::OrthMethod::cholesky); EXPECT_EQ(result.fallbacks, 1); @@ -141,8 +143,9 @@ TYPED_TEST(OrthonormalTest, RejectsInvalidOrbitalsWithoutOverwritingThem) const std::vector original(psi); const hsolver::OrthResult result = orth.apply(psi.data(), ld, dim, bands, method, false); EXPECT_EQ(result.status, hsolver::OrthStatus::failed); - const hsolver::OrthFailure nonfinite - = method == hsolver::OrthMethod::none ? hsolver::OrthFailure::nonfinite_norm : hsolver::OrthFailure::nonfinite_gram; + const hsolver::OrthFailure nonfinite = method == hsolver::OrthMethod::none + ? hsolver::OrthFailure::nonfinite_norm + : hsolver::OrthFailure::nonfinite_gram; const hsolver::OrthFailure expected = value == 0.0 ? hsolver::OrthFailure::nonpositive_norm : nonfinite; EXPECT_EQ(result.failure, expected); EXPECT_EQ(std::memcmp(psi.data(), original.data(), psi.size() * sizeof(T)), 0); @@ -152,7 +155,8 @@ TYPED_TEST(OrthonormalTest, RejectsInvalidOrbitalsWithoutOverwritingThem) duplicate[0] = T(1); duplicate[ld] = T(1); const std::vector original(duplicate); - const hsolver::OrthResult result = orth.apply(duplicate.data(), ld, dim, bands, hsolver::OrthMethod::cholesky, false); + const hsolver::OrthResult result + = orth.apply(duplicate.data(), ld, dim, bands, hsolver::OrthMethod::cholesky, false); EXPECT_EQ(result.status, hsolver::OrthStatus::failed); EXPECT_EQ(result.failure, hsolver::OrthFailure::factorization_failed); EXPECT_EQ(std::memcmp(duplicate.data(), original.data(), duplicate.size() * sizeof(T)), 0); diff --git a/source/source_hsolver/test/test_para_linear_trans.cpp b/source/source_hsolver/test/test_para_linear_trans.cpp index 56dbbbfaf8f..045b6d9a233 100644 --- a/source/source_hsolver/test/test_para_linear_trans.cpp +++ b/source/source_hsolver/test/test_para_linear_trans.cpp @@ -1,5 +1,4 @@ #include "../para_lin_tf.h" - #include "source_base/kernels/math_kernel_op.h" #include