Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2024 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : module m_nvlamath
7 : use m_types_solver
8 : use m_types_mat
9 : use m_judft
10 : #ifdef CPP_GPU_NVLAMATH
11 : use nvlamath
12 : #endif
13 : implicit none
14 : private
15 : type, extends(t_solver)::t_solver_nvlamath
16 : contains
17 : !procedure :: solve_std_dp => nvlamath_diag !solver for standard eigenvalue problem
18 : !procedure :: solve_std_sp => nvlamath_diag_sp !solver for standard eigenvalue problem
19 : procedure :: solve_gev => nvlamath_gev !solver for generalized eigenvalue problem
20 : !procedure :: to_std => nvlamath_reduction !transform the H of the generalized problem to a std problem
21 : !procedure :: backtrans => nvlamath_recover !transform the Eigenvalue back to the generalized problem
22 : end type
23 : public :: t_solver_nvlamath, get_solver_nvlamath
24 :
25 : contains
26 :
27 97 : function get_solver_nvlamath() result(solver)
28 : type(t_solver_nvlamath), pointer::solver
29 97 : allocate (solver)
30 97 : solver%name = "nvlamath"
31 : #ifdef CPP_GPU_NVLAMATH
32 : solver%available = .true.
33 : #else
34 97 : solver%available = .false.
35 : #endif
36 97 : solver%serial = .true.
37 97 : solver%generalized = .true.
38 97 : solver%GPU = .true.
39 97 : end function
40 :
41 0 : subroutine nvlamath_gev(self, hmat, smat, ne, eig, zmat, ikpt)
42 : !Simple driver to solve Generalized Eigenvalue Problem using nvlamath routine
43 : implicit none
44 : class(t_solver_nvlamath) :: self
45 : class(t_mat), intent(INOUT) :: hmat, smat
46 : integer, intent(INOUT) :: ne
47 : class(t_mat), allocatable, intent(OUT) :: zmat
48 : real, intent(OUT) :: eig(:)
49 : integer, intent(IN) :: ikpt
50 :
51 : integer :: lwork, info, m
52 0 : integer, allocatable:: ifail(:), iwork(:)
53 0 : complex, allocatable:: work(:)
54 0 : real, allocatable :: rwork(:)
55 : real :: dumrwork(1), abstol
56 : complex :: dumwork(1)
57 : real, external :: dlamch
58 0 : real :: eigTemp(hmat%matsize1)
59 :
60 0 : allocate (t_mat::zmat)
61 0 : call zmat%alloc(hmat%l_real, hmat%matsize1, ne)
62 0 : abstol = 2*dlamch('S')
63 0 : if (hmat%l_real) then
64 0 : allocate (iwork(5*hmat%matsize1), ifail(hmat%matsize1))
65 : !$acc data create(iwork,ifail,eigTemp,zmat%data_r)
66 : !$acc host_data use_device(hmat%data_r,smat%data_r,eigTemp,zmat%data_r,dumrwork,iwork,ifail)
67 : call dsygvx(1, 'V', 'I', 'U', hmat%matsize1, hmat%data_r, size(hmat%data_r, 1), smat%data_r, size(smat%data_r, 1), &
68 0 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_r, size(zmat%data_r, 1), dumrwork, -1, iwork, ifail, info)
69 : !$acc end host_data
70 : !$acc update self(dumrwork)
71 0 : lwork = dumrwork(1)
72 0 : allocate (rwork(lwork))
73 : !$acc data create(rwork)
74 0 : if (info .ne. 0) call judft_error("Diagonalization via nvlamath failed (Workspace)", no=info)
75 : !$acc host_data use_device(hmat%data_r,smat%data_r,eigTemp,zmat%data_r,rwork,iwork,ifail)
76 : call dsygvx(1, 'V', 'I', 'U', hmat%matsize1, hmat%data_r, size(hmat%data_r, 1), smat%data_r, size(smat%data_r, 1), &
77 0 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_r, size(zmat%data_r, 1), rwork, lwork, iwork, ifail, info)
78 : !$acc end host_data
79 : !$acc end data
80 : !$acc update self(eigTemp,zmat%data_r)
81 : !$acc end data
82 : else
83 :
84 0 : allocate (rwork(7*hmat%matsize1), iwork(5*hmat%matsize1), ifail(hmat%matsize1))
85 : !$acc data create(iwork,ifail,eigTemp,zmat%data_c,rwork)
86 : !Do a workspace query
87 : !$acc host_data use_device(hmat%data_c,smat%data_c,eigTemp,zmat%data_c,dumwork,rwork,iwork,ifail)
88 :
89 : call zhegvx(1, 'V', 'I', 'U', hmat%matsize1, hmat%data_c, size(hmat%data_c, 1), smat%data_c, size(smat%data_c, 1), &
90 0 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_c, size(zmat%data_c, 1), dumwork, -1, rwork, iwork, ifail, info)
91 : !$acc end host_data
92 : !$acc update self(dumwork)
93 0 : lwork = dumwork(1)
94 0 : allocate (work(lwork))
95 : !$acc data create(work)
96 0 : if (info .ne. 0) call judft_error("Diagonalization via nvlamath failed (Workspace)", no=info)
97 : !Perform diagonalization
98 : !$acc host_data use_device(hmat%data_c,smat%data_c,eigTemp,zmat%data_c,work,rwork,iwork,ifail)
99 : call zhegvx(1, 'V', 'I', 'U', hmat%matsize1, hmat%data_c, size(hmat%data_c, 1), smat%data_c, size(smat%data_c, 1), &
100 0 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_c, size(zmat%data_c, 1), work, lwork, rwork, iwork, ifail, info)
101 : !$acc end host_data
102 : !$acc end data
103 : !$acc update self(eigTemp,zmat%data_c)
104 : !$acc end data
105 : end if
106 0 : eig(:min(size(eig), size(eigTemp))) = eigTemp(:min(size(eig), size(eigTemp)))
107 0 : if (info .ne. 0) call judft_error("Diagonalization via nvlamath failed(zhegvx/dsygvx)", no=info)
108 0 : if (m .ne. ne) call judft_error("Diagonalization via nvlamath failed failed without explicit errorcode.")
109 0 : end subroutine nvlamath_gev
110 :
111 97 : end module m_nvlamath
|
