DLA-Future complex eigensolvers

This commit is contained in:
Rocco Meli 2024-12-17 15:09:03 +01:00 committed by GitHub
parent 5b22e8b240
commit 8512fd5998
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
24 changed files with 768 additions and 210 deletions

2
.gitignore vendored
View file

@ -99,6 +99,7 @@ build*/
# Spack
build-*
.spack_patched
### CMake Patch ###
# External projects
@ -313,6 +314,7 @@ venv/
ENV/
env.bak/
venv.bak/
.envrc
# Spyder project settings
.spyderproject

View file

@ -244,16 +244,22 @@ tools for the solution of dense linear systems and eigenvalue problems.
- Add `-D__CUSOLVERMP` to `DFLAGS`
- Add `-lcusolverMp -lcusolver -lcal -lnvidia-ml` to `LIBS`
### 2m. DLA-Future (optional, experimental, improved performance for diagonalization on Nvidia and AMD GPUs)
### 2m. DLA-Future (optional, improved performance for diagonalization on Nvidia and AMD GPUs)
[DLA-Future](https://github.com/eth-cscs/DLA-Future) is a high-performance, distributed-memory,
GPU-accelerated library that provides tools for the solution of eigenvalue problems, based on the
[pika](https://pikacpp.org/) runtime.
[DLA-Future-Fortran](https://github.com/eth-cscs/DLA-Future-Fortran) provides a Fortran interface to
DLA-Future.
- DLA-Future replaces the ScaLAPACK `SYEVD` to improve performance of the diagonalization
- DLA-Future-Fortran replaces the ScaLAPACK `SYEVD`, `HEEVD`, and `HEGVD` to improve performance of
the diagonalization
- DLA-Future is available at <https://github.com/eth-cscs/DLA-Future>
- DLA-Future-Fortran is available at <https://github.com/eth-cscs/DLA-Future-Fortran>
- DLA-Future is available via the [Spack](https://packages.spack.io/package.html?name=dla-future)
package manager
- DLA-Future-Fortran is available via the
[Spack](https://packages.spack.io/package.html?name=dla-future-fortran) package manager
- `-D__DLAF` is defined by CMake when `-DCP2K_USE_DLAF=ON`
### 2n. PEXSI (optional, low scaling SCF method)

View file

@ -1015,6 +1015,7 @@ list(
CP2K_SRCS_F
fm/cp_blacs_env.F
fm/cp_blacs_types.F
fm/cp_cfm_cholesky.F
fm/cp_cfm_basic_linalg.F
fm/cp_cfm_diag.F
fm/cp_cfm_basic_linalg.F
@ -1032,7 +1033,8 @@ list(
list(APPEND CP2K_SRCS_C fm/cp_fm_cusolver.c)
list(APPEND CP2K_SRCS_F fm/cp_dlaf_utils_api.F fm/cp_fm_dlaf_api.F)
list(APPEND CP2K_SRCS_F fm/cp_dlaf_utils_api.F fm/cp_fm_dlaf_api.F
fm/cp_cfm_dlaf_api.F)
list(APPEND CP2K_SRCS_F hfxbase/hfx_compression_core_methods.F
hfxbase/hfx_contract_block.F hfxbase/hfx_contraction_methods.F)

View file

@ -19,12 +19,12 @@ MODULE admm_methods
USE atomic_kind_types, ONLY: atomic_kind_type
USE bibliography, ONLY: Merlot2014,&
cite_reference
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_cholesky_invert,&
cp_cfm_scale,&
USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
cp_cfm_scale_and_add,&
cp_cfm_scale_and_add_fm,&
cp_cfm_transpose
USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_cholesky_invert
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_get_info,&
cp_cfm_release,&

View file

@ -14,8 +14,8 @@ MODULE rt_propagation_methods
Kuhne2007,&
cite_reference
USE cell_types, ONLY: cell_type
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_triangular_multiply
USE cp_cfm_basic_linalg, ONLY: cp_cfm_triangular_multiply
USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_release,&
cp_cfm_type

View file

@ -36,8 +36,7 @@ MODULE cp_cfm_basic_linalg
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_basic_linalg'
PUBLIC :: cp_cfm_cholesky_decompose, &
cp_cfm_column_scale, &
PUBLIC :: cp_cfm_column_scale, &
cp_cfm_gemm, &
cp_cfm_lu_decompose, &
cp_cfm_lu_invert, &
@ -53,8 +52,8 @@ MODULE cp_cfm_basic_linalg
cp_cfm_triangular_multiply, &
cp_cfm_rot_rows, &
cp_cfm_rot_cols, &
cp_cfm_cholesky_invert, &
cp_cfm_det ! determinant of a complex matrix with correct sign
cp_cfm_det, & ! determinant of a complex matrix with correct sign
cp_cfm_upper_to_full
REAL(kind=dp), EXTERNAL :: zlange, pzlange
@ -906,118 +905,6 @@ CONTAINS
END SUBROUTINE cp_cfm_lu_invert
! **************************************************************************************************
!> \brief Used to replace a symmetric positive definite matrix M with its Cholesky
!> decomposition U: M = U^T * U, with U upper triangular.
!> \param matrix the matrix to replace with its Cholesky decomposition
!> \param n the number of row (and columns) of the matrix &
!> (defaults to the min(size(matrix)))
!> \param info_out if present, outputs info from (p)zpotrf
!> \par History
!> 05.2002 created [JVdV]
!> 12.2002 updated, added n optional parm [fawzi]
!> 09.2021 removed CPASSERT(info == 0) since there is already check of info [Jan Wilhelm]
!> \author Joost
! **************************************************************************************************
SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out)
TYPE(cp_cfm_type), INTENT(IN) :: matrix
INTEGER, INTENT(in), OPTIONAL :: n
INTEGER, INTENT(out), OPTIONAL :: info_out
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_decompose'
COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
INTEGER :: handle, info, my_n
#if defined(__parallel)
INTEGER, DIMENSION(9) :: desca
#else
INTEGER :: lda
#endif
CALL timeset(routineN, handle)
my_n = MIN(matrix%matrix_struct%nrow_global, &
matrix%matrix_struct%ncol_global)
IF (PRESENT(n)) THEN
CPASSERT(n <= my_n)
my_n = n
END IF
a => matrix%local_data
#if defined(__parallel)
desca(:) = matrix%matrix_struct%descriptor(:)
CALL pzpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
#else
lda = SIZE(a, 1)
CALL zpotrf('U', my_n, a(1, 1), lda, info)
#endif
IF (PRESENT(info_out)) THEN
info_out = info
ELSE
IF (info /= 0) &
CALL cp_abort(__LOCATION__, &
"Cholesky decompose failed: matrix is not positive definite or ill-conditioned")
END IF
CALL timestop(handle)
END SUBROUTINE cp_cfm_cholesky_decompose
! **************************************************************************************************
!> \brief Used to replace Cholesky decomposition by the inverse.
!> \param matrix : the matrix to invert (must be an upper triangular matrix),
!> and is the output of Cholesky decomposition
!> \param n : size of the matrix to invert (defaults to the min(size(matrix)))
!> \param info_out : if present, outputs info of (p)zpotri
!> \par History
!> 05.2002 created Lianheng Tong, based on cp_fm_cholesky_invert
!> \author Lianheng Tong
! **************************************************************************************************
SUBROUTINE cp_cfm_cholesky_invert(matrix, n, info_out)
TYPE(cp_cfm_type), INTENT(IN) :: matrix
INTEGER, INTENT(in), OPTIONAL :: n
INTEGER, INTENT(out), OPTIONAL :: info_out
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_invert'
COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
INTEGER :: info, handle
INTEGER :: my_n
#if defined(__parallel)
INTEGER, DIMENSION(9) :: desca
#endif
CALL timeset(routineN, handle)
my_n = MIN(matrix%matrix_struct%nrow_global, &
matrix%matrix_struct%ncol_global)
IF (PRESENT(n)) THEN
CPASSERT(n <= my_n)
my_n = n
END IF
aa => matrix%local_data
#if defined(__parallel)
desca = matrix%matrix_struct%descriptor
CALL pzpotri('U', my_n, aa(1, 1), 1, 1, desca, info)
#else
CALL zpotri('U', my_n, aa(1, 1), SIZE(aa, 1), info)
#endif
IF (PRESENT(info_out)) THEN
info_out = info
ELSE
IF (info /= 0) &
CALL cp_abort(__LOCATION__, &
"Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
END IF
CALL timestop(handle)
END SUBROUTINE cp_cfm_cholesky_invert
! **************************************************************************************************
!> \brief Returns the trace of matrix_a^T matrix_b, i.e
!> sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
@ -1457,4 +1344,64 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE cp_cfm_rot_cols
! **************************************************************************************************
!> \brief ...
!> \param matrix ...
!> \param workspace ...
!> \par History
!> 12.2024 Added optional workspace as input [Rocco Meli]
!> \author Jan Wilhelm
! **************************************************************************************************
SUBROUTINE cp_cfm_upper_to_full(matrix, workspace)
TYPE(cp_cfm_type), INTENT(IN) :: matrix
TYPE(cp_cfm_type), OPTIONAL :: workspace
CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_upper_to_full'
INTEGER :: handle, i_global, iiB, j_global, jjB, &
ncol_local, nrow_local
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
TYPE(cp_cfm_type) :: work
CALL timeset(routineN, handle)
IF (.NOT. PRESENT(workspace)) THEN
CALL cp_cfm_create(work, matrix%matrix_struct)
ELSE
work = workspace
END IF
! get info of fm_mat_Q
CALL cp_cfm_get_info(matrix=matrix, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)
IF (j_global < i_global) THEN
matrix%local_data(iiB, jjB) = z_zero
END IF
IF (j_global == i_global) THEN
matrix%local_data(iiB, jjB) = matrix%local_data(iiB, jjB)/(2.0_dp, 0.0_dp)
END IF
END DO
END DO
CALL cp_cfm_transpose(matrix, 'C', work)
CALL cp_cfm_scale_and_add(z_one, matrix, z_one, work)
IF (.NOT. PRESENT(workspace)) THEN
CALL cp_cfm_release(work)
END IF
CALL timestop(handle)
END SUBROUTINE cp_cfm_upper_to_full
END MODULE cp_cfm_basic_linalg

154
src/fm/cp_cfm_cholesky.F Normal file
View file

@ -0,0 +1,154 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief various cholesky decomposition related routines
!> \par History
!> 12.2002 Moved routines from cp_cfm_basic_linalg to this new module [Rocco Meli]
! **************************************************************************************************
MODULE cp_cfm_cholesky
USE cp_cfm_types, ONLY: cp_cfm_type
USE kinds, ONLY: dp
#if defined(__DLAF)
USE cp_cfm_dlaf_api, ONLY: cp_cfm_pzpotrf_dlaf
#endif
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_cholesky'
PUBLIC :: cp_cfm_cholesky_decompose, &
cp_cfm_cholesky_invert
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
!> \brief Used to replace a symmetric positive definite matrix M with its Cholesky
!> decomposition U: M = U^T * U, with U upper triangular.
!> \param matrix the matrix to replace with its Cholesky decomposition
!> \param n the number of row (and columns) of the matrix &
!> (defaults to the min(size(matrix)))
!> \param info_out if present, outputs info from (p)zpotrf
!> \par History
!> 05.2002 created [JVdV]
!> 12.2002 updated, added n optional parm [fawzi]
!> 09.2021 removed CPASSERT(info == 0) since there is already check of info [Jan Wilhelm]
!> 12.2024 Added DLA-Future support [Rocco Meli]
!> \author Joost
! **************************************************************************************************
SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out)
TYPE(cp_cfm_type), INTENT(IN) :: matrix
INTEGER, INTENT(in), OPTIONAL :: n
INTEGER, INTENT(out), OPTIONAL :: info_out
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_decompose'
COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a
INTEGER :: handle, info, my_n
#if defined(__parallel)
INTEGER, DIMENSION(9) :: desca
#else
INTEGER :: lda
#endif
CALL timeset(routineN, handle)
my_n = MIN(matrix%matrix_struct%nrow_global, &
matrix%matrix_struct%ncol_global)
IF (PRESENT(n)) THEN
CPASSERT(n <= my_n)
my_n = n
END IF
a => matrix%local_data
#if defined(__parallel)
desca(:) = matrix%matrix_struct%descriptor(:)
#if defined(__DLAF)
CALL cp_cfm_pzpotrf_dlaf('U', my_n, a, 1, 1, desca, info)
#else
CALL pzpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
#endif
#else
lda = SIZE(a, 1)
CALL zpotrf('U', my_n, a(1, 1), lda, info)
#endif
IF (PRESENT(info_out)) THEN
info_out = info
ELSE
IF (info /= 0) &
CALL cp_abort(__LOCATION__, &
"Cholesky decompose failed: matrix is not positive definite or ill-conditioned")
END IF
CALL timestop(handle)
END SUBROUTINE cp_cfm_cholesky_decompose
! **************************************************************************************************
!> \brief Used to replace Cholesky decomposition by the inverse.
!> \param matrix : the matrix to invert (must be an upper triangular matrix),
!> and is the output of Cholesky decomposition
!> \param n : size of the matrix to invert (defaults to the min(size(matrix)))
!> \param info_out : if present, outputs info of (p)zpotri
!> \par History
!> 05.2002 created Lianheng Tong, based on cp_fm_cholesky_invert
!> \author Lianheng Tong
! **************************************************************************************************
SUBROUTINE cp_cfm_cholesky_invert(matrix, n, info_out)
TYPE(cp_cfm_type), INTENT(IN) :: matrix
INTEGER, INTENT(in), OPTIONAL :: n
INTEGER, INTENT(out), OPTIONAL :: info_out
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_invert'
COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
INTEGER :: info, handle
INTEGER :: my_n
#if defined(__parallel)
INTEGER, DIMENSION(9) :: desca
#endif
CALL timeset(routineN, handle)
my_n = MIN(matrix%matrix_struct%nrow_global, &
matrix%matrix_struct%ncol_global)
IF (PRESENT(n)) THEN
CPASSERT(n <= my_n)
my_n = n
END IF
aa => matrix%local_data
#if defined(__parallel)
desca = matrix%matrix_struct%descriptor
CALL pzpotri('U', my_n, aa(1, 1), 1, 1, desca, info)
#else
CALL zpotri('U', my_n, aa(1, 1), SIZE(aa, 1), info)
#endif
IF (PRESENT(info_out)) THEN
info_out = info
ELSE
IF (info /= 0) &
CALL cp_abort(__LOCATION__, &
"Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
END IF
CALL timestop(handle)
END SUBROUTINE cp_cfm_cholesky_invert
END MODULE cp_cfm_cholesky

View file

@ -12,8 +12,8 @@
!> \author Joost VandeVondele (2003-09)
! **************************************************************************************************
MODULE cp_cfm_diag
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose, &
cp_cfm_gemm, &
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, &
@ -22,6 +22,11 @@ MODULE cp_cfm_diag
cp_cfm_set_element, &
cp_cfm_to_cfm, &
cp_cfm_type
#if defined(__DLAF)
USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
cp_cfm_diag_dlaf
USE cp_fm_diag, ONLY: diag_type, dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
#endif
USE kinds, ONLY: dp
#if defined (__HAS_IEEE_EXCEPTIONS)
USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
@ -46,15 +51,49 @@ CONTAINS
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \par History
!> - (De)Allocation checks updated (15.02.2011,MK)
!> 12.2024 Added DLA-Future support [Rocco Meli]
!> \author Joost VandeVondele
! **************************************************************************************************
SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
INTEGER :: handle
CALL timeset(routineN, handle)
#if defined(__DLAF)
IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
ELSE
#endif
CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
#if defined(__DLAF)
END IF
#endif
CALL timestop(handle)
END SUBROUTINE cp_cfm_heevd
! **************************************************************************************************
!> \brief Perform a diagonalisation of a complex matrix
!> \param matrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \par History
!> - (De)Allocation checks updated (15.02.2011,MK)
!> \author Joost VandeVondele
! **************************************************************************************************
SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
COMPLEX(KIND=dp), DIMENSION(:, :), &
@ -131,7 +170,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE cp_cfm_heevd
END SUBROUTINE cp_cfm_heevd_base
! **************************************************************************************************
!> \brief General Eigenvalue Problem AX = BXE
@ -141,6 +180,8 @@ CONTAINS
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \param work ...
!> \par History
!> 12.2024 Added DLA-Future support [Rocco Meli]
! **************************************************************************************************
SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
@ -157,18 +198,29 @@ CONTAINS
CALL cp_cfm_get_info(amatrix, nrow_global=nao)
ALLOCATE (evals(nao))
! Cholesky decompose S=U(T)U
CALL cp_cfm_cholesky_decompose(bmatrix)
! Invert to get U^(-1)
CALL cp_cfm_triangular_invert(bmatrix)
! Reduce to get U^(-T) * H * U^(-1)
CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
! Diagonalize
CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
! Restore vectors C = U^(-1) * C*
CALL cp_cfm_triangular_multiply(bmatrix, work)
nmo = SIZE(eigenvalues)
#if defined(__DLAF)
IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
! Use DLA-Future generalized eigenvalue solver for large matrices
CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
ELSE
#endif
! Cholesky decompose S=U(T)U
CALL cp_cfm_cholesky_decompose(bmatrix)
! Invert to get U^(-1)
CALL cp_cfm_triangular_invert(bmatrix)
! Reduce to get U^(-T) * H * U^(-1)
CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
! Diagonalize
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)

287
src/fm/cp_cfm_dlaf_api.F Normal file
View file

@ -0,0 +1,287 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
MODULE cp_cfm_dlaf_api
USE cp_cfm_basic_linalg, ONLY: cp_cfm_upper_to_full
USE cp_cfm_types, ONLY: cp_cfm_type
#if defined(__DLAF)
USE dlaf_fortran, ONLY: dlaf_pzheevd, &
dlaf_pzhegvd, &
dlaf_pzpotrf
#endif
USE kinds, ONLY: dp
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_dlaf_api'
PUBLIC :: cp_cfm_pzpotrf_dlaf
PUBLIC :: cp_cfm_diag_dlaf, cp_cfm_diag_gen_dlaf
CONTAINS
!***************************************************************************************************
!> \brief Cholesky factorization using DLA-Future
!> \param uplo ...
!> \param n Matrix size
!> \param a Local matrix
!> \param ia Row index of first row (has to be 1)
!> \param ja Col index of first column (has to be 1)
!> \param desca ScaLAPACK matrix descriptor
!> \param info 0 if factorization completed normally
!> \author Rocco Meli
! **************************************************************************************************
SUBROUTINE cp_cfm_pzpotrf_dlaf(uplo, n, a, ia, ja, desca, info)
CHARACTER, INTENT(IN) :: uplo
INTEGER, INTENT(IN) :: n
COMPLEX(KIND=dp), DIMENSION(:, :), TARGET :: a
INTEGER, INTENT(IN) :: ia, ja
INTEGER, DIMENSION(9) :: desca
INTEGER, TARGET :: info
CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_pzpotrf_dlaf'
INTEGER :: handle
CALL timeset(routineN, handle)
#if defined(__DLAF)
CALL dlaf_pzpotrf(uplo, n, a, ia, ja, desca, info)
#else
MARK_USED(uplo)
MARK_USED(n)
MARK_USED(a)
MARK_USED(ia)
MARK_USED(ja)
MARK_USED(desca)
MARK_USED(info)
CPABORT("CP2K compiled without the DLA-Future library.")
#endif
CALL timestop(handle)
END SUBROUTINE cp_cfm_pzpotrf_dlaf
! **************************************************************************************************
!> \brief DLA-Future eigensolver for complex Hermitian matrices
!> \param matrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \author Rocco Meli
! **************************************************************************************************
SUBROUTINE cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_diag_dlaf'
INTEGER :: handle, n, nmo
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: eig
CALL timeset(routineN, handle)
n = matrix%matrix_struct%nrow_global
ALLOCATE (eig(n))
CALL cp_cfm_diag_dlaf_base(matrix, eigenvectors, eig)
nmo = SIZE(eigenvalues, 1)
IF (nmo > n) THEN
eigenvalues(1:n) = eig(1:n)
ELSE
eigenvalues(1:nmo) = eig(1:nmo)
END IF
DEALLOCATE (eig)
CALL timestop(handle)
END SUBROUTINE cp_cfm_diag_dlaf
! **************************************************************************************************
!> \brief DLA-Future generalized eigensolver for complex Hermitian matrices
!> \param amatrix ...
!> \param bmatrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \author Rocco Meli
! **************************************************************************************************
SUBROUTINE cp_cfm_diag_gen_dlaf(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_diag_gen_dlaf'
INTEGER :: handle, n, nmo
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: eig
CALL timeset(routineN, handle)
n = amatrix%matrix_struct%nrow_global
ALLOCATE (eig(n))
CALL cp_cfm_diag_gen_dlaf_base(amatrix, bmatrix, eigenvectors, eig)
nmo = SIZE(eigenvalues, 1)
IF (nmo > n) THEN
eigenvalues(1:n) = eig(1:n)
ELSE
eigenvalues(1:nmo) = eig(1:nmo)
END IF
DEALLOCATE (eig)
CALL timestop(handle)
END SUBROUTINE cp_cfm_diag_gen_dlaf
!***************************************************************************************************
!> \brief DLA-Future standard eigensolver for complex Hermitian matrices
!> \param matrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \author Rocco Meli
! **************************************************************************************************
SUBROUTINE cp_cfm_diag_dlaf_base(matrix, eigenvectors, eigenvalues)
TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
REAL(kind=dp), DIMENSION(:), INTENT(OUT), TARGET :: eigenvalues
CHARACTER(len=*), PARAMETER :: dlaf_name = 'pzheevd_dlaf', routineN = 'cp_cfm_diag_dlaf_base'
CHARACTER, PARAMETER :: uplo = 'L'
CHARACTER(LEN=100) :: message
COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, z
INTEGER :: dlaf_handle, handle, n
INTEGER, DIMENSION(9) :: desca, descz
INTEGER, TARGET :: info
CALL timeset(routineN, handle)
#if defined(__DLAF)
! DLAF needs the lower triangular part
! Use eigenvectors matrix as workspace
CALL cp_cfm_upper_to_full(matrix, eigenvectors)
n = matrix%matrix_struct%nrow_global
a => matrix%local_data
z => eigenvectors%local_data
desca(:) = matrix%matrix_struct%descriptor(:)
descz(:) = eigenvectors%matrix_struct%descriptor(:)
info = -1
CALL timeset(dlaf_name, dlaf_handle)
CALL dlaf_pzheevd(uplo, n, a, 1, 1, desca, eigenvalues, z, 1, 1, descz, info)
CALL timestop(dlaf_handle)
IF (info /= 0) THEN
WRITE (message, "(A,I0,A)") "ERROR in DLAF_PZHEEVD: Eigensolver failed (INFO = ", info, ")"
CPABORT(TRIM(message))
END IF
#else
MARK_USED(a)
MARK_USED(z)
MARK_USED(desca)
MARK_USED(descz)
MARK_USED(matrix)
MARK_USED(eigenvectors)
MARK_USED(eigenvalues)
MARK_USED(uplo)
MARK_USED(n)
MARK_USED(info)
MARK_USED(dlaf_handle)
MARK_USED(dlaf_name)
MARK_USED(message)
CPABORT("CP2K compiled without DLAF library.")
#endif
CALL timestop(handle)
END SUBROUTINE cp_cfm_diag_dlaf_base
!***************************************************************************************************
!> \brief DLA-Future generalized eigensolver for complex Hermitian matrices
!> \param amatrix ...
!> \param bmatrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \author Rocco Meli
! **************************************************************************************************
SUBROUTINE cp_cfm_diag_gen_dlaf_base(amatrix, bmatrix, eigenvectors, eigenvalues)
TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
REAL(kind=dp), DIMENSION(:), INTENT(OUT), TARGET :: eigenvalues
CHARACTER(len=*), PARAMETER :: dlaf_name = 'pzhegvd_dlaf', &
routineN = 'cp_cfm_diag_gen_dlaf_base'
CHARACTER, PARAMETER :: uplo = 'L'
CHARACTER(LEN=100) :: message
COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
INTEGER :: dlaf_handle, handle, n
INTEGER, DIMENSION(9) :: desca, descb, descz
INTEGER, TARGET :: info
CALL timeset(routineN, handle)
#if defined(__DLAF)
! DLAF needs the lower triangular part
! Use eigenvectors matrix as workspace
CALL cp_cfm_upper_to_full(amatrix, eigenvectors)
CALL cp_cfm_upper_to_full(bmatrix, eigenvectors)
n = amatrix%matrix_struct%nrow_global
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(:)
info = -1
CALL timeset(dlaf_name, dlaf_handle)
CALL dlaf_pzhegvd(uplo, n, a, 1, 1, desca, b, 1, 1, descb, eigenvalues, z, 1, 1, descz, info)
CALL timestop(dlaf_handle)
IF (info /= 0) THEN
WRITE (message, "(A,I0,A)") "ERROR in DLAF_PZHEGVD: Eigensolver failed (INFO = ", info, ")"
CPABORT(TRIM(message))
END IF
#else
MARK_USED(a)
MARK_USED(b)
MARK_USED(z)
MARK_USED(desca)
MARK_USED(descb)
MARK_USED(descz)
MARK_USED(amatrix)
MARK_USED(bmatrix)
MARK_USED(eigenvectors)
MARK_USED(eigenvalues)
MARK_USED(uplo)
MARK_USED(n)
MARK_USED(info)
MARK_USED(dlaf_handle)
MARK_USED(dlaf_name)
MARK_USED(message)
CPABORT("CP2K compiled without DLAF library.")
#endif
CALL timestop(handle)
END SUBROUTINE cp_cfm_diag_gen_dlaf_base
END MODULE cp_cfm_dlaf_api

View file

@ -87,14 +87,14 @@ MODULE cp_fm_diag
! The following saved variables are diagonalization global
! Stores the default library for diagonalization
INTEGER, SAVE :: diag_type = 0
INTEGER, SAVE, PUBLIC :: diag_type = 0
! 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
#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 :: dlaf_neigvec_min = 0
INTEGER, SAVE, PUBLIC :: dlaf_neigvec_min = 0
#endif
! Threshold value for the orthonormality check of the eigenvectors obtained
! after a diagonalization. A negative value disables the check.

View file

@ -16,7 +16,8 @@ MODULE gw_large_cell_gamma
get_cell,&
pbc
USE constants_operator, ONLY: operator_coulomb
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
USE cp_cfm_basic_linalg, ONLY: cp_cfm_upper_to_full
USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_cholesky_invert
USE cp_cfm_diag, ONLY: cp_cfm_geeig
USE cp_cfm_types, ONLY: cp_cfm_create,&
@ -80,8 +81,7 @@ MODULE gw_large_cell_gamma
qs_environment_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_tensors, ONLY: build_3c_integrals
USE rpa_gw_kpoints_util, ONLY: cp_cfm_power,&
cp_cfm_upper_to_full
USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
#include "./base/base_uses.f90"
IMPLICIT NONE

View file

@ -18,8 +18,8 @@ MODULE post_scf_bandstructure_utils
get_cell,&
pbc
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_scale
USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale
USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
cp_cfm_geeig_canon,&
cp_cfm_heevd

View file

@ -15,12 +15,11 @@ MODULE rpa_gw_kpoints_util
get_cell,&
pbc
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_cholesky_invert,&
cp_cfm_column_scale,&
cp_cfm_scale_and_add,&
USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
cp_cfm_scale_and_add_fm,&
cp_cfm_transpose
cp_cfm_upper_to_full
USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose,&
cp_cfm_cholesky_invert
USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
cp_cfm_geeig_canon,&
cp_cfm_heevd
@ -88,7 +87,7 @@ MODULE rpa_gw_kpoints_util
PUBLIC :: invert_eps_compute_W_and_Erpa_kp, cp_cfm_power, real_space_to_kpoint_transform_rpa, &
get_mat_cell_T_from_mat_gamma, get_bandstruc_and_k_dependent_MOs, &
compute_wkp_W, mat_kp_from_mat_gamma, cp_cfm_upper_to_full
compute_wkp_W, mat_kp_from_mat_gamma
CONTAINS
@ -1422,55 +1421,6 @@ CONTAINS
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param cfm_mat_Q ...
! **************************************************************************************************
SUBROUTINE cp_cfm_upper_to_full(cfm_mat_Q)
TYPE(cp_cfm_type), INTENT(IN) :: cfm_mat_Q
CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_upper_to_full'
INTEGER :: handle, i_global, iiB, j_global, jjB, &
ncol_local, nrow_local
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
TYPE(cp_cfm_type) :: cfm_mat_work
CALL timeset(routineN, handle)
CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct)
! get info of fm_mat_Q
CALL cp_cfm_get_info(matrix=cfm_mat_Q, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)
DO jjB = 1, ncol_local
j_global = col_indices(jjB)
DO iiB = 1, nrow_local
i_global = row_indices(iiB)
IF (j_global < i_global) THEN
cfm_mat_Q%local_data(iiB, jjB) = z_zero
END IF
IF (j_global == i_global) THEN
cfm_mat_Q%local_data(iiB, jjB) = cfm_mat_Q%local_data(iiB, jjB)/(2.0_dp, 0.0_dp)
END IF
END DO
END DO
CALL cp_cfm_transpose(cfm_mat_Q, 'C', cfm_mat_work)
CALL cp_cfm_scale_and_add(z_one, cfm_mat_Q, z_one, cfm_mat_work)
CALL cp_cfm_release(cfm_mat_work)
CALL timestop(handle)
END SUBROUTINE cp_cfm_upper_to_full
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...

View file

@ -0,0 +1,67 @@
&GLOBAL
DLAF_NEIGVEC_MIN 3
PREFERRED_DIAG_LIBRARY DLAF
PRINT_LEVEL LOW
PROJECT Ar-HF-2p-SOC-os
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.0
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
AUTO_BASIS RI_XAS MEDIUM
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
LSD
POTENTIAL_FILE_NAME POTENTIAL
&POISSON
PERIODIC NONE
PSOLVER MT
&END POISSON
&QS
METHOD GAPW
&END QS
&XAS_TDP
EXCITATIONS OS_SPIN_CONS
EXCITATIONS OS_SPIN_FLIP
GRID Ar 150 150
N_EXCITED 4
SPIN_ORBIT_COUPLING
TAMM_DANCOFF
&DONOR_STATES
AT_LIST 1
DEFINE_EXCITED BY_INDEX
N_SEARCH 5
STATE_TYPES 2p
&END DONOR_STATES
&KERNEL
&EXACT_EXCHANGE
OPERATOR COULOMB
&END EXACT_EXCHANGE
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&END KERNEL
&END XAS_TDP
&XC
&HF
&END HF
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
Ar 2.5 2.5 2.5
&END COORD
&KIND Ar
BASIS_SET 6-31G*
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL

View file

@ -1,2 +1,4 @@
# Tests from regtest-xastdp
Ne-pbc-shortrange.inp 88 1e-08 890.608044
Ar-HF-2p-SOC-os.inp 88 1e-08 261.347679
#EOF

View file

@ -4,6 +4,9 @@
PREFERRED_DIAG_LIBRARY DLAF
PRINT_LEVEL low
PROJECT H2O-6
&TIMINGS
THRESHOLD 0.0
&END TIMINGS
&END GLOBAL
&FORCE_EVAL

View file

@ -1,2 +1,4 @@
H2O-6.inp 1 2e-14 -17.14603641519601
# Test from regtest-kp-1
c_2.inp 1 1.0E-14 -45.68042106170509
#EOF

View file

@ -0,0 +1,73 @@
@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 4 4 4
SYMMETRY ON
VERBOSE F
&END KPOINTS
&MGRID
CUTOFF 120
REL_CUTOFF 30
&END MGRID
&QS
EPS_DEFAULT 1.0E-12
EXTRAPOLATION USE_GUESS
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-6
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 5
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 0.750000 0.750000 0.250000
&END COORD
&KIND C
BASIS_SET DZVP-GTH
POTENTIAL GTH-PADE-q4
&END KIND
&TOPOLOGY
MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
DLAF_NEIGVEC_MIN 3
PREFERRED_DIAG_LIBRARY DLAF
PRINT_LEVEL LOW
PROJECT C
&END GLOBAL

View file

@ -306,6 +306,9 @@ class Config:
env["CUDA_VISIBLE_DEVICES"] = ",".join(visible_gpu_devices)
env["HIP_VISIBLE_DEVICES"] = ",".join(visible_gpu_devices)
env["OMP_NUM_THREADS"] = str(self.ompthreads)
env["PIKA_COMMANDLINE_OPTIONS"] = (
f"--pika:bind=none --pika:threads={self.ompthreads}"
)
exe_name = f"{exe_stem}.{self.version}"
cmd = [str(self.binary_dir / exe_name)]
if self.valgrind:

View file

@ -53,6 +53,13 @@ ENV PATH="/opt/spack/bin:${PATH}"
RUN spack compiler find
RUN spack external find --all --not-buildable
# Enable Spack build cache
ARG SPACK_BUILD_CACHE=develop-2024-12-15
RUN spack mirror add ${SPACK_BUILD_CACHE} https://binaries.spack.io/${SPACK_BUILD_CACHE} && \
spack mirror add develop https://binaries.spack.io/develop && \
spack buildcache keys --install --trust --force && \
spack mirror rm develop
# Install CP2K's dependencies via Spack.
WORKDIR /
COPY ./tools/spack/cp2k-dependencies.yaml .

View file

@ -851,6 +851,13 @@ ENV PATH="/opt/spack/bin:${{PATH}}"
RUN spack compiler find
RUN spack external find --all --not-buildable
# Enable Spack build cache
ARG SPACK_BUILD_CACHE=develop-2024-12-15
RUN spack mirror add ${{SPACK_BUILD_CACHE}} https://binaries.spack.io/${{SPACK_BUILD_CACHE}} && \
spack mirror add develop https://binaries.spack.io/develop && \
spack buildcache keys --install --trust --force && \
spack mirror rm develop
# Install CP2K's dependencies via Spack.
WORKDIR /
COPY ./tools/spack/cp2k-dependencies.yaml .

View file

@ -38,7 +38,7 @@ if [[ "${PROFILE}" == "spack" ]] && [[ "${VERSION}" == "psmp" ]]; then
-DCP2K_USE_DLAF=ON \
-DCP2K_BLAS_VENDOR=OpenBLAS \
-DCP2K_USE_SPGLIB=ON \
-DCP2K_USE_LIBINT2=OFF \
-DCP2K_USE_LIBINT2=ON \
-DCP2K_USE_LIBXC=ON \
-DCP2K_USE_LIBTORCH=OFF \
-DCP2K_USE_MPI=ON \

View file

@ -33,9 +33,6 @@ export OMP_STACKSIZE=64m
# Improve code coverage on COSMA.
export COSMA_DIM_THRESHOLD=0
# Bind pika threads to first two cores. This is a hack. Do not use for production!
export PIKA_PROCESS_MASK="0x3"
# Load Spack or Toolchain environment.
if [[ "${PROFILE}" == "spack" ]]; then
eval "$(spack env activate myenv --sh)"

View file

@ -1,5 +1,4 @@
# This is a Spack Environment file for developing CP2K.
spack:
specs:
- mpich@4.1.2 device=ch3 netmod=tcp
@ -11,12 +10,12 @@ spack:
- "libxc @6.2.2"
- "spglib @2.3.0"
- "dla-future-fortran@0.2.0"
- "libint@2.9.0 +fortran tune=cp2k-lmax-5"
# Unfortunately, ScaLAPACK 2.2.1 has not yet been packaged by Spack.
# https://github.com/Reference-ScaLAPACK/scalapack/tree/v2.2.1
# which contains https://github.com/Reference-ScaLAPACK/scalapack/pull/26
# TODO:
# - "libint @2.6.0 +fortran tune=cp2k-lmax-5" # build takes 16 min
# - "sirius @7.4.3 +fortran +shared +openmp" # build fails
# - "plumed @2.9.0 +shared +mpi optional_modules=all" # requires cython
# - "libvori @220621" # is not recompile with -fPIC
@ -31,7 +30,6 @@ spack:
# libvdwxc-0.4.0.tar.gz: OK
# libtorch-cxx11-abi-shared-with-deps-1.12.1+cpu.zip: OK
# SpLA-1.5.5.tar.gz: OK
view:
default:
root: ./spack-env/view
@ -40,5 +38,4 @@ spack:
- gcc-runtime
concretizer:
unify: true
#EOF