diff --git a/CMakeLists.txt b/CMakeLists.txt index ea1373267b..4162e0fefc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -255,7 +255,7 @@ set_property( "NVHPCBlas" "CUSTOM") -set(CP2K_SCALAPACK_VENDOR_LIST "MKL" "SCI" "GENERIC" "auto") +set(CP2K_SCALAPACK_VENDOR_LIST "MKL" "NVPL" "SCI" "GENERIC" "auto") set(CP2K_SCALAPACK_VENDOR "auto" CACHE STRING "ScaLAPACK vendor/generic backend") @@ -268,6 +268,17 @@ if(DEFINED CP2K_SCALAPACK_VENDOR) endif() endif() +set(CP2K_NVPL_SCALAPACK_MPI + "auto" + CACHE STRING "NVPL BLACS MPI interface") +set(CP2K_NVPL_SCALAPACK_MPI_LIST "auto" "mpich" "openmpi3" "openmpi4" + "openmpi5") +set_property(CACHE CP2K_NVPL_SCALAPACK_MPI + PROPERTY STRINGS ${CP2K_NVPL_SCALAPACK_MPI_LIST}) +if(NOT CP2K_NVPL_SCALAPACK_MPI IN_LIST CP2K_NVPL_SCALAPACK_MPI_LIST) + message(FATAL_ERROR "An invalid NVPL ScaLAPACK MPI interface was specified") +endif() + set(CP2K_DATA_DIR "default" CACHE STRING "Set the location for CP2K data") @@ -285,6 +296,7 @@ set(CP2K_SUPPORTED_CUDA_ARCHITECTURES V100 A100 H100 + GB10 A40) set(CP2K_SUPPORTED_HIP_ARCHITECTURES Mi50 @@ -299,6 +311,7 @@ set(CP2K_SUPPORTED_HIP_ARCHITECTURES V100 A100 H100 + GB10 A40) set(CP2K_WITH_GPU @@ -460,6 +473,7 @@ if((CP2K_USE_ACCEL MATCHES CUDA) OR (CP2K_USE_ACCEL MATCHES HIP)) set(CP2K_GPU_ARCH_NUMBER_V100 70) set(CP2K_GPU_ARCH_NUMBER_A100 80) set(CP2K_GPU_ARCH_NUMBER_H100 90) + set(CP2K_GPU_ARCH_NUMBER_GB10 121) set(CP2K_GPU_ARCH_NUMBER_A40 86) set(CP2K_GPU_ARCH_NUMBER_Mi50 gfx906) set(CP2K_GPU_ARCH_NUMBER_Mi100 gfx908) @@ -700,8 +714,39 @@ endif() if(CP2K_USE_DLAF) find_package(DLAFFortran 0.4.0 REQUIRED) + if(TARGET DLAF::Fortran) + set(CP2K_DLAF_FORTRAN_TARGET DLAF::Fortran) + elseif(TARGET DLAF::DLAF_Fortran) + set(CP2K_DLAF_FORTRAN_TARGET DLAF::DLAF_Fortran) + else() + message(FATAL_ERROR "Could not find a DLA-Future-Fortran CMake target") + endif() + get_target_property(CP2K_DLAF_INCLUDE_DIRS DLAF::DLAF INTERFACE_INCLUDE_DIRECTORIES) + get_target_property(CP2K_DLAF_FORTRAN_INCLUDE_DIRS + ${CP2K_DLAF_FORTRAN_TARGET} INTERFACE_INCLUDE_DIRECTORIES) + if(CP2K_DLAF_FORTRAN_INCLUDE_DIRS) + list(APPEND CP2K_DLAF_INCLUDE_DIRS ${CP2K_DLAF_FORTRAN_INCLUDE_DIRS}) + else() + get_target_property(CP2K_DLAF_FORTRAN_LIBRARY ${CP2K_DLAF_FORTRAN_TARGET} + IMPORTED_LOCATION_RELEASE) + if(NOT CP2K_DLAF_FORTRAN_LIBRARY) + get_target_property(CP2K_DLAF_FORTRAN_LIBRARY ${CP2K_DLAF_FORTRAN_TARGET} + IMPORTED_LOCATION) + endif() + if(CP2K_DLAF_FORTRAN_LIBRARY) + get_filename_component(CP2K_DLAF_FORTRAN_PREFIX + "${CP2K_DLAF_FORTRAN_LIBRARY}" DIRECTORY) + get_filename_component(CP2K_DLAF_FORTRAN_PREFIX + "${CP2K_DLAF_FORTRAN_PREFIX}" DIRECTORY) + if(EXISTS "${CP2K_DLAF_FORTRAN_PREFIX}/include/dlaf_fortran.mod") + list(APPEND CP2K_DLAF_INCLUDE_DIRS + "${CP2K_DLAF_FORTRAN_PREFIX}/include") + endif() + endif() + endif() + list(REMOVE_DUPLICATES CP2K_DLAF_INCLUDE_DIRS) get_target_property(CP2K_DLAF_LINK_LIBRARIES DLAF::dlaf.prop INTERFACE_LINK_LIBRARIES) message("${CP2K_DLAF_INCLUDE_DIRS} ${CP2K_DLAF_LINK_LIBRARIES}") @@ -988,15 +1033,19 @@ if((CP2K_USE_ACCEL MATCHES "CUDA") OR (CP2K_USE_ACCEL MATCHES "HIP")) endif() if(CP2K_USE_CUSOLVER_MP) - message( - " - CUSolverMP: \n" - " - library: ${CP2K_CUSOLVER_MP_LINK_LIBRARIES} \n" - " - include: ${CP2K_CUSOLVER_MP_INCLUDE_DIRS} \n" - " - CAL library: ${CP2K_CAL_LINK_LIBRARIES} \n" - " - CAL include: ${CP2K_CAL_INCLUDE_DIRS} \n" - " - ucc library: ${CP2K_UCC_LINK_LIBRARIES} \n" - " - ucx library: ${CP2K_UCX_LINK_LIBRARIES} \n" - " - ucc include: ${CP2K_UCC_INCLUDE_DIRS} \n") + message(" - CUSolverMP: \n" + " - library: ${CP2K_CUSOLVER_MP_LINK_LIBRARIES} \n" + " - include: ${CP2K_CUSOLVER_MP_INCLUDE_DIRS} \n") + if(CP2K_CUSOLVERMP_USE_NCCL) + message(" - NCCL library: ${CP2K_NCCL_LINK_LIBRARIES} \n" + " - NCCL include: ${CP2K_NCCL_INCLUDE_DIRS} \n") + else() + message(" - CAL library: ${CP2K_CAL_LINK_LIBRARIES} \n" + " - CAL include: ${CP2K_CAL_INCLUDE_DIRS} \n") + endif() + message(" - ucc library: ${CP2K_UCC_LINK_LIBRARIES} \n" + " - ucx library: ${CP2K_UCX_LINK_LIBRARIES} \n" + " - ucc include: ${CP2K_UCC_INCLUDE_DIRS} \n") endif() if(CP2K_USE_LIBXC) diff --git a/cmake/CompilerConfiguration.cmake b/cmake/CompilerConfiguration.cmake index 6926d842c5..12efc126ee 100644 --- a/cmake/CompilerConfiguration.cmake +++ b/cmake/CompilerConfiguration.cmake @@ -5,8 +5,9 @@ #! SPDX-License-Identifier: GPL-2.0-or-later ! #!-------------------------------------------------------------------------------------------------! -set(CP2K_C_COMPILER_LIST "GNU;Intel;IntelLLVM;NAG;Cray;PGI;Clang;AppleClang") -set(CP2K_Fortran_COMPILER_LIST "GNU;Intel;IntelLLVM;NAG;Cray;PGI") +set(CP2K_C_COMPILER_LIST + "GNU;Intel;IntelLLVM;NAG;Cray;PGI;NVHPC;Clang;AppleClang") +set(CP2K_Fortran_COMPILER_LIST "GNU;Intel;IntelLLVM;NAG;Cray;PGI;NVHPC") if(NOT CMAKE_C_COMPILER_ID IN_LIST CP2K_C_COMPILER_LIST) message( @@ -157,7 +158,7 @@ add_compile_options( # Baseline add_compile_options( - "$<$:-Mfreeform;-Mextend;-Mallocatable=03>" + "$<$:-Mfreeform;-Mextend;-Mallocatable=03>" "$<$:-f2008;-free;-Warn=reallocation;-Warn=subnormal>" "$<$:-f;free;-M3105;-ME7212;-hnoacc;-M1234>" ) @@ -166,11 +167,11 @@ add_compile_options( # Release add_compile_options( - "$<$,$>:-fast>" + "$<$,$>:-fast>" "$<$,$>:-O2;-G2>" "$<$,$>:-gline>") add_compile_options( - "$<$,$>:-fast>" + "$<$,$>:-fast>" "$<$,$>:-O3;-g>" "$<$,$>:-O3>" "$<$,$>:-gline>" @@ -181,10 +182,10 @@ add_compile_options( # Debug add_compile_options( - "$<$,$>:-g>" + "$<$,$>:-g>" "$<$,$>:-G2>") add_compile_options( - "$<$,$>:-fast>" + "$<$,$>:-fast>" "$<$,$>:-O2;-g>" "$<$,$>:-G2>" "$<$,$>:-g;-C>" diff --git a/cmake/modules/FindSCALAPACK.cmake b/cmake/modules/FindSCALAPACK.cmake index 2b9c643657..1f33b08da6 100644 --- a/cmake/modules/FindSCALAPACK.cmake +++ b/cmake/modules/FindSCALAPACK.cmake @@ -30,6 +30,41 @@ if(NOT CP2K_CONFIG_PACKAGE) CP2K_SCALAPACK_LINK_LIBRARIES cp2k::BLAS::SCI::scalapack_link INTERFACE_LINK_LIBRARIES) set(CP2K_SCALAPACK_FOUND yes) + elseif(CP2K_SCALAPACK_VENDOR MATCHES "NVPL") + if(CP2K_BLAS_INTERFACE MATCHES "64bits") + set(_nvpl_int_type "ilp64") + else() + set(_nvpl_int_type "lp64") + endif() + + set(_nvpl_mpi_type "${CP2K_NVPL_SCALAPACK_MPI}") + if(_nvpl_mpi_type STREQUAL "auto") + if(MPI_Fortran_LIBRARY_VERSION_STRING MATCHES "Open MPI[^0-9]*([0-9]+)") + set(_nvpl_mpi_type "openmpi${CMAKE_MATCH_1}") + elseif(MPI_Fortran_LIBRARY_VERSION_STRING MATCHES "MPICH|HYDRA") + set(_nvpl_mpi_type "mpich") + else() + message( + FATAL_ERROR + "Could not determine the NVPL BLACS MPI interface. Set " + "CP2K_NVPL_SCALAPACK_MPI to mpich, openmpi3, openmpi4, or openmpi5." + ) + endif() + endif() + + find_package(nvpl REQUIRED COMPONENTS scalapack) + set(_nvpl_scalapack_target "nvpl::scalapack_${_nvpl_int_type}") + set(_nvpl_blacs_target "nvpl::blacs_${_nvpl_int_type}_${_nvpl_mpi_type}") + if(NOT TARGET "${_nvpl_scalapack_target}") + message( + FATAL_ERROR "NVPL ScaLAPACK target ${_nvpl_scalapack_target} not found") + endif() + if(NOT TARGET "${_nvpl_blacs_target}") + message(FATAL_ERROR "NVPL BLACS target ${_nvpl_blacs_target} not found") + endif() + set(CP2K_SCALAPACK_LINK_LIBRARIES + "${_nvpl_scalapack_target};${_nvpl_blacs_target}") + set(CP2K_SCALAPACK_FOUND yes) else() # if(CP2K_SCALAPACK_VENDOR MATCHES "GENERIC|auto") if(TARGET cp2k::BLAS::MKL::scalapack_link) message( diff --git a/cmake/modules/Finducc.cmake b/cmake/modules/Finducc.cmake index a1b7a1c433..ea73d478da 100644 --- a/cmake/modules/Finducc.cmake +++ b/cmake/modules/Finducc.cmake @@ -11,6 +11,7 @@ include(FindPackageHandleStandardArgs) include(cp2k_utils) cp2k_set_default_paths(UCC "ucc") +cp2k_set_default_paths(UCX "ucx") cp2k_find_libraries(UCC "ucc") cp2k_find_libraries(UCX "ucs") diff --git a/cmake/modules/cp2k_utils.cmake b/cmake/modules/cp2k_utils.cmake index 7de4e13f88..25c54d953b 100644 --- a/cmake/modules/cp2k_utils.cmake +++ b/cmake/modules/cp2k_utils.cmake @@ -13,8 +13,12 @@ function(cp2k_set_default_paths _varname _package_name) # find_library should work when ${PACKAGE_ROOT} is given to cmake # (-DPACKAGE_ROOT=bla) but I use only one variable syntax CP2K_PACKAGE_PREFIX set(CP2K_${_varname}_PREFIX_TMP "") - if(DEFINED ${_package_name}_ROOT) - set(CP2K_${_varname}_PREFIX_TMP "${${_varname}_ROOT}") + set(_cp2k_root_var "CP2K_${_varname}_ROOT") + set(_package_root_var "${_package_name}_ROOT") + if(DEFINED ${_cp2k_root_var}) + set(CP2K_${_varname}_PREFIX_TMP "${${_cp2k_root_var}}") + elseif(DEFINED ${_package_root_var}) + set(CP2K_${_varname}_PREFIX_TMP "${${_package_root_var}}") endif() # search common environment variables names diff --git a/docs/technologies/eigensolvers/cusolvermp.md b/docs/technologies/eigensolvers/cusolvermp.md index be115c1488..eb46dc0167 100644 --- a/docs/technologies/eigensolvers/cusolvermp.md +++ b/docs/technologies/eigensolvers/cusolvermp.md @@ -27,6 +27,16 @@ The `FindCuSolverMP.cmake` module tries to automatically deduce the [cuSOLVERmp] `cusolvermp.h` header file, and enables the `CP2K_CUSOLVERMP_USE_NCCL` CMake option in case it finds cuSOLVERmp >= 0.7. +## Generalized diagonalization + +The `DIRECT_GENERALIZED_DIAGONALIZATION` global input keyword enables direct generalized +diagonalization paths that avoid a CP2K-side Cholesky reduction where the selected eigensolver +supports them. With [cuSOLVERMp], this currently covers real symmetric and complex Hermitian +generalized eigenproblems through `cusolverMpSygvd`. + +When [cuSOLVERMp] uses the NCCL backend, CP2K requires at most one local MPI rank per visible GPU. +Use one MPI rank per GPU or reduce the number of local MPI ranks when only one GPU is visible. + [cal]: https://developer.download.nvidia.com/compute/cublasmp/redist/libcal/ [cusolvermp]: https://docs.nvidia.com/cuda/cusolvermp/ [nccl]: https://developer.nvidia.com/nccl diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a56a637fd5..f9624387b0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1712,6 +1712,10 @@ endif() # mix the target and variables pointing to the include directories. +if(CP2K_USE_DLAF) + include_directories(${CP2K_DLAF_INCLUDE_DIRS}) +endif() + include_directories( ${CMAKE_MPI_INCLUDE_DIRECTORIES} $<$:${CUDAToolkit_INCLUDE_DIRS}> @@ -1752,7 +1756,7 @@ target_link_libraries( $<$:CP2K::MIMIC::mcl> $<$:cp2k::Libint2::int2> $<$:cp2k::cosma> - $<$:DLAF::Fortran> + $<$:${CP2K_DLAF_FORTRAN_TARGET}> $<$:cp2k::greenx> $<$:cp2k::trexio::trexio> $<$:HDF5::HDF5 diff --git a/src/environment.F b/src/environment.F index 18486c3a02..15c2fe6193 100644 --- a/src/environment.F +++ b/src/environment.F @@ -34,6 +34,7 @@ MODULE environment FM_DIAG_TYPE_DLAF,& FM_DIAG_TYPE_ELPA,& FM_DIAG_TYPE_SCALAPACK,& + cusolver_n_min,& diag_finalize,& diag_init,& eps_check_diag_default @@ -545,7 +546,8 @@ CONTAINS CALL section_vals_val_get(global_section, "BLACS_GRID", i_val=globenv%blacs_grid_layout) CALL section_vals_val_get(global_section, "BLACS_REPEATABLE", l_val=globenv%blacs_repeatable) CALL section_vals_val_get(global_section, "PREFERRED_DIAG_LIBRARY", i_val=i_diag) - CALL section_vals_val_get(global_section, "CUSOLVER_GENERALIZED", l_val=globenv%cusolver_generalized) + CALL section_vals_val_get(global_section, "DIRECT_GENERALIZED_DIAGONALIZATION", & + l_val=globenv%direct_generalized_diagonalization) CALL section_vals_val_get(global_section, "PREFERRED_CHOLESKY_LIBRARY", i_val=i_cholesky) CALL section_vals_val_get(global_section, "PREFERRED_DGEMM_LIBRARY", i_val=i_dgemm) CALL section_vals_val_get(global_section, "EPS_CHECK_DIAG", r_val=globenv%eps_check_diag) @@ -830,6 +832,19 @@ CONTAINS globenv%dlaf_neigvec_min END IF + IF (globenv%diag_library == "cuSOLVER" .OR. globenv%diag_library == "ScaLAPACK" .OR. & + globenv%diag_library == "DLAF") THEN + WRITE (UNIT=output_unit, FMT="(T2,A,T71,L10)") & + start_section_label//"| Direct generalized diagonalization requested", & + globenv%direct_generalized_diagonalization + END IF + + IF (globenv%diag_library == "cuSOLVER") THEN + WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & + start_section_label//"| Minimum matrix size for cuSOLVER diagonalization", & + cusolver_n_min + END IF + IF (globenv%cholesky_library == "DLAF") THEN WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & start_section_label//"| Minimum matrix size for Cholesky decomposition with DLAF", & @@ -1119,7 +1134,8 @@ CONTAINS elpa_one_stage=globenv%elpa_one_stage, & dlaf_neigvec_min_input=globenv%dlaf_neigvec_min, & eps_check_diag_input=globenv%eps_check_diag, & - cusolver_generalized_input=globenv%cusolver_generalized) + direct_generalized_diagonalization_input= & + globenv%direct_generalized_diagonalization) IF (fallback_applied) THEN message = "Diagonalization library "//TRIM(globenv%diag_library)// & diff --git a/src/fm/cp_cfm_diag.F b/src/fm/cp_cfm_diag.F index 0f52050e68..d858ef39c1 100644 --- a/src/fm/cp_cfm_diag.F +++ b/src/fm/cp_cfm_diag.F @@ -12,23 +12,37 @@ !> \author Joost VandeVondele (2003-09) ! ************************************************************************************************** MODULE cp_cfm_diag + USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, & cp_cfm_column_scale, & cp_cfm_scale, & cp_cfm_triangular_invert, & cp_cfm_triangular_multiply - USE cp_cfm_types, ONLY: cp_cfm_get_info, & + USE cp_cfm_types, ONLY: cp_cfm_create, & + cp_cfm_get_info, & + cp_cfm_release, & cp_cfm_set_element, & cp_cfm_to_cfm, & cp_cfm_type + USE cp_fm_diag, ONLY: diag_check_requested, & + diag_check_warning_threshold, & + diag_type, & + direct_generalized_diagonalization, & + cusolver_n_min, & + FM_DIAG_TYPE_CUSOLVER, & + FM_DIAG_TYPE_SCALAPACK + USE cp_fm_cusolver_api, ONLY: cp_cfm_general_cusolver #if defined(__DLAF) USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, & cp_cfm_diag_dlaf USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid - USE cp_fm_diag, ONLY: diag_type, dlaf_neigvec_min, FM_DIAG_TYPE_DLAF + USE cp_fm_diag, ONLY: dlaf_neigvec_min, FM_DIAG_TYPE_DLAF #endif - USE kinds, ONLY: dp + USE cp_log_handling, ONLY: cp_to_string + USE kinds, ONLY: default_string_length, & + dp + USE machine, ONLY: default_output_unit USE mathconstants, ONLY: z_one, & z_zero #if defined (__HAS_IEEE_EXCEPTIONS) @@ -177,6 +191,109 @@ CONTAINS END SUBROUTINE cp_cfm_heevd_base +! ************************************************************************************************** +!> \brief Check C^H*S*C = I for a generalized complex eigenvalue problem. +!> \param overlap original overlap matrix S; used as work matrix and overwritten +!> \param eigenvectors eigenvectors C to be checked +!> \param scratch work matrix +!> \param nvec ... +! ************************************************************************************************** + SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec) + + TYPE(cp_cfm_type), INTENT(IN) :: eigenvectors + TYPE(cp_cfm_type), INTENT(INOUT) :: overlap, scratch + INTEGER, INTENT(IN) :: nvec + + CHARACTER(LEN=*), PARAMETER :: routineN = 'check_generalized_diag' + + CHARACTER(LEN=default_string_length) :: diag_type_name + COMPLEX(KIND=dp) :: gold, test + INTEGER :: handle, i, j, ncol, nrow, output_unit + REAL(KIND=dp) :: eps, eps_abort, eps_warning +#if defined(__parallel) + TYPE(cp_blacs_env_type), POINTER :: context + INTEGER :: il, jl, ipcol, iprow, & + mypcol, myprow, npcol, nprow + INTEGER, DIMENSION(9) :: desca +#endif + + CALL timeset(routineN, handle) + + IF (.NOT. diag_check_requested()) THEN + CALL timestop(handle) + RETURN + END IF + + output_unit = default_output_unit + eps_warning = diag_check_warning_threshold() + eps_abort = 10.0_dp*eps_warning + + nrow = eigenvectors%matrix_struct%nrow_global + ncol = MIN(eigenvectors%matrix_struct%ncol_global, nvec) + + CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch) + CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap) + + gold = z_zero + test = z_zero + eps = 0.0_dp + +#if defined(__parallel) + context => overlap%matrix_struct%context + myprow = context%mepos(1) + mypcol = context%mepos(2) + nprow = context%num_pe(1) + npcol = context%num_pe(2) + desca(:) = overlap%matrix_struct%descriptor(:) + outer: DO j = 1, ncol + DO i = 1, ncol + CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol) + IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN + gold = MERGE(z_zero, z_one, i /= j) + test = overlap%local_data(il, jl) + eps = ABS(test - gold) + IF (eps > eps_warning) EXIT outer + END IF + END DO + END DO outer +#else + outer: DO j = 1, ncol + DO i = 1, ncol + gold = MERGE(z_zero, z_one, i /= j) + test = overlap%local_data(i, j) + eps = ABS(test - gold) + IF (eps > eps_warning) EXIT outer + END DO + END DO outer +#endif + + IF (eps > eps_warning) THEN + IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN + diag_type_name = "HEGVX" + ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER) THEN + diag_type_name = "CUSOLVER" +#if defined(__DLAF) + ELSE IF (diag_type == FM_DIAG_TYPE_DLAF) THEN + diag_type_name = "DLAF" +#endif + ELSE + diag_type_name = "generalized eigensolver" + END IF + WRITE (UNIT=output_unit, FMT="(/,T2,A,/,T2,A,I0,A,I0,A,ES10.3,/,T2,A,F0.0,A,ES10.3)") & + "The generalized eigenvectors returned by "//TRIM(diag_type_name)//" are not S-orthonormal", & + "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, & + "The deviation from the expected value ", REAL(gold, KIND=dp), " is", eps + IF (eps > eps_abort) THEN + CPABORT("ERROR in "//routineN//": Check of generalized matrix diagonalization failed") + ELSE + CPWARN("Check of generalized matrix diagonalization failed in routine "//routineN) + END IF + END IF + + CALL timestop(handle) + + END SUBROUTINE check_generalized_diag + ! ************************************************************************************************** !> \brief General Eigenvalue Problem AX = BXE !> Single option version: Cholesky decomposition of B @@ -197,16 +314,35 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig' INTEGER :: handle, nao, nmo + LOGICAL :: check_eigenvectors REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals + TYPE(cp_cfm_type) :: overlap_check, scratch_check CALL timeset(routineN, handle) CALL cp_cfm_get_info(amatrix, nrow_global=nao) ALLOCATE (evals(nao)) nmo = SIZE(eigenvalues) + check_eigenvectors = diag_check_requested() + IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. direct_generalized_diagonalization .AND. & + nao >= cusolver_n_min) THEN + ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side + ! Cholesky reduction. + IF (check_eigenvectors) THEN + CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_cfm_to_cfm(bmatrix, overlap_check) + END IF + CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals) + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_cfm_release(scratch_check) + CALL cp_cfm_release(overlap_check) + END IF #if defined(__DLAF) - IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN + ELSE IF (diag_type == FM_DIAG_TYPE_DLAF .AND. direct_generalized_diagonalization .AND. & + nao >= dlaf_neigvec_min) THEN ! Initialize DLA-Future on-demand; if already initialized, does nothing CALL cp_dlaf_initialize() @@ -216,9 +352,35 @@ CONTAINS CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle()) ! Use DLA-Future generalized eigenvalue solver for large matrices + IF (check_eigenvectors) THEN + CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_cfm_to_cfm(bmatrix, overlap_check) + END IF CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals) - ELSE + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_cfm_release(scratch_check) + CALL cp_cfm_release(overlap_check) + END IF #endif +#if defined(__parallel) + ELSE IF (diag_type == FM_DIAG_TYPE_SCALAPACK .AND. direct_generalized_diagonalization) THEN + ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side + ! Cholesky reduction. + IF (check_eigenvectors) THEN + CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_cfm_to_cfm(bmatrix, overlap_check) + END IF + CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals) + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_cfm_release(scratch_check) + CALL cp_cfm_release(overlap_check) + END IF +#endif + ELSE ! Cholesky decompose S=U(T)U CALL cp_cfm_cholesky_decompose(bmatrix) ! Invert to get U^(-1) @@ -230,9 +392,7 @@ CONTAINS CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals) ! Restore vectors C = U^(-1) * C* CALL cp_cfm_triangular_multiply(bmatrix, work) -#if defined(__DLAF) END IF -#endif CALL cp_cfm_to_cfm(work, eigenvectors, nmo) eigenvalues(1:nmo) = evals(1:nmo) @@ -243,6 +403,123 @@ CONTAINS END SUBROUTINE cp_cfm_geeig +! ************************************************************************************************** +!> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PZHEGVX. +!> \param amatrix ... +!> \param bmatrix ... +!> \param eigenvectors ... +!> \param eigenvalues ... +! ************************************************************************************************** + SUBROUTINE cp_cfm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues) + + TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues + + CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_scalapack' + +#if defined(__parallel) + REAL(KIND=dp), PARAMETER :: orfac = -1.0_dp, & + vl = 0.0_dp, & + vu = 0.0_dp + + COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work + COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z + INTEGER :: handle, info, liwork, lwork, lrwork, & + m, n, nb, neig, npcol, nprow, nz + INTEGER, DIMENSION(9) :: desca, descb, descz + INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork + REAL(KIND=dp) :: abstol + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: gap, rwork, w + + INTEGER :: mq0, nn, np0, npe + INTEGER, EXTERNAL :: iceil, numroc + REAL(KIND=dp), EXTERNAL :: dlamch +#if defined (__HAS_IEEE_EXCEPTIONS) + LOGICAL, DIMENSION(5) :: halt +#endif +#else + INTEGER :: handle +#endif + + CALL timeset(routineN, handle) + +#if defined(__parallel) + n = amatrix%matrix_struct%nrow_global + neig = MIN(SIZE(eigenvalues), n) + + IF (neig == 0) THEN + CALL timestop(handle) + RETURN + END IF + + IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN + CPABORT("ERROR in "//routineN//": Invalid blocksize (no square blocks) found") + END IF + + a => amatrix%local_data + b => bmatrix%local_data + z => eigenvectors%local_data + desca(:) = amatrix%matrix_struct%descriptor(:) + descb(:) = bmatrix%matrix_struct%descriptor(:) + descz(:) = eigenvectors%matrix_struct%descriptor(:) + + nprow = amatrix%matrix_struct%context%num_pe(1) + npcol = amatrix%matrix_struct%context%num_pe(2) + npe = nprow*npcol + nb = amatrix%matrix_struct%nrow_block + nn = MAX(n, nb, 2) + np0 = numroc(nn, nb, 0, 0, nprow) + mq0 = MAX(numroc(nn, nb, 0, 0, npcol), nb) + + lwork = n + (np0 + mq0 + nb)*nb + lrwork = 4*n + MAX(5*nn, np0*mq0) + iceil(neig, npe)*nn + MAX(0, neig - 1)*n + liwork = 6*MAX(n, npe + 1, 4) + + ALLOCATE (gap(npe)) + gap = 0.0_dp + ALLOCATE (iclustr(2*npe)) + iclustr = 0 + ALLOCATE (ifail(n)) + ifail = 0 + ALLOCATE (iwork(liwork)) + ALLOCATE (rwork(lrwork)) + ALLOCATE (w(n)) + ALLOCATE (work(lwork)) + + abstol = 2.0_dp*dlamch("S") + +#if defined (__HAS_IEEE_EXCEPTIONS) + CALL ieee_get_halting_mode(IEEE_ALL, halt) + CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.) +#endif + CALL pzhegvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, & + vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, & + work(1), lwork, rwork(1), lrwork, iwork(1), liwork, ifail(1), & + iclustr(1), gap(1), info) +#if defined (__HAS_IEEE_EXCEPTIONS) + CALL ieee_set_halting_mode(IEEE_ALL, halt) +#endif + + IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN + CPABORT("ERROR in PZHEGVX (ScaLAPACK), info="//TRIM(cp_to_string(info))) + END IF + + eigenvalues(:) = 0.0_dp + eigenvalues(1:neig) = w(1:neig) + + DEALLOCATE (gap, iclustr, ifail, iwork, rwork, w, work) +#else + MARK_USED(amatrix) + MARK_USED(bmatrix) + MARK_USED(eigenvectors) + MARK_USED(eigenvalues) + CPABORT("ERROR in "//routineN//": PZHEGVX requested without ScaLAPACK support") +#endif + + CALL timestop(handle) + + END SUBROUTINE cp_cfm_geeig_scalapack + ! ************************************************************************************************** !> \brief General Eigenvalue Problem AX = BXE !> Use canonical orthogonalization diff --git a/src/fm/cp_fm_cusolver.c b/src/fm/cp_fm_cusolver.c index 55d7339714..b4e6ce7497 100644 --- a/src/fm/cp_fm_cusolver.c +++ b/src/fm/cp_fm_cusolver.c @@ -9,10 +9,12 @@ #include "../offload/offload_library.h" #include +#include #include #include #include #include +#include #include #include @@ -51,6 +53,60 @@ } \ } while (0) +/******************************************************************************* + * \brief Check that local MPI ranks can be mapped one-to-one to GPUs. + ******************************************************************************/ +static void check_nccl_rank_device_mapping(MPI_Comm comm, + const int local_device) { + int device_count = 0; + CUDA_CHECK(cudaGetDeviceCount(&device_count)); + + MPI_Comm node_comm; + MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &node_comm); + + int rank, local_rank, local_nranks; + MPI_Comm_rank(comm, &rank); + MPI_Comm_rank(node_comm, &local_rank); + MPI_Comm_size(node_comm, &local_nranks); + + if (device_count <= 0) { + if (local_rank == 0) { + fprintf(stderr, + "ERROR: cuSOLVERMp with NCCL requires at least one visible GPU " + "on every node. No CUDA device is visible on the node containing " + "MPI rank %d.\n", + rank); + } + MPI_Comm_free(&node_comm); + MPI_Abort(comm, EXIT_FAILURE); + } + + if (local_device < 0 || local_device >= device_count) { + fprintf(stderr, + "ERROR: cuSOLVERMp with NCCL selected invalid CUDA device %d on " + "MPI rank %d; %d CUDA device(s) are visible on this node.\n", + local_device, rank, device_count); + MPI_Comm_free(&node_comm); + MPI_Abort(comm, EXIT_FAILURE); + } + + if (local_nranks > device_count) { + if (local_rank == 0) { + fprintf(stderr, + "ERROR: cuSOLVERMp with NCCL in CP2K requires at most one local " + "MPI rank per visible GPU. This node has %d local MPI ranks and " + "%d visible GPU(s). Use fewer MPI ranks per node, expose more " + "GPUs with CUDA_VISIBLE_DEVICES, or choose another " + "PREFERRED_DIAG_LIBRARY.\n", + local_nranks, device_count); + } + MPI_Comm_free(&node_comm); + MPI_Abort(comm, EXIT_FAILURE); + } + + MPI_Comm_free(&node_comm); +} + #else /******************************************************************************* * \brief Decode given cal error. @@ -215,7 +271,6 @@ void cp_fm_diag_cusolver(const int fortran_comm, const int matrix_desc[9], const int mypcol, const int n, const double *matrix, double *eigenvectors, double *eigenvalues) { - offload_activate_chosen_device(); const int local_device = offload_get_chosen_device(); MPI_Comm comm = MPI_Comm_f2c(fortran_comm); @@ -223,6 +278,11 @@ void cp_fm_diag_cusolver(const int fortran_comm, const int matrix_desc[9], MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &nranks); +#if defined(__CUSOLVERMP_NCCL) + check_nccl_rank_device_mapping(comm, local_device); +#endif + offload_activate_chosen_device(); + #if defined(__CUSOLVERMP_NCCL) // Create NCCL communicator. ncclUniqueId nccl_id; @@ -376,7 +436,6 @@ void cp_fm_diag_cusolver_sygvd(const int fortran_comm, const double *aMatrix, const double *bMatrix, double *eigenvectors, double *eigenvalues) { - offload_activate_chosen_device(); const int local_device = offload_get_chosen_device(); MPI_Comm comm = MPI_Comm_f2c(fortran_comm); @@ -384,6 +443,11 @@ void cp_fm_diag_cusolver_sygvd(const int fortran_comm, MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &nranks); +#if defined(__CUSOLVERMP_NCCL) + check_nccl_rank_device_mapping(comm, local_device); +#endif + offload_activate_chosen_device(); + #if defined(__CUSOLVERMP_NCCL) // Create NCCL communicator. ncclUniqueId nccl_id; @@ -443,7 +507,7 @@ void cp_fm_diag_cusolver_sygvd(const int fortran_comm, // Ensure consistency in block sizes, sources, and leading dimensions assert(mb_a == mb_b && nb_a == nb_b); assert(rsrc_a == rsrc_b && csrc_a == csrc_b); - (void)ldB; // Suppress unused variable warning + assert(ldA == ldB); const int np_a = cusolverMpNUMROC(n, mb_a, myprow, rsrc_a, nprow); const int nq_a = cusolverMpNUMROC(n, nb_a, mypcol, csrc_a, npcol); @@ -563,6 +627,190 @@ void cp_fm_diag_cusolver_sygvd(const int fortran_comm, MPI_Barrier(comm); // Synchronize MPI ranks } + +/******************************************************************************* + * \brief Driver routine to solve A*x = lambda*B*x with cuSOLVERMp hegvd. + ******************************************************************************/ +void cp_cfm_diag_cusolver_hegvd( + const int fortran_comm, const int a_matrix_desc[9], + const int b_matrix_desc[9], const int nprow, const int npcol, + const int myprow, const int mypcol, const int n, + const cuDoubleComplex *aMatrix, const cuDoubleComplex *bMatrix, + cuDoubleComplex *eigenvectors, double *eigenvalues) { + + const int local_device = offload_get_chosen_device(); + + MPI_Comm comm = MPI_Comm_f2c(fortran_comm); + int rank, nranks; + MPI_Comm_rank(comm, &rank); + MPI_Comm_size(comm, &nranks); + +#if defined(__CUSOLVERMP_NCCL) + check_nccl_rank_device_mapping(comm, local_device); +#endif + offload_activate_chosen_device(); + +#if defined(__CUSOLVERMP_NCCL) + ncclUniqueId nccl_id; + if (rank == 0) { + NCCL_CHECK(ncclGetUniqueId(&nccl_id)); + } + MPI_Bcast(&nccl_id, sizeof(nccl_id), MPI_BYTE, 0, comm); + + ncclComm_t nccl_comm; + NCCL_CHECK(ncclCommInitRank(&nccl_comm, nranks, nccl_id, rank)); +#else + cal_comm_t cal_comm = NULL; + cal_comm_create_params_t params; + params.allgather = &allgather; + params.req_test = &req_test; + params.req_free = &req_free; + params.data = &comm; + params.rank = rank; + params.nranks = nranks; + params.local_device = local_device; + CAL_CHECK(cal_comm_create(params, &cal_comm)); +#endif + + cudaStream_t stream = NULL; + CUDA_CHECK(cudaStreamCreate(&stream)); + + cusolverMpHandle_t cusolvermp_handle = NULL; + CUSOLVER_CHECK(cusolverMpCreate(&cusolvermp_handle, local_device, stream)); + + cusolverMpGrid_t grid = NULL; +#if defined(__CUSOLVERMP_NCCL) + CUSOLVER_CHECK(cusolverMpCreateDeviceGrid(cusolvermp_handle, &grid, nccl_comm, + nprow, npcol, + CUSOLVERMP_GRID_MAPPING_ROW_MAJOR)); +#else + CUSOLVER_CHECK(cusolverMpCreateDeviceGrid(cusolvermp_handle, &grid, cal_comm, + nprow, npcol, + CUSOLVERMP_GRID_MAPPING_ROW_MAJOR)); +#endif + + const int mb_a = a_matrix_desc[4]; + const int nb_a = a_matrix_desc[5]; + const int rsrc_a = a_matrix_desc[6]; + const int csrc_a = a_matrix_desc[7]; + const int ldA = a_matrix_desc[8]; + + const int mb_b = b_matrix_desc[4]; + const int nb_b = b_matrix_desc[5]; + const int rsrc_b = b_matrix_desc[6]; + const int csrc_b = b_matrix_desc[7]; + const int ldB = b_matrix_desc[8]; + + assert(mb_a == mb_b && nb_a == nb_b); + assert(rsrc_a == rsrc_b && csrc_a == csrc_b); + assert(ldA == ldB); + + const int np_a = cusolverMpNUMROC(n, mb_a, myprow, rsrc_a, nprow); + const int nq_a = cusolverMpNUMROC(n, nb_a, mypcol, csrc_a, npcol); + assert(np_a == ldA); + + const cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER; + const cusolverEigType_t itype = CUSOLVER_EIG_TYPE_1; + const cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; + const cudaDataType_t data_type = CUDA_C_64F; + + cusolverMpMatrixDescriptor_t descrA = NULL; + cusolverMpMatrixDescriptor_t descrB = NULL; + cusolverMpMatrixDescriptor_t descrZ = NULL; + + CUSOLVER_CHECK(cusolverMpCreateMatrixDesc(&descrA, grid, data_type, n, n, + mb_a, nb_a, rsrc_a, csrc_a, ldA)); + CUSOLVER_CHECK(cusolverMpCreateMatrixDesc(&descrB, grid, data_type, n, n, + mb_b, nb_b, rsrc_b, csrc_b, ldA)); + CUSOLVER_CHECK(cusolverMpCreateMatrixDesc(&descrZ, grid, data_type, n, n, + mb_a, nb_a, rsrc_a, csrc_a, ldA)); + + cuDoubleComplex *dev_A = NULL, *dev_B = NULL; + size_t matrix_local_size = ldA * nq_a * sizeof(cuDoubleComplex); + CUDA_CHECK(cudaMalloc((void **)&dev_A, matrix_local_size)); + CUDA_CHECK(cudaMalloc((void **)&dev_B, matrix_local_size)); + + CUDA_CHECK(cudaMemcpyAsync(dev_A, aMatrix, matrix_local_size, + cudaMemcpyHostToDevice, stream)); + CUDA_CHECK(cudaMemcpyAsync(dev_B, bMatrix, matrix_local_size, + cudaMemcpyHostToDevice, stream)); + + cuDoubleComplex *dev_Z = NULL; + double *eigenvalues_dev = NULL; + CUDA_CHECK(cudaMalloc((void **)&dev_Z, matrix_local_size)); + CUDA_CHECK(cudaMalloc((void **)&eigenvalues_dev, n * sizeof(double))); + + size_t work_dev_size = 0, work_host_size = 0; + const int64_t ia = 1, ja = 1, ib = 1, jb = 1, iz = 1, jz = 1; + const int64_t m = (int64_t)n; + + cusolverStatus_t status_bufsize = cusolverMpSygvd_bufferSize( + cusolvermp_handle, itype, jobz, uplo, m, ia, ja, descrA, ib, jb, descrB, + iz, jz, descrZ, data_type, &work_dev_size, &work_host_size); + if (status_bufsize != CUSOLVER_STATUS_SUCCESS) { + fprintf(stderr, "ERROR: cusolverMpSygvd_bufferSize failed with status=%d\n", + (int)status_bufsize); + abort(); + } + + void *work_dev = NULL, *work_host = NULL; + CUDA_CHECK(cudaMalloc(&work_dev, work_dev_size)); + CUDA_CHECK(cudaMallocHost(&work_host, work_host_size)); + + int *info_dev = NULL; + CUDA_CHECK(cudaMalloc((void **)&info_dev, sizeof(int))); + CUDA_CHECK(cudaMemset(info_dev, 0, sizeof(int))); + + cusolverStatus_t status_sygvd = cusolverMpSygvd( + cusolvermp_handle, itype, jobz, uplo, m, dev_A, ia, ja, descrA, dev_B, ib, + jb, descrB, eigenvalues_dev, dev_Z, iz, jz, descrZ, data_type, work_dev, + work_dev_size, work_host, work_host_size, info_dev); + if (status_sygvd != CUSOLVER_STATUS_SUCCESS) { + fprintf(stderr, "ERROR: cusolverMpSygvd failed with status=%d\n", + (int)status_sygvd); + abort(); + } + + CUDA_CHECK(cudaStreamSynchronize(stream)); +#if !defined(__CUSOLVERMP_NCCL) + CAL_CHECK(cal_stream_sync(cal_comm, stream)); +#endif + + int info; + CUDA_CHECK(cudaMemcpy(&info, info_dev, sizeof(int), cudaMemcpyDeviceToHost)); + if (info != 0) { + fprintf(stderr, "ERROR: cusolverMpSygvd failed with info = %d\n", info); + abort(); + } + + CUDA_CHECK(cudaMemcpyAsync(eigenvectors, dev_Z, matrix_local_size, + cudaMemcpyDeviceToHost, stream)); + CUDA_CHECK(cudaMemcpyAsync(eigenvalues, eigenvalues_dev, n * sizeof(double), + cudaMemcpyDeviceToHost, stream)); + + CUDA_CHECK(cudaStreamSynchronize(stream)); + + CUDA_CHECK(cudaFree(dev_A)); + CUDA_CHECK(cudaFree(dev_B)); + CUDA_CHECK(cudaFree(dev_Z)); + CUDA_CHECK(cudaFree(eigenvalues_dev)); + CUDA_CHECK(cudaFree(info_dev)); + CUDA_CHECK(cudaFree(work_dev)); + CUDA_CHECK(cudaFreeHost(work_host)); + CUSOLVER_CHECK(cusolverMpDestroyMatrixDesc(descrA)); + CUSOLVER_CHECK(cusolverMpDestroyMatrixDesc(descrB)); + CUSOLVER_CHECK(cusolverMpDestroyMatrixDesc(descrZ)); + CUSOLVER_CHECK(cusolverMpDestroyGrid(grid)); + CUSOLVER_CHECK(cusolverMpDestroy(cusolvermp_handle)); + CUDA_CHECK(cudaStreamDestroy(stream)); +#if defined(__CUSOLVERMP_NCCL) + NCCL_CHECK(ncclCommDestroy(nccl_comm)); +#else + CAL_CHECK(cal_comm_destroy(cal_comm)); +#endif + + MPI_Barrier(comm); +} #endif // EOF diff --git a/src/fm/cp_fm_cusolver_api.F b/src/fm/cp_fm_cusolver_api.F index 34a7cf4879..79a2acc842 100644 --- a/src/fm/cp_fm_cusolver_api.F +++ b/src/fm/cp_fm_cusolver_api.F @@ -11,8 +11,10 @@ ! ************************************************************************************************** MODULE cp_fm_cusolver_api USE ISO_C_BINDING, ONLY: C_DOUBLE,& + C_DOUBLE_COMPLEX,& C_INT USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_cfm_types, ONLY: cp_cfm_type USE cp_fm_types, ONLY: cp_fm_type USE kinds, ONLY: dp #include "../base/base_uses.f90" @@ -22,6 +24,7 @@ MODULE cp_fm_cusolver_api PRIVATE PUBLIC :: cp_fm_diag_cusolver + PUBLIC :: cp_cfm_general_cusolver PUBLIC :: cp_fm_general_cusolver CONTAINS @@ -178,4 +181,83 @@ CONTAINS CALL timestop(handle) END SUBROUTINE cp_fm_general_cusolver +! ************************************************************************************************** +!> \brief Driver routine to solve generalized complex eigenvalue problem A*x = lambda*B*x with +!> cuSOLVERMp. +!> \param aMatrix the first matrix for the generalized eigenvalue problem +!> \param bMatrix the second matrix for the generalized eigenvalue problem +!> \param eigenvectors eigenvectors of the input matrix +!> \param eigenvalues eigenvalues of the input matrix +! ************************************************************************************************** + SUBROUTINE cp_cfm_general_cusolver(aMatrix, bMatrix, eigenvectors, eigenvalues) + USE ISO_C_BINDING, ONLY: C_DOUBLE, C_DOUBLE_COMPLEX, C_INT + TYPE(cp_cfm_type), INTENT(IN) :: aMatrix, bMatrix, eigenvectors + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues + + CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_general_cusolver' + + INTEGER(kind=C_INT) :: handle, n, nmo + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_buffer + TYPE(cp_blacs_env_type), POINTER :: context + INTERFACE + SUBROUTINE cp_cfm_general_cusolver_c(fortran_comm, a_matrix_desc, b_matrix_desc, & + nprow, npcol, myprow, mypcol, & + n, aMatrix, bMatrix, eigenvectors, eigenvalues) & + BIND(C, name="cp_cfm_diag_cusolver_hegvd") + IMPORT :: C_DOUBLE, C_DOUBLE_COMPLEX, C_INT + INTEGER(kind=C_INT), VALUE :: fortran_comm + INTEGER(kind=C_INT), DIMENSION(*) :: a_matrix_desc, b_matrix_desc + INTEGER(kind=C_INT), VALUE :: nprow + INTEGER(kind=C_INT), VALUE :: npcol + INTEGER(kind=C_INT), VALUE :: myprow + INTEGER(kind=C_INT), VALUE :: mypcol + INTEGER(kind=C_INT), VALUE :: n + COMPLEX(kind=C_DOUBLE_COMPLEX), DIMENSION(*) :: aMatrix + COMPLEX(kind=C_DOUBLE_COMPLEX), DIMENSION(*) :: bMatrix + COMPLEX(kind=C_DOUBLE_COMPLEX), DIMENSION(*) :: eigenvectors + REAL(kind=C_DOUBLE), DIMENSION(*) :: eigenvalues + END SUBROUTINE cp_cfm_general_cusolver_c + END INTERFACE + + CALL timeset(routineN, handle) + +#if defined(__CUSOLVERMP) + n = INT(aMatrix%matrix_struct%nrow_global, C_INT) + context => aMatrix%matrix_struct%context + + ALLOCATE (eigenvalues_buffer(n)) + + CALL cp_cfm_general_cusolver_c( & + fortran_comm=INT(aMatrix%matrix_struct%para_env%get_handle(), C_INT), & + a_matrix_desc=INT(aMatrix%matrix_struct%descriptor, C_INT), & + b_matrix_desc=INT(bMatrix%matrix_struct%descriptor, C_INT), & + nprow=INT(context%num_pe(1), C_INT), & + npcol=INT(context%num_pe(2), C_INT), & + myprow=INT(context%mepos(1), C_INT), & + mypcol=INT(context%mepos(2), C_INT), & + n=n, & + aMatrix=aMatrix%local_data, & + bMatrix=bMatrix%local_data, & + eigenvectors=eigenvectors%local_data, & + eigenvalues=eigenvalues_buffer) + + nmo = SIZE(eigenvalues) + eigenvalues(1:nmo) = eigenvalues_buffer(1:nmo) + + DEALLOCATE (eigenvalues_buffer) +#else + MARK_USED(aMatrix) + MARK_USED(bMatrix) + MARK_USED(eigenvectors) + eigenvalues = 0.0_dp + MARK_USED(n) + MARK_USED(nmo) + MARK_USED(eigenvalues_buffer) + MARK_USED(context) + CPABORT("CP2K compiled without the cuSOLVERMp library.") +#endif + + CALL timestop(handle) + END SUBROUTINE cp_cfm_general_cusolver + END MODULE cp_fm_cusolver_api diff --git a/src/fm/cp_fm_diag.F b/src/fm/cp_fm_diag.F index 37474a53f2..9ebbded7ac 100644 --- a/src/fm/cp_fm_diag.F +++ b/src/fm/cp_fm_diag.F @@ -65,6 +65,7 @@ MODULE cp_fm_diag dp USE machine, ONLY: default_output_unit, & m_memory + USE parallel_gemm_api, ONLY: parallel_gemm #if defined (__parallel) USE message_passing, ONLY: mp_comm_type #endif @@ -88,12 +89,15 @@ MODULE cp_fm_diag ! Minimum number of eigenvectors for the use of the ELPA eigensolver. ! The ScaLAPACK eigensolver is used as fallback for all smaller cases. INTEGER, SAVE :: elpa_neigvec_min = 0 + ! Minimum matrix size for the use of the cuSOLVERMp eigensolver. + ! Smaller matrices use the ScaLAPACK fallback to avoid GPU launch overheads. + INTEGER, PARAMETER, PUBLIC :: cusolver_n_min = 64 #if defined(__DLAF) ! Minimum number of eigenvectors for the use of the DLAF eigensolver. ! The ScaLAPACK eigensolver is used as fallback for all smaller cases. INTEGER, SAVE, PUBLIC :: dlaf_neigvec_min = 0 #endif - LOGICAL, SAVE, PUBLIC :: cusolver_generalized = .TRUE. + LOGICAL, SAVE, PUBLIC :: direct_generalized_diagonalization = .FALSE. ! Threshold value for the orthonormality check of the eigenvectors obtained ! after a diagonalization. A negative value disables the check. REAL(KIND=dp), SAVE :: eps_check_diag = -1.0_dp @@ -120,6 +124,8 @@ MODULE cp_fm_diag cp_fm_svd, & cp_fm_geeig, & cp_fm_geeig_canon, & + diag_check_requested, & + diag_check_warning_threshold, & diag_init, & diag_finalize @@ -139,21 +145,21 @@ CONTAINS !> \param elpa_one_stage logical that enables the one-stage solver !> \param dlaf_neigvec_min_input ... !> \param eps_check_diag_input ... -!> \param cusolver_generalized_input ... +!> \param direct_generalized_diagonalization_input ... !> \par History !> - Add support for DLA-Future (05.09.2023, RMeli) !> \author MI 11.2013 ! ************************************************************************************************** SUBROUTINE diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_neigvec_min_input, elpa_qr, & elpa_print, elpa_one_stage, dlaf_neigvec_min_input, eps_check_diag_input, & - cusolver_generalized_input) + direct_generalized_diagonalization_input) CHARACTER(LEN=*), INTENT(IN) :: diag_lib LOGICAL, INTENT(OUT) :: fallback_applied INTEGER, INTENT(IN) :: elpa_kernel, elpa_neigvec_min_input LOGICAL, INTENT(IN) :: elpa_qr, elpa_print, elpa_one_stage INTEGER, INTENT(IN) :: dlaf_neigvec_min_input REAL(KIND=dp), INTENT(IN) :: eps_check_diag_input - LOGICAL, INTENT(IN), OPTIONAL :: cusolver_generalized_input + LOGICAL, INTENT(IN), OPTIONAL :: direct_generalized_diagonalization_input LOGICAL, SAVE :: initialized = .FALSE. @@ -200,10 +206,10 @@ CONTAINS elpa_neigvec_min = elpa_neigvec_min_input eps_check_diag = eps_check_diag_input - IF (PRESENT(cusolver_generalized_input)) THEN - cusolver_generalized = cusolver_generalized_input + IF (PRESENT(direct_generalized_diagonalization_input)) THEN + direct_generalized_diagonalization = direct_generalized_diagonalization_input ELSE - cusolver_generalized = .TRUE. + direct_generalized_diagonalization = .FALSE. END IF END SUBROUTINE diag_init @@ -260,7 +266,7 @@ CONTAINS END IF ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER) THEN - IF (matrix%matrix_struct%nrow_global < 64) THEN + IF (matrix%matrix_struct%nrow_global < cusolver_n_min) THEN ! We don't trust cuSolver with very small matrices. CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info) ELSE @@ -285,6 +291,35 @@ CONTAINS END SUBROUTINE choose_eigv_solver +! ************************************************************************************************** +!> \brief Return whether diagonalization checks should be performed. +!> \return ... +! ************************************************************************************************** + FUNCTION diag_check_requested() RESULT(check_requested) + LOGICAL :: check_requested + +#if defined(__CHECK_DIAG) + check_requested = .TRUE. +#else + check_requested = eps_check_diag >= 0.0_dp +#endif + + END FUNCTION diag_check_requested + +! ************************************************************************************************** +!> \brief Return the warning threshold for diagonalization checks. +!> \return ... +! ************************************************************************************************** + FUNCTION diag_check_warning_threshold() RESULT(eps_warning) + REAL(KIND=dp) :: eps_warning + + eps_warning = eps_check_diag_default + IF (eps_check_diag >= 0.0_dp) THEN + eps_warning = eps_check_diag + END IF + + END FUNCTION diag_check_warning_threshold + ! ************************************************************************************************** !> \brief Check result of diagonalization, i.e. the orthonormality of the eigenvectors !> \param matrix Work matrix @@ -312,20 +347,8 @@ CONTAINS CALL timeset(routineN, handle) output_unit = default_output_unit - eps_warning = eps_check_diag_default -#if defined(__CHECK_DIAG) - check_eigenvectors = .TRUE. - IF (eps_check_diag >= 0.0_dp) THEN - eps_warning = eps_check_diag - END IF -#else - IF (eps_check_diag >= 0.0_dp) THEN - check_eigenvectors = .TRUE. - eps_warning = eps_check_diag - ELSE - check_eigenvectors = .FALSE. - END IF -#endif + check_eigenvectors = diag_check_requested() + eps_warning = diag_check_warning_threshold() eps_abort = 10.0_dp*eps_warning gold = 0.0_dp @@ -398,6 +421,108 @@ CONTAINS END SUBROUTINE check_diag +! ************************************************************************************************** +!> \brief Check C^T*S*C = I for a generalized eigenvalue problem. +!> \param overlap original overlap matrix S; used as work matrix and overwritten +!> \param eigenvectors eigenvectors C to be checked +!> \param scratch work matrix +!> \param nvec ... +! ************************************************************************************************** + SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec) + + TYPE(cp_fm_type), INTENT(IN) :: eigenvectors + TYPE(cp_fm_type), INTENT(INOUT) :: overlap, scratch + INTEGER, INTENT(IN) :: nvec + + CHARACTER(LEN=*), PARAMETER :: routineN = 'check_generalized_diag' + + CHARACTER(LEN=default_string_length) :: diag_type_name + REAL(KIND=dp) :: eps, eps_abort, eps_warning, gold, test + INTEGER :: handle, i, j, ncol, nrow, output_unit +#if defined(__parallel) + TYPE(cp_blacs_env_type), POINTER :: context + INTEGER :: il, jl, ipcol, iprow, & + mypcol, myprow, npcol, nprow + INTEGER, DIMENSION(9) :: desca +#endif + + CALL timeset(routineN, handle) + + IF (.NOT. diag_check_requested()) THEN + CALL timestop(handle) + RETURN + END IF + + output_unit = default_output_unit + eps_warning = diag_check_warning_threshold() + eps_abort = 10.0_dp*eps_warning + + nrow = eigenvectors%matrix_struct%nrow_global + ncol = MIN(eigenvectors%matrix_struct%ncol_global, nvec) + + CALL parallel_gemm("N", "N", nrow, ncol, nrow, 1.0_dp, overlap, eigenvectors, 0.0_dp, scratch) + CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, eigenvectors, scratch, 0.0_dp, overlap) + + gold = 0.0_dp + test = 0.0_dp + eps = 0.0_dp + +#if defined(__parallel) + context => overlap%matrix_struct%context + myprow = context%mepos(1) + mypcol = context%mepos(2) + nprow = context%num_pe(1) + npcol = context%num_pe(2) + desca(:) = overlap%matrix_struct%descriptor(:) + outer: DO j = 1, ncol + DO i = 1, ncol + CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol) + IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN + gold = MERGE(0.0_dp, 1.0_dp, i /= j) + test = overlap%local_data(il, jl) + eps = ABS(test - gold) + IF (eps > eps_warning) EXIT outer + END IF + END DO + END DO outer +#else + outer: DO j = 1, ncol + DO i = 1, ncol + gold = MERGE(0.0_dp, 1.0_dp, i /= j) + test = overlap%local_data(i, j) + eps = ABS(test - gold) + IF (eps > eps_warning) EXIT outer + END DO + END DO outer +#endif + + IF (eps > eps_warning) THEN + IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN + diag_type_name = "SYGVX" + ELSE IF (diag_type == FM_DIAG_TYPE_ELPA) THEN + diag_type_name = "ELPA" + ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER) THEN + diag_type_name = "CUSOLVER" + ELSE IF (diag_type == FM_DIAG_TYPE_DLAF) THEN + diag_type_name = "DLAF" + ELSE + CPABORT("Unknown diag_type") + END IF + WRITE (UNIT=output_unit, FMT="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,F0.0,A,ES10.3)") & + "The generalized eigenvectors returned by "//TRIM(diag_type_name)//" are not S-orthonormal", & + "Matrix element (", i, ", ", j, ") = ", test, & + "The deviation from the expected value ", gold, " is", eps + IF (eps > eps_abort) THEN + CPABORT("ERROR in "//routineN//": Check of generalized matrix diagonalization failed") + ELSE + CPWARN("Check of generalized matrix diagonalization failed in routine "//routineN) + END IF + END IF + + CALL timestop(handle) + + END SUBROUTINE check_generalized_diag + ! ************************************************************************************************** !> \brief Issues an error messages and exits (optionally only warns). !> \param mesg message to be issued @@ -1348,8 +1473,9 @@ CONTAINS END SUBROUTINE cp_fm_block_jacobi ! ************************************************************************************************** -!> \brief General Eigenvalue Problem AX = BXE -!> Single option version: Cholesky decomposition of B +!> \brief General Eigenvalue Problem AX = BXE. +!> Use cuSOLVERMp directly when requested and large enough; otherwise +!> reduce the problem through a Cholesky decomposition of B. !> \param amatrix ... !> \param bmatrix ... !> \param eigenvectors ... @@ -1365,21 +1491,65 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geeig' INTEGER :: handle, nao, nmo + LOGICAL :: check_eigenvectors + TYPE(cp_fm_type) :: overlap_check, scratch_check CALL timeset(routineN, handle) CALL cp_fm_get_info(amatrix, nrow_global=nao) nmo = SIZE(eigenvalues) + check_eigenvectors = diag_check_requested() - IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. cusolver_generalized .AND. nao >= 64) THEN - ! Use cuSolverMP generalized eigenvalue solver for large matrices + IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. direct_generalized_diagonalization .AND. & + nao >= cusolver_n_min) THEN + ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side + ! Cholesky reduction. ! Use work as intermediate buffer since eigenvectors may be smaller (nao x nmo) + IF (check_eigenvectors) THEN + CALL cp_fm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_fm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_fm_to_fm(bmatrix, overlap_check) + END IF CALL cp_fm_general_cusolver(amatrix, bmatrix, work, eigenvalues) + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_fm_release(scratch_check) + CALL cp_fm_release(overlap_check) + END IF CALL cp_fm_to_fm(work, eigenvectors, nmo) +#if defined(__parallel) + ELSE IF (diag_type == FM_DIAG_TYPE_SCALAPACK .AND. direct_generalized_diagonalization) THEN + ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side + ! Cholesky reduction. + IF (check_eigenvectors) THEN + CALL cp_fm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_fm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_fm_to_fm(bmatrix, overlap_check) + END IF + CALL cp_fm_geeig_scalapack(amatrix, bmatrix, work, eigenvalues) + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_fm_release(scratch_check) + CALL cp_fm_release(overlap_check) + END IF + CALL cp_fm_to_fm(work, eigenvectors, nmo) +#endif #if defined(__DLAF) - ELSE IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN + ELSE IF (diag_type == FM_DIAG_TYPE_DLAF .AND. direct_generalized_diagonalization .AND. & + nao >= dlaf_neigvec_min) THEN ! Use DLA-Future generalized eigenvalue solver for large matrices - CALL cp_fm_diag_gen_dlaf(amatrix, bmatrix, eigenvectors, eigenvalues) + IF (check_eigenvectors) THEN + CALL cp_fm_create(overlap_check, bmatrix%matrix_struct) + CALL cp_fm_create(scratch_check, bmatrix%matrix_struct) + CALL cp_fm_to_fm(bmatrix, overlap_check) + END IF + CALL cp_fm_diag_gen_dlaf(amatrix, bmatrix, work, eigenvalues) + IF (check_eigenvectors) THEN + CALL check_generalized_diag(overlap_check, work, scratch_check, nmo) + CALL cp_fm_release(scratch_check) + CALL cp_fm_release(overlap_check) + END IF + CALL cp_fm_to_fm(work, eigenvectors, nmo) #endif ELSE ! Cholesky decompose S=U(T)U @@ -1401,6 +1571,120 @@ CONTAINS END SUBROUTINE cp_fm_geeig +! ************************************************************************************************** +!> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PDSYGVX. +!> \param amatrix ... +!> \param bmatrix ... +!> \param eigenvectors ... +!> \param eigenvalues ... +! ************************************************************************************************** + SUBROUTINE cp_fm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues) + + TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues + + CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geeig_scalapack' + +#if defined(__parallel) + REAL(KIND=dp), PARAMETER :: orfac = -1.0_dp, & + vl = 0.0_dp, & + vu = 0.0_dp + + INTEGER :: handle, info, liwork, lwork, m, n, nb, & + neig, npcol, nprow, nz + INTEGER, DIMENSION(9) :: desca, descb, descz + INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork + REAL(KIND=dp) :: abstol + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: gap, w, work + REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z + + INTEGER :: mq0, nn, np0, npe + INTEGER, EXTERNAL :: iceil, numroc + REAL(KIND=dp), EXTERNAL :: dlamch +#if defined (__HAS_IEEE_EXCEPTIONS) + LOGICAL, DIMENSION(5) :: halt +#endif +#else + INTEGER :: handle +#endif + + CALL timeset(routineN, handle) + +#if defined(__parallel) + n = amatrix%matrix_struct%nrow_global + neig = MIN(SIZE(eigenvalues), n) + + IF (neig == 0) THEN + CALL timestop(handle) + RETURN + END IF + + IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN + CPABORT("ERROR in "//routineN//": Invalid blocksize (no square blocks) found") + END IF + + a => amatrix%local_data + b => bmatrix%local_data + z => eigenvectors%local_data + desca(:) = amatrix%matrix_struct%descriptor(:) + descb(:) = bmatrix%matrix_struct%descriptor(:) + descz(:) = eigenvectors%matrix_struct%descriptor(:) + + nprow = amatrix%matrix_struct%context%num_pe(1) + npcol = amatrix%matrix_struct%context%num_pe(2) + npe = nprow*npcol + nb = amatrix%matrix_struct%nrow_block + nn = MAX(n, nb, 2) + np0 = numroc(nn, nb, 0, 0, nprow) + mq0 = MAX(numroc(nn, nb, 0, 0, npcol), nb) + + lwork = 5*n + MAX(5*nn, np0*mq0 + 2*nb*nb) + iceil(neig, npe)*nn + & + MAX(0, neig - 1)*n + liwork = 6*MAX(n, npe + 1, 4) + + ALLOCATE (gap(npe)) + gap = 0.0_dp + ALLOCATE (iclustr(2*npe)) + iclustr = 0 + ALLOCATE (ifail(n)) + ifail = 0 + ALLOCATE (iwork(liwork)) + ALLOCATE (w(n)) + ALLOCATE (work(lwork)) + + abstol = 2.0_dp*dlamch("S") + +#if defined (__HAS_IEEE_EXCEPTIONS) + CALL ieee_get_halting_mode(IEEE_ALL, halt) + CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.) +#endif + CALL pdsygvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, & + vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, & + work(1), lwork, iwork(1), liwork, ifail(1), iclustr(1), gap(1), info) +#if defined (__HAS_IEEE_EXCEPTIONS) + CALL ieee_set_halting_mode(IEEE_ALL, halt) +#endif + + IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN + CPABORT("ERROR in PDSYGVX (ScaLAPACK), info="//TRIM(cp_to_string(info))) + END IF + + eigenvalues(:) = 0.0_dp + eigenvalues(1:neig) = w(1:neig) + + DEALLOCATE (gap, iclustr, ifail, iwork, w, work) +#else + MARK_USED(amatrix) + MARK_USED(bmatrix) + MARK_USED(eigenvectors) + MARK_USED(eigenvalues) + CPABORT("ERROR in "//routineN//": PDSYGVX requested without ScaLAPACK support") +#endif + + CALL timestop(handle) + + END SUBROUTINE cp_fm_geeig_scalapack + ! ************************************************************************************************** !> \brief General Eigenvalue Problem AX = BXE !> Use canonical diagonalization : U*s**(-1/2) diff --git a/src/fm/cp_fm_dlaf_api.F b/src/fm/cp_fm_dlaf_api.F index 9a5a546ef1..27a64f3723 100644 --- a/src/fm/cp_fm_dlaf_api.F +++ b/src/fm/cp_fm_dlaf_api.F @@ -351,7 +351,7 @@ CONTAINS TYPE(cp_fm_type), INTENT(IN) :: a_matrix, b_matrix, eigenvectors REAL(kind=dp), DIMENSION(:), INTENT(OUT), TARGET :: eigenvalues - CHARACTER(len=*), PARAMETER :: dlaf_name = 'pdsyevd_dlaf', & + CHARACTER(len=*), PARAMETER :: dlaf_name = 'pdsygvd_dlaf', & routineN = 'cp_fm_diag_gen_dlaf_base' CHARACTER, PARAMETER :: uplo = 'L' diff --git a/src/global_types.F b/src/global_types.F index 1a3bad68b0..9d076d7514 100644 --- a/src/global_types.F +++ b/src/global_types.F @@ -85,7 +85,7 @@ MODULE global_types LOGICAL :: elpa_one_stage = .FALSE. ! enable one-stage ELPA solver INTEGER :: dlaf_neigvec_min = 0 ! Minimum number of eigenvectors for DLAF eigensolver usage INTEGER :: dlaf_cholesky_n_min = 0 ! Minimum matrix size for DLAF Cholesky decomposition usage - LOGICAL :: cusolver_generalized = .TRUE. ! use cuSOLVERMp generalized solver + LOGICAL :: direct_generalized_diagonalization = .FALSE. ! use a direct generalized eigensolver LOGICAL :: blacs_repeatable = .FALSE. ! will store the user preference for the repeatability of BLACS collectives REAL(KIND=dp) :: cp2k_start_time = 0.0_dp REAL(KIND=dp) :: cp2k_target_time = HUGE(0.0_dp) ! Maximum run time in seconds diff --git a/src/input_cp2k_global.F b/src/input_cp2k_global.F index 6baf7c17b8..4b67463a76 100644 --- a/src/input_cp2k_global.F +++ b/src/input_cp2k_global.F @@ -141,10 +141,12 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="CUSOLVER_GENERALIZED", & - description="Use cuSOLVERMp to solve the generalized eigenvalue problem directly on the GPU", & - usage="CUSOLVER_GENERALIZED", & - default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) + CALL keyword_create(keyword, __LOCATION__, name="DIRECT_GENERALIZED_DIAGONALIZATION", & + description="Request direct generalized eigenvalue problem diagonalization "// & + "without a CP2K-side Cholesky reduction in supported dense matrix paths. "// & + "The eigensolver is still selected by PREFERRED_DIAG_LIBRARY.", & + usage="DIRECT_GENERALIZED_DIAGONALIZATION", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/qs_scf_diagonalization.F b/src/qs_scf_diagonalization.F index eefc0e8022..1b68fa757c 100644 --- a/src/qs_scf_diagonalization.F +++ b/src/qs_scf_diagonalization.F @@ -46,8 +46,9 @@ MODULE qs_scf_diagonalization choose_eigv_solver,& cp_fm_geeig,& cp_fm_geeig_canon,& - cusolver_generalized,& - diag_type + cusolver_n_min,& + diag_type,& + direct_generalized_diagonalization USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& fm_pool_create_fm,& fm_pool_give_back_fm @@ -168,7 +169,7 @@ CONTAINS TYPE(section_vals_type), POINTER :: scf_section LOGICAL, INTENT(INOUT) :: diis_step - INTEGER :: ispin, nspin + INTEGER :: ispin, nao, nspin LOGICAL :: do_level_shift, owns_ortho, use_jacobi REAL(KIND=dp) :: diis_error, eps_diis TYPE(cp_fm_type), POINTER :: ortho @@ -255,7 +256,9 @@ CONTAINS END IF DO ispin = 1, nspin - IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. cusolver_generalized .AND. .NOT. do_level_shift) THEN + CALL cp_fm_get_info(scf_env%scf_work1(ispin), nrow_global=nao) + IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. direct_generalized_diagonalization .AND. & + nao >= cusolver_n_min .AND. .NOT. do_level_shift) THEN CALL eigensolver_generalized(matrix_ks_fm=scf_env%scf_work1(ispin), & matrix_s=matrix_s(ispin)%matrix, & mo_set=mos(ispin), & diff --git a/src/qs_scf_initialization.F b/src/qs_scf_initialization.F index 0548669cb0..66ea4c8079 100644 --- a/src/qs_scf_initialization.F +++ b/src/qs_scf_initialization.F @@ -29,8 +29,9 @@ MODULE qs_scf_initialization USE cp_fm_diag, ONLY: FM_DIAG_TYPE_CUSOLVER,& choose_eigv_solver,& cp_fm_power,& - cusolver_generalized,& - diag_type + cusolver_n_min,& + diag_type,& + direct_generalized_diagonalization USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& fm_pool_get_el_struct USE cp_fm_struct, ONLY: cp_fm_struct_create,& @@ -805,10 +806,13 @@ CONTAINS scf_env%method = general_diag_method_nr scf_env%needs_ortho = (.NOT. has_unit_metric) .AND. (.NOT. do_kpoints) IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. & - cusolver_generalized .AND. & + direct_generalized_diagonalization .AND. & scf_control%level_shift == 0.0_dp .AND. & scf_env%cholesky_method /= cholesky_off) THEN - scf_env%needs_ortho = .FALSE. + CALL get_mo_set(mos(1), nao=nao) + IF (nao >= cusolver_n_min) THEN + scf_env%needs_ortho = .FALSE. + END IF END IF IF (has_unit_metric) THEN scf_env%method = special_diag_method_nr diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F index 651ae9d00a..5ab71400d6 100644 --- a/src/qs_vxc_atom.F +++ b/src/qs_vxc_atom.F @@ -805,12 +805,12 @@ CONTAINS REAL(KIND=dp) :: agr, alpha, density_cut, gradient_cut, & rtot, tau_cut REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), & - POINTER :: vxc_h, vxc_s + POINTER :: vtau_h, vtau_s, vxc_h, vxc_s REAL(KIND=dp), DIMENSION(1, 1, 1) :: rtau REAL(KIND=dp), DIMENSION(1, 1, 1, 1) :: rrho REAL(KIND=dp), DIMENSION(:, :), POINTER :: weight_h, weight_s REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, & - tau1_s, tau_h, tau_s, vtau_h, vtau_s + tau1_s, tau_h, tau_s REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho1_h, drho1_s, drho_h, drho_s, vxg_h, & vxg_s TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -1572,13 +1572,13 @@ CONTAINS paw_atom, tau_f REAL(dp) :: agr, alpha, beta, density_cut, & gradient_cut, oeps1, tau_cut - REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s + REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s, vxc_h, vxc_s REAL(dp), DIMENSION(1, 1, 1) :: tau_d REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, & rho_s, tau0_h, tau0_s, tau1_h, tau1_s, & - tau_h, tau_s, vtau_h, vtau_s + tau_h, tau_s REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, & drho_h, drho_s, vxg_h, vxg_s REAL(KIND=dp), DIMENSION(-4:4) :: ak diff --git a/src/xc/xc_atom.F b/src/xc/xc_atom.F index bea8930cba..bf569d119c 100644 --- a/src/xc/xc_atom.F +++ b/src/xc/xc_atom.F @@ -414,7 +414,8 @@ CONTAINS REAL(dp), DIMENSION(:, :), INTENT(IN) :: w REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg - REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: vtau + REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), & + OPTIONAL, POINTER :: vtau LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r' diff --git a/tests/QS/regtest-cusolver/Si8-generalized-complex.inp b/tests/QS/regtest-cusolver/Si8-generalized-complex.inp new file mode 100644 index 0000000000..00fb055b40 --- /dev/null +++ b/tests/QS/regtest-cusolver/Si8-generalized-complex.inp @@ -0,0 +1,79 @@ +&GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T + EPS_CHECK_DIAG 1.0E-12 + PREFERRED_DIAG_LIBRARY CUSOLVER + PRINT_LEVEL LOW + PROJECT Si8-generalized-complex + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &KPOINTS + FULL_GRID ON + PARALLEL_GROUP_SIZE 0 + SCHEME MONKHORST-PACK 2 1 1 + SYMMETRY ON + VERBOSE F + WAVEFUNCTIONS COMPLEX + &END KPOINTS + &MGRID + CUTOFF 40 + NGRIDS 4 + &END MGRID + &QS + EXTRAPOLATION PS + EXTRAPOLATION_ORDER 2 + METHOD GPW + &END QS + &SCF + ADDED_MOS 20 20 + EPS_DIIS 1.0E-7 + EPS_SCF 1.0E-5 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 5 + SCF_GUESS MOPAC + &MIXING + ALPHA 0.4 + METHOD DIRECT_P_MIXING + &END MIXING + &OT OFF + ENERGY_GAP 0.001 + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF OFF + EPS_SCF 1.0E-6 + MAX_SCF 1 + &END OUTER_SCF + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &SMEAR + METHOD GAUSSIAN + SIGMA [eV] 0.05 + &END SMEAR + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.42858871335 5.42858871335 5.42858871335 + &END CELL + &KIND Si + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q4 + &END KIND + &TOPOLOGY + CONNECTIVITY OFF + COORDINATE XYZ + COORD_FILE_NAME ../sample_xyz/SI_8.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-cusolver/Si8-generalized.inp b/tests/QS/regtest-cusolver/Si8-generalized.inp new file mode 100644 index 0000000000..eddf09d263 --- /dev/null +++ b/tests/QS/regtest-cusolver/Si8-generalized.inp @@ -0,0 +1,67 @@ +&GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T + EPS_CHECK_DIAG 1.0E-14 + PREFERRED_DIAG_LIBRARY CUSOLVER + PRINT_LEVEL LOW + PROJECT Si8-generalized + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 40 + NGRIDS 4 + &END MGRID + &QS + EXTRAPOLATION PS + EXTRAPOLATION_ORDER 2 + METHOD GPW + &END QS + &SCF + ADDED_MOS 20 20 + EPS_DIIS 1.0E-7 + EPS_SCF 1.0E-5 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 5 + SCF_GUESS MOPAC + &MIXING + ALPHA 0.4 + METHOD DIRECT_P_MIXING + &END MIXING + &OT OFF + ENERGY_GAP 0.001 + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF OFF + EPS_SCF 1.0E-6 + MAX_SCF 1 + &END OUTER_SCF + &SMEAR + METHOD GAUSSIAN + SIGMA [eV] 0.05 + &END SMEAR + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.42858871335 5.42858871335 5.42858871335 + &END CELL + &KIND Si + BASIS_SET DZVP-GTH + POTENTIAL GTH-PBE-q4 + &END KIND + &TOPOLOGY + CONNECTIVITY OFF + COORDINATE XYZ + COORD_FILE_NAME ../sample_xyz/SI_8.xyz + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-cusolver/TEST_FILES.toml b/tests/QS/regtest-cusolver/TEST_FILES.toml index 6dba512161..e5fb2202bd 100644 --- a/tests/QS/regtest-cusolver/TEST_FILES.toml +++ b/tests/QS/regtest-cusolver/TEST_FILES.toml @@ -1,2 +1,4 @@ "H2O-6.inp" = [{matcher="E_total", tol=2e-14, ref=-17.14603641519601}] +"Si8-generalized.inp" = [{matcher="E_total", tol=1e-11, ref=-31.187602969867214}] +"Si8-generalized-complex.inp" = [{matcher="E_total", tol=1e-10, ref=-31.461089087225638}] #EOF diff --git a/tests/QS/regtest-dlaf/H2O-6.inp b/tests/QS/regtest-dlaf/H2O-6.inp index 247ec67684..57e2555019 100644 --- a/tests/QS/regtest-dlaf/H2O-6.inp +++ b/tests/QS/regtest-dlaf/H2O-6.inp @@ -1,4 +1,5 @@ &GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION F DLAF_CHOLESKY_N_MIN 3 DLAF_NEIGVEC_MIN 3 EPS_CHECK_DIAG 1.0E-14 diff --git a/tests/QS/regtest-dlaf/TEST_FILES.toml b/tests/QS/regtest-dlaf/TEST_FILES.toml index 677793404e..9de8ed2054 100644 --- a/tests/QS/regtest-dlaf/TEST_FILES.toml +++ b/tests/QS/regtest-dlaf/TEST_FILES.toml @@ -1,4 +1,5 @@ "H2O-6.inp" = [{matcher="E_total", tol=2e-14, ref=-17.14603641519601}] # Test from regtest-kp-1 "c_2.inp" = [{matcher="E_total", tol=1.0E-14, ref=-45.68042106170509}] +"real_kp.inp" = [{matcher="E_total", tol=1e-10, ref=-46.843539755568322}] #EOF diff --git a/tests/QS/regtest-dlaf/c_2.inp b/tests/QS/regtest-dlaf/c_2.inp index a27f9d7ede..a9c403fbe8 100644 --- a/tests/QS/regtest-dlaf/c_2.inp +++ b/tests/QS/regtest-dlaf/c_2.inp @@ -66,7 +66,9 @@ &END FORCE_EVAL &GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T DLAF_NEIGVEC_MIN 3 + EPS_CHECK_DIAG 1.0E-12 PREFERRED_DIAG_LIBRARY DLAF PRINT_LEVEL LOW PROJECT C diff --git a/tests/QS/regtest-dlaf/real_kp.inp b/tests/QS/regtest-dlaf/real_kp.inp new file mode 100644 index 0000000000..d9480b3c2a --- /dev/null +++ b/tests/QS/regtest-dlaf/real_kp.inp @@ -0,0 +1,88 @@ +@SET NREP 1 +&FORCE_EVAL + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &KPOINTS + EPS_GEO 1.e-8 + FULL_GRID ON + PARALLEL_GROUP_SIZE 0 + SCHEME MONKHORST-PACK 1 1 1 + SYMMETRY ON + VERBOSE F + WAVEFUNCTIONS REAL + &END KPOINTS + &MGRID + CUTOFF 120 + REL_CUTOFF 30 + &END MGRID + &PRINT + &OVERLAP_CONDITION + 1-NORM + DIAGONALIZATION + ARNOLDI + &END OVERLAP_CONDITION + &END PRINT + &QS + EPS_DEFAULT 1.0E-14 + EXTRAPOLATION USE_GUESS + METHOD GPW + &END QS + &SCF + EPS_EIGVAL 1.0E-5 + EPS_SCF 1.0E-6 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS ATOMIC + &MIXING + ALPHA 0.70 + METHOD DIRECT_P_MIXING + &END MIXING + &PRINT + &RESTART off + &END RESTART + &END PRINT + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 3.56683 3.56683 3.56683 + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END CELL + &COORD + SCALED + C 0.000000 0.000000 0.000000 + C 0.500000 0.500000 0.000000 + C 0.500000 0.000000 0.500000 + C 0.000000 0.500000 0.500000 + C 0.250000 0.250000 0.250000 + C 0.250000 0.750000 0.750000 + C 0.750000 0.250000 0.750000 + C 3/4 3/4 1/4 + &END COORD + &KIND C + BASIS_SET TZV2P-GTH + POTENTIAL GTH-PADE-q4 + &END KIND + &TOPOLOGY + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + +&GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T + DLAF_NEIGVEC_MIN 0 + EPS_CHECK_DIAG 1.0E-14 + PREFERRED_DIAG_LIBRARY DLAF + PRINT_LEVEL LOW + PROJECT C + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.0 + &END TIMINGS +&END GLOBAL diff --git a/tests/QS/regtest-scalapack-generalized/TEST_FILES.toml b/tests/QS/regtest-scalapack-generalized/TEST_FILES.toml new file mode 100644 index 0000000000..6761b68f2c --- /dev/null +++ b/tests/QS/regtest-scalapack-generalized/TEST_FILES.toml @@ -0,0 +1,2 @@ +"real_kp.inp" = [{matcher="E_total", tol=1e-10, ref=-46.843539755568315}] +"complex_kp.inp" = [{matcher="E_total", tol=1e-10, ref=-41.897080363575512}] diff --git a/tests/QS/regtest-scalapack-generalized/complex_kp.inp b/tests/QS/regtest-scalapack-generalized/complex_kp.inp new file mode 100644 index 0000000000..a6831678bf --- /dev/null +++ b/tests/QS/regtest-scalapack-generalized/complex_kp.inp @@ -0,0 +1,71 @@ +&GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T + EPS_CHECK_DIAG 1.0E-14 + PREFERRED_DIAG_LIBRARY ScaLAPACK + PRINT_LEVEL medium + PROJECT NO2 + RUN_TYPE energy +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_SET + LSD + POTENTIAL_FILE_NAME POTENTIAL + &KPOINTS + FULL_GRID ON + PARALLEL_GROUP_SIZE 0 + SCHEME MONKHORST-PACK 2 2 2 + SYMMETRY ON + VERBOSE F + WAVEFUNCTIONS COMPLEX + &END KPOINTS + &MGRID + CUTOFF 160 + &END MGRID + &PRINT + &LOWDIN ON + PRINT_ALL + &END LOWDIN + &MULLIKEN OFF + &END MULLIKEN + &END PRINT + &QS + EPS_DEFAULT 1.0E-10 + &END QS + &SCF + EPS_DIIS 0.1 + EPS_SCF 1.0E-6 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS ATOMIC + &MIXING + ALPHA 0.2 + METHOD DIRECT_P_MIXING + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL Pade + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + O 3.71590 2.5 2.5 + O 1.65536 3.37465 2.5 + N 2.5 2.5 2.5 + &END COORD + &KIND O + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q6 + &END KIND + &KIND N + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-scalapack-generalized/real_kp.inp b/tests/QS/regtest-scalapack-generalized/real_kp.inp new file mode 100644 index 0000000000..b49c70b072 --- /dev/null +++ b/tests/QS/regtest-scalapack-generalized/real_kp.inp @@ -0,0 +1,84 @@ +@SET NREP 1 +&FORCE_EVAL + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &KPOINTS + EPS_GEO 1.e-8 + FULL_GRID ON + PARALLEL_GROUP_SIZE 0 + SCHEME MONKHORST-PACK 1 1 1 + SYMMETRY ON + VERBOSE F + WAVEFUNCTIONS REAL + &END KPOINTS + &MGRID + CUTOFF 120 + REL_CUTOFF 30 + &END MGRID + &PRINT + &OVERLAP_CONDITION + 1-NORM + DIAGONALIZATION + ARNOLDI + &END OVERLAP_CONDITION + &END PRINT + &QS + EPS_DEFAULT 1.0E-14 + EXTRAPOLATION USE_GUESS + METHOD GPW + &END QS + &SCF + EPS_EIGVAL 1.0E-5 + EPS_SCF 1.0E-6 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 1 + SCF_GUESS ATOMIC + &MIXING + ALPHA 0.70 + METHOD DIRECT_P_MIXING + &END MIXING + &PRINT + &RESTART off + &END RESTART + &END PRINT + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 3.56683 3.56683 3.56683 + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END CELL + &COORD + SCALED + C 0.000000 0.000000 0.000000 + C 0.500000 0.500000 0.000000 + C 0.500000 0.000000 0.500000 + C 0.000000 0.500000 0.500000 + C 0.250000 0.250000 0.250000 + C 0.250000 0.750000 0.750000 + C 0.750000 0.250000 0.750000 + C 3/4 3/4 1/4 + &END COORD + &KIND C + BASIS_SET TZV2P-GTH + POTENTIAL GTH-PADE-q4 + &END KIND + &TOPOLOGY + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + +&GLOBAL + DIRECT_GENERALIZED_DIAGONALIZATION T + EPS_CHECK_DIAG 1.0E-14 + PREFERRED_DIAG_LIBRARY ScaLAPACK + PRINT_LEVEL LOW + PROJECT C + RUN_TYPE ENERGY +&END GLOBAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 4a7586115e..fe086ff234 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -21,6 +21,7 @@ QS/regtest-double-hybrid-stress-laplace libint libxc !ifx QS/regtest-rpa-sigma libint greenx !ifx TMC/regtest_ana_on_the_fly parallel mpiranks>2 QS/regtest-cusolver cusolvermp +QS/regtest-scalapack-generalized parallel scalapack QS/regtest-dlaf dlaf QS/regtest-dlaf-2 dlaf libint !ifx QS/regtest-pao-2 diff --git a/tools/toolchain/install_cp2k_toolchain.sh b/tools/toolchain/install_cp2k_toolchain.sh index 4bdce9e8d3..ce2b6188de 100755 --- a/tools/toolchain/install_cp2k_toolchain.sh +++ b/tools/toolchain/install_cp2k_toolchain.sh @@ -141,7 +141,7 @@ OPTIONS: Default = native --gpu-ver Select the target GPU architecture for compiling. Available options are: K20X, K40, K80, P100, V100, - A100, H100, A40, Mi50, Mi100, Mi250, and no. + A100, H100, GB10, A40, Mi50, Mi100, Mi250, and no. This option determines the value of nvcc -arch flag. Default = no --libint-lmax Maximum supported angular momentum by libint if the @@ -640,7 +640,8 @@ while [ $# -ge 1 ]; do for ii in ${package_list}; do if [ "${ii}" != "intel" ] && [ "${ii}" != "intelmpi" ] && - [ "${ii}" != "amd" ]; then + [ "${ii}" != "amd" ] && + [ "${ii}" != "cusolvermp" ]; then eval "with_${ii}=__INSTALL__" fi done @@ -697,13 +698,13 @@ Otherwise use option no." --gpu-ver=*) user_input="${1#*=}" case "${user_input}" in - K20X | K40 | K80 | P100 | V100 | A100 | H100 | A40 | Mi50 | Mi100 | Mi250 | no) + K20X | K40 | K80 | P100 | V100 | A100 | H100 | GB10 | A40 | Mi50 | Mi100 | Mi250 | no) export GPUVER="${user_input}" ;; *) report_error ${LINENO} "Invalid value for --gpu-ver found." echo "Currently only one of the following options is supported: - K20X, K40, K80, P100, V100, A100, H100, A40, Mi50, Mi100, Mi250. + K20X, K40, K80, P100, V100, A100, H100, GB10, A40, Mi50, Mi100, Mi250. Otherwise use option no." exit 1 ;; @@ -845,7 +846,7 @@ Otherwise use option no." with_elpa=$(read_with "${1}") ;; --with-cusolvermp*) - with_cusolvermp=$(read_with "${1}") + with_cusolvermp=$(read_with "${1}" "__SYSTEM__") ;; --with-deepmd*) with_deepmd=$(read_with "${1}") @@ -1166,6 +1167,9 @@ case ${GPUVER} in H100) export ARCH_NUM="90" ;; + GB10) + export ARCH_NUM="121" + ;; Mi50) # TODO: export ARCH_NUM= ;; diff --git a/tools/toolchain/scripts/stage4/install_cusolvermp.sh b/tools/toolchain/scripts/stage4/install_cusolvermp.sh new file mode 100755 index 0000000000..9944522f33 --- /dev/null +++ b/tools/toolchain/scripts/stage4/install_cusolvermp.sh @@ -0,0 +1,143 @@ +#!/bin/bash -e + +# TODO: Review and if possible fix shellcheck errors. +# shellcheck disable=all + +[ "${BASH_SOURCE[0]}" ] && SCRIPT_NAME="${BASH_SOURCE[0]}" || SCRIPT_NAME=$0 +SCRIPT_DIR="$(cd "$(dirname "$SCRIPT_NAME")/.." && pwd -P)" + +source "${SCRIPT_DIR}"/common_vars.sh +source "${SCRIPT_DIR}"/tool_kit.sh +source "${SCRIPT_DIR}"/signal_trap.sh +source "${INSTALLDIR}"/toolchain.conf +source "${INSTALLDIR}"/toolchain.env + +[ -f "${BUILDDIR}/setup_cusolvermp" ] && rm "${BUILDDIR}/setup_cusolvermp" + +! [ -d "${BUILDDIR}" ] && mkdir -p "${BUILDDIR}" +cd "${BUILDDIR}" + +case "${with_cusolvermp}" in + __INSTALL__) + report_error ${LINENO} "Installing cuSOLVERMp is not supported. Use --with-cusolvermp=system or --with-cusolvermp=." + exit 1 + ;; + __SYSTEM__) + echo "==================== Finding cuSOLVERMp from system paths ====================" + add_include_from_paths CUSOLVERMP_CFLAGS "cusolverMp.h" $INCLUDE_PATHS + add_lib_from_paths CUSOLVERMP_LDFLAGS "libcusolverMp.*" $LIB_PATHS + cusolvermp_header="$(find_in_paths "cusolverMp.h" $INCLUDE_PATHS)" + cusolvermp_lib="$(find_in_paths "libcusolverMp.*" $LIB_PATHS)" + ;; + __DONTUSE__) ;; + *) + echo "==================== Linking cuSOLVERMp to user paths ====================" + pkg_install_dir="${with_cusolvermp}" + CUSOLVERMP_LIBDIR="${pkg_install_dir}/lib" + [ -d "${pkg_install_dir}/lib64" ] && CUSOLVERMP_LIBDIR="${pkg_install_dir}/lib64" + + check_dir "${CUSOLVERMP_LIBDIR}" + check_dir "${pkg_install_dir}/include" + CUSOLVERMP_CFLAGS="-I'${pkg_install_dir}/include'" + CUSOLVERMP_LDFLAGS="-L'${CUSOLVERMP_LIBDIR}' -Wl,-rpath,'${CUSOLVERMP_LIBDIR}'" + cusolvermp_header="${pkg_install_dir}/include/cusolverMp.h" + cusolvermp_lib="$(find_in_paths "libcusolverMp.*" CUSOLVERMP_LIBDIR)" + ;; +esac + +if [ "${with_cusolvermp}" != "__DONTUSE__" ]; then + if [ "${cusolvermp_header}" = "__FALSE__" ] || ! [ -f "${cusolvermp_header}" ]; then + report_error ${LINENO} "Could not find cusolverMp.h." + exit 1 + fi + if [ "${cusolvermp_lib}" = "__FALSE__" ] || ! [ -f "${cusolvermp_lib}" ]; then + report_error ${LINENO} "Could not find libcusolverMp." + exit 1 + fi + + pkg_install_dir="$(dirname "$(dirname "${cusolvermp_lib}")")" + cusolvermp_major="$(grep -E '^#define CUSOLVERMP_VER_MAJOR' "${cusolvermp_header}" | awk '{print $3}')" + cusolvermp_minor="$(grep -E '^#define CUSOLVERMP_VER_MINOR' "${cusolvermp_header}" | awk '{print $3}')" + cusolvermp_major="${cusolvermp_major:-0}" + cusolvermp_minor="${cusolvermp_minor:-0}" + + CUSOLVERMP_LIBS="-lcusolverMp" + CUSOLVERMP_DFLAGS="-D__CUSOLVERMP" + if [ "${cusolvermp_major}" -gt 0 ] || [ "${cusolvermp_minor}" -ge 7 ]; then + add_include_from_paths NCCL_CFLAGS "nccl.h" $INCLUDE_PATHS + add_lib_from_paths NCCL_LDFLAGS "libnccl.*" $LIB_PATHS + nccl_lib="$(find_in_paths "libnccl.*" $LIB_PATHS)" + if [ "${nccl_lib}" = "__FALSE__" ]; then + report_error ${LINENO} "Could not find NCCL required by cuSOLVERMp ${cusolvermp_major}.${cusolvermp_minor}." + exit 1 + fi + NCCL_ROOT="$(dirname "$(dirname "${nccl_lib}")")" + CUSOLVERMP_CFLAGS="${CUSOLVERMP_CFLAGS} ${NCCL_CFLAGS}" + CUSOLVERMP_LIBS="${CUSOLVERMP_LIBS} -lnccl" + CUSOLVERMP_LDFLAGS="${CUSOLVERMP_LDFLAGS} ${NCCL_LDFLAGS}" + CUSOLVERMP_DFLAGS="${CUSOLVERMP_DFLAGS} -D__CUSOLVERMP_NCCL" + else + add_include_from_paths CAL_CFLAGS "cal.h" $INCLUDE_PATHS + add_lib_from_paths CAL_LDFLAGS "libcal.*" $LIB_PATHS + cal_lib="$(find_in_paths "libcal.*" $LIB_PATHS)" + if [ "${cal_lib}" = "__FALSE__" ]; then + report_error ${LINENO} "Could not find CAL required by cuSOLVERMp ${cusolvermp_major}.${cusolvermp_minor}." + exit 1 + fi + CAL_ROOT="$(dirname "$(dirname "${cal_lib}")")" + CUSOLVERMP_CFLAGS="${CUSOLVERMP_CFLAGS} ${CAL_CFLAGS}" + CUSOLVERMP_LIBS="${CUSOLVERMP_LIBS} -lcal" + CUSOLVERMP_LDFLAGS="${CUSOLVERMP_LDFLAGS} ${CAL_LDFLAGS}" + fi + + add_include_from_paths UCC_CFLAGS "ucc/api/ucc.h" $INCLUDE_PATHS + add_lib_from_paths UCC_LDFLAGS "libucc.*" $LIB_PATHS + add_lib_from_paths UCX_LDFLAGS "libucs.*" $LIB_PATHS + ucc_lib="$(find_in_paths "libucc.*" $LIB_PATHS)" + ucx_lib="$(find_in_paths "libucs.*" $LIB_PATHS)" + if [ "${ucc_lib}" = "__FALSE__" ] || [ "${ucx_lib}" = "__FALSE__" ]; then + report_error ${LINENO} "Could not find UCC/UCX required by cuSOLVERMp." + exit 1 + fi + UCC_ROOT="$(dirname "$(dirname "${ucc_lib}")")" + UCX_ROOT="$(dirname "$(dirname "${ucx_lib}")")" + CUSOLVERMP_LDFLAGS="${CUSOLVERMP_LDFLAGS} ${UCC_LDFLAGS} ${UCX_LDFLAGS}" + CUSOLVERMP_LIBS="${CUSOLVERMP_LIBS} -lucc -lucs" + + cat << EOF > "${BUILDDIR}/setup_cusolvermp" +export CUSOLVERMP_VER="${cusolvermp_major}.${cusolvermp_minor}" +export CUSOLVERMP_CFLAGS="${CUSOLVERMP_CFLAGS}" +export CUSOLVERMP_LDFLAGS="${CUSOLVERMP_LDFLAGS}" +export CUSOLVERMP_LIBS="${CUSOLVERMP_LIBS}" +export CUSOLVER_MP_ROOT="${pkg_install_dir}" +export UCC_ROOT="${UCC_ROOT}" +export UCX_ROOT="${UCX_ROOT}" +export CP_DFLAGS="\${CP_DFLAGS} IF_CUDA(${CUSOLVERMP_DFLAGS}|)" +export CP_CFLAGS="\${CP_CFLAGS} IF_CUDA(${CUSOLVERMP_CFLAGS}|)" +export CP_LDFLAGS="\${CP_LDFLAGS} IF_CUDA(${CUSOLVERMP_LDFLAGS}|)" +export CP_LIBS="IF_CUDA(${CUSOLVERMP_LIBS}|) \${CP_LIBS}" +prepend_path CMAKE_PREFIX_PATH "${pkg_install_dir}" +prepend_path CMAKE_PREFIX_PATH "${UCC_ROOT}" +prepend_path CMAKE_PREFIX_PATH "${UCX_ROOT}" +EOF + if [ -n "${NCCL_ROOT}" ]; then + cat << EOF >> "${BUILDDIR}/setup_cusolvermp" +export NCCL_ROOT="${NCCL_ROOT}" +prepend_path CMAKE_PREFIX_PATH "${NCCL_ROOT}" +EOF + fi + if [ -n "${CAL_ROOT}" ]; then + cat << EOF >> "${BUILDDIR}/setup_cusolvermp" +export CAL_ROOT="${CAL_ROOT}" +prepend_path CMAKE_PREFIX_PATH "${CAL_ROOT}" +EOF + fi + filter_setup "${BUILDDIR}/setup_cusolvermp" "${SETUPFILE}" +fi + +load "${BUILDDIR}/setup_cusolvermp" +write_toolchain_env "${INSTALLDIR}" + +report_timing "cusolvermp" + +#EOF diff --git a/tools/toolchain/scripts/stage4/install_stage4.sh b/tools/toolchain/scripts/stage4/install_stage4.sh index b73e6462de..ea1a82fe88 100755 --- a/tools/toolchain/scripts/stage4/install_stage4.sh +++ b/tools/toolchain/scripts/stage4/install_stage4.sh @@ -5,6 +5,7 @@ ./scripts/stage4/install_libxsmm.sh ./scripts/stage4/install_scalapack.sh +./scripts/stage4/install_cusolvermp.sh ./scripts/stage4/install_cosma.sh #EOF diff --git a/tools/toolchain/scripts/stage9/install_dbcsr.sh b/tools/toolchain/scripts/stage9/install_dbcsr.sh index d941a61900..d412c8842c 100755 --- a/tools/toolchain/scripts/stage9/install_dbcsr.sh +++ b/tools/toolchain/scripts/stage9/install_dbcsr.sh @@ -32,6 +32,14 @@ case "${with_dbcsr}" in [ -d dbcsr-${dbcsr_ver} ] && rm -rf dbcsr-${dbcsr_ver} tar -xzf dbcsr-${dbcsr_ver}.tar.gz cd dbcsr-${dbcsr_ver} + if [ "${ENABLE_CUDA}" == "__TRUE__" ] && [ "${GPUVER}" == "GB10" ]; then + # DBCSR 2.9.1 predates GB10. Build native sm_121 code while + # reusing the closest available libsmm_acc parameters. + sed -i "s/ H100)/ H100\\n GB10)/" CMakeLists.txt + sed -i "/ set(GPU_ARCH_NUMBER_H100 90)/a\\ set(GPU_ARCH_NUMBER_GB10 121)" CMakeLists.txt + cp src/acc/libsmm_acc/parameters/parameters_H100.json \ + src/acc/libsmm_acc/parameters/parameters_GB10.json + fi mkdir build-cpu cd build-cpu CMAKE_OPTIONS="-DBUILD_TESTING=NO -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_VERBOSE_MAKEFILE=ON" @@ -56,7 +64,14 @@ case "${with_dbcsr}" in if [ "${ENABLE_CUDA}" == "__TRUE__" ]; then mkdir build-cuda cd build-cuda - CMAKE_OPTIONS="${CMAKE_OPTIONS} -DUSE_ACCEL=cuda -DWITH_GPU=P100" + CMAKE_OPTIONS="${CMAKE_OPTIONS} -DUSE_ACCEL=cuda" + # CUDA 13 deprecates APIs still used by DBCSR 2.9.1 under -Werror. + CMAKE_OPTIONS="${CMAKE_OPTIONS} -DCMAKE_CXX_FLAGS=-Wno-error=deprecated-declarations" + if [ "${GPUVER}" == "GB10" ]; then + CMAKE_OPTIONS="${CMAKE_OPTIONS} -DWITH_GPU=GB10 -DWITH_GPU_PARAMS=GB10" + else + CMAKE_OPTIONS="${CMAKE_OPTIONS} -DWITH_GPU=P100" + fi cmake \ -DCMAKE_INSTALL_PREFIX=${pkg_install_dir}-cuda \ ${CMAKE_OPTIONS} .. \