comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 21:55:16 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Diagonalization: add opt-in direct generalized eigensolvers (#5108)

Browse source 
Co-authored-by: DCM-Uni-Paderborn <DCM-Uni-Paderborn@users.noreply.github.com>
Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-05 16:26:45 +02:00 • committed by GitHub
parent 4a77a3b84e
commit c7270bbdb5
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
34 changed files with 1667 additions and 85 deletions
Download patch file Download diff file Expand all files Collapse all files

69
CMakeLists.txt
View file

@ -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)

15
cmake/CompilerConfiguration.cmake
View file

@ -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(
"$<$<COMPILE_LANG_AND_ID:Fortran,PGI>:-Mfreeform;-Mextend;-Mallocatable=03>"
"$<$<COMPILE_LANG_AND_ID:Fortran,PGI,NVHPC>:-Mfreeform;-Mextend;-Mallocatable=03>"
"$<$<COMPILE_LANG_AND_ID:Fortran,NAG>:-f2008;-free;-Warn=reallocation;-Warn=subnormal>"
"$<$<COMPILE_LANG_AND_ID:Fortran,Cray>:-f;free;-M3105;-ME7212;-hnoacc;-M1234>"
)
@ -166,11 +167,11 @@ add_compile_options(
# Release
add_compile_options(
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:Fortran,PGI>>:-fast>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:Fortran,PGI,NVHPC>>:-fast>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:Fortran,Cray>>:-O2;-G2>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:Fortran,NAG>>:-gline>")
add_compile_options(
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:C,PGI>>:-fast>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:C,PGI,NVHPC>>:-fast>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:C,Intel>>:-O3;-g>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:C,Cray>>:-O3>"
"$<$<AND:$<CONFIG:RELEASE>,$<COMPILE_LANG_AND_ID:C,NAG>>:-gline>"
@ -181,10 +182,10 @@ add_compile_options(
# Debug
add_compile_options(
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:Fortran,PGI>>:-g>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:Fortran,PGI,NVHPC>>:-g>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:Fortran,Cray>>:-G2>")
add_compile_options(
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:C,PGI>>:-fast>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:C,PGI,NVHPC>>:-fast>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:C,Intel>>:-O2;-g>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:C,Cray>>:-G2>"
"$<$<AND:$<CONFIG:DEBUG>,$<COMPILE_LANG_AND_ID:C,NAG>>:-g;-C>"

35
cmake/modules/FindSCALAPACK.cmake
View file

@ -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(

1
cmake/modules/Finducc.cmake
View file

@ -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")

8
cmake/modules/cp2k_utils.cmake
View file

@ -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

10
docs/technologies/eigensolvers/cusolvermp.md
View file

@ -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

6
src/CMakeLists.txt
View file

@ -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}
$<$<BOOL:${CUDAToolkit_FOUND}>:${CUDAToolkit_INCLUDE_DIRS}>
@ -1752,7 +1756,7 @@ target_link_libraries(
$<$<BOOL:${CP2K_USE_MIMIC}>:CP2K::MIMIC::mcl>
$<$<BOOL:${CP2K_USE_LIBINT2}>:cp2k::Libint2::int2>
$<$<BOOL:${CP2K_USE_COSMA}>:cp2k::cosma>
$<$<BOOL:${CP2K_USE_DLAF}>:DLAF::Fortran>
$<$<BOOL:${CP2K_USE_DLAF}>:${CP2K_DLAF_FORTRAN_TARGET}>
$<$<BOOL:${CP2K_USE_GREENX}>:cp2k::greenx>
$<$<BOOL:${CP2K_USE_TREXIO}>:cp2k::trexio::trexio>
$<$<BOOL:${CP2K_USE_HDF5}>:HDF5::HDF5

20
src/environment.F
View file

@ -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)// &

291
src/fm/cp_cfm_diag.F
View file

@ -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

254
src/fm/cp_fm_cusolver.c
View file

@ -9,10 +9,12 @@
#include "../offload/offload_library.h"
#include <assert.h>
#include <cuComplex.h>
#include <cuda_runtime.h>
#include <cusolverMp.h>
#include <math.h>
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
@ -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

82
src/fm/cp_fm_cusolver_api.F
View file

@ -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

340
src/fm/cp_fm_diag.F
View file

@ -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)

2
src/fm/cp_fm_dlaf_api.F
View file

@ -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'

2
src/global_types.F
View file

@ -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

10
src/input_cp2k_global.F
View file

@ -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)

11
src/qs_scf_diagonalization.F
View file

@ -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), &

12
src/qs_scf_initialization.F
View file

@ -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

8
src/qs_vxc_atom.F
View file

@ -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

3
src/xc/xc_atom.F
View file

@ -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'

79
tests/QS/regtest-cusolver/Si8-generalized-complex.inp Normal file
View file

@ -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

67
tests/QS/regtest-cusolver/Si8-generalized.inp Normal file
View file

@ -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

2
tests/QS/regtest-cusolver/TEST_FILES.toml
View file

@ -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

1
tests/QS/regtest-dlaf/H2O-6.inp
View file

@ -1,4 +1,5 @@
&GLOBAL
DIRECT_GENERALIZED_DIAGONALIZATION F
DLAF_CHOLESKY_N_MIN 3
DLAF_NEIGVEC_MIN 3
EPS_CHECK_DIAG 1.0E-14

1
tests/QS/regtest-dlaf/TEST_FILES.toml
View file

@ -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

2
tests/QS/regtest-dlaf/c_2.inp
View file

@ -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

88
tests/QS/regtest-dlaf/real_kp.inp Normal file
View file

@ -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

2
tests/QS/regtest-scalapack-generalized/TEST_FILES.toml Normal file
View file

@ -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}]

71
tests/QS/regtest-scalapack-generalized/complex_kp.inp Normal file
View file

@ -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

84
tests/QS/regtest-scalapack-generalized/real_kp.inp Normal file
View file

@ -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

1
tests/TEST_DIRS
View file

@ -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

14
tools/toolchain/install_cp2k_toolchain.sh
View file

@ -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=
;;

143
tools/toolchain/scripts/stage4/install_cusolvermp.sh Executable file
View file

@ -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=<prefix>."
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

1
tools/toolchain/scripts/stage4/install_stage4.sh
View file

@ -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

17
tools/toolchain/scripts/stage9/install_dbcsr.sh
View file

@ -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} .. \