Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 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 :
7 : module m_magma
8 : !! Implementation of a t_solver using the MAGMA library from ICL
9 : use m_juDFT
10 : use m_types_solver
11 : use m_types_mat
12 : #ifdef CPP_MAGMA
13 : use magma
14 : use openacc
15 : #endif
16 : private
17 : logical ,save :: initialized=.false.
18 : !integer, save :: Magma_NumGPU = 1
19 : type, extends(t_solver)::t_solver_magma
20 : !! provides all solvers& transforms for a "serial" case on the GPU
21 : contains
22 : procedure :: solve_gev => magma_GEV
23 : procedure :: solve_std_dp => magma_diag !solver for standard eigenvalue problem
24 : procedure :: solve_std_sp => magma_diag_sp !solver for standard eigenvalue problem
25 : procedure :: to_std => magma_reduction !transform the H of the generalized problem to a std problem
26 : procedure :: backtrans => magma_recover !transform the Eigenvalue back to the generalized problem
27 : end type
28 : public :: get_solver_magma
29 :
30 : contains
31 :
32 91 : function get_solver_magma() result(solver)
33 : type(t_solver_magma), pointer::solver
34 91 : allocate (solver)
35 91 : solver%name = "magma"
36 : #ifdef CPP_MAGMA
37 : solver%available = .true.
38 : #else
39 91 : solver%available = .false.
40 : #endif
41 91 : solver%parallel = .false.
42 91 : solver%serial = .true.
43 91 : solver%generalized = .true.
44 91 : solver%standard = .true.
45 91 : solver%single_precision = .true.
46 91 : solver%transform = .true.
47 91 : solver%GPU = .true.
48 91 : end function
49 :
50 : subroutine init()
51 : #ifdef CPP_MAGMA
52 : if (.not. initialized) then
53 : initialized = .true.
54 : call magmaf_init()
55 : call magmaf_setdevice(acc_get_device_num(acc_device_nvidia))
56 : print *, acc_get_device_num(acc_device_nvidia)
57 : end if
58 : #endif
59 : end subroutine
60 :
61 0 : subroutine magma_gev(self, hmat, smat, ne, eig, zmat, ikpt)
62 :
63 : use m_types_mat
64 : implicit none
65 :
66 : ! ... Arguments ...
67 : class(t_solver_magma) :: self
68 : class(t_mat), intent(INOUT) :: hmat, smat
69 : integer, intent(INOUT) :: ne
70 : class(t_mat), allocatable, intent(OUT) :: zmat
71 : real, intent(OUT) :: eig(:)
72 : integer, intent(IN) :: ikpt
73 :
74 : #ifdef CPP_MAGMA
75 :
76 : ! ... Local Variables ..
77 : integer :: lwork, liwork, lrwork, error, mout(1), i
78 : real :: eigTemp(hmat%matsize1)
79 : logical :: initialized = .false.
80 :
81 : real, allocatable :: rwork(:)
82 : integer, allocatable :: iwork(:)
83 : complex, allocatable :: work(:)
84 :
85 : call init()
86 :
87 : if (hmat%l_real) then
88 : allocate (rwork(1), iwork(1))
89 : !CALL magmaf_dsygvdx_m(Magma_numGPU,1,'v','i','U',hmat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,&
90 : ! SIZE(smat%data_r,1),0.0,0.0,1,ne,mout,eigTemp,rwork,-1,iwork,-1,error)
91 : call magmaf_dsygvdx(1, 'v', 'i', 'U', hmat%matsize1, hmat%data_r, size(hmat%data_r, 1), smat%data_r, &
92 : size(smat%data_r, 1), 0.0, 0.0, 1, ne, mout, eigTemp, rwork, -1, iwork, -1, error)
93 : if (error /= 0) then
94 : write (*, *) 'magmaf_dsygvdx error code: ', error
95 : call juDFT_error("Failed to query workspaces (1)", calledby="magma.F90")
96 : end if
97 : print *, "Magma1"
98 : lrwork = rwork(1)
99 : liwork = iwork(1)
100 : deallocate (rwork, iwork)
101 : allocate (rwork(lrwork), iwork(liwork))
102 : ! CALL magmaf_dsygvdx_m(Magma_numGPU,1,'v','i','U',hmat%matsize1,hmat%data_r,SIZE(hmat%data_r,1),smat%data_r,&
103 : ! SIZE(smat%data_r,1),0.0,0.0,1,ne,mout,eigTemp,rwork,lrwork,iwork,liwork,error)
104 : call magmaf_dsygvdx(1, 'v', 'i', 'U', hmat%matsize1, hmat%data_r, size(hmat%data_r, 1), smat%data_r, &
105 : size(smat%data_r, 1), 0.0, 0.0, 1, ne, mout, eigTemp, rwork, lrwork, iwork, liwork, error)
106 : print *, "Magma2"
107 : if (error /= 0) then
108 : write (*, *) 'magmaf_dsygvdx error code: ', error
109 : call juDFT_error("Magma failed to diagonalize Hamiltonian (1)", calledby="magma.F90")
110 : end if
111 : else
112 : !Query the workspace size
113 : allocate (work(1), rwork(1), iwork(1))
114 : !CALL magmaf_zhegvdx_2stage_m(NGPU_CONST,&
115 : !CALL magmaf_zhegvdx_m(Magma_numGPU,1,'v','i','U',hmat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,&
116 : ! SIZE(smat%data_c,1),0.0,0.0,1,ne,mout,eigTemp,work,-1,rwork,-1,iwork,-1,error)
117 : call magmaf_zhegvdx(1, 'v', 'i', 'U', hmat%matsize1, hmat%data_c, size(hmat%data_c, 1), smat%data_c, &
118 : size(smat%data_c, 1), 0.0, 0.0, 1, ne, mout, eigTemp, work, -1, rwork, -1, iwork, -1, error)
119 : if (error /= 0) then
120 : write (*, *) 'magmaf_zhegvdx error code: ', error
121 : call juDFT_error("Failed to query workspaces (2)", calledby="magma.F90")
122 : end if
123 : print *, "Magma1"
124 :
125 : lwork = work(1)
126 : lrwork = rwork(1)
127 : liwork = iwork(1)
128 : deallocate (work, rwork, iwork)
129 : allocate (work(lwork), rwork(lrwork), iwork(liwork))
130 : !Now the diagonalization
131 : !CALL magmaf_zhegvdx_2stage_m(NGPU_CONST,&
132 : !CALL magmaf_zhegvdx_m(Magma_numGPU,1,'v','i','U',hmat%matsize1,hmat%data_c,SIZE(hmat%data_c,1),smat%data_c,&
133 : ! SIZE(smat%data_c,1),0.0,0.0,1,ne,mout,eigTemp,work,lwork,rwork,lrwork,iwork,liwork,error)
134 : call magmaf_zhegvdx(1, 'v', 'i', 'U', hmat%matsize1, hmat%data_c, size(hmat%data_c, 1), smat%data_c, &
135 : size(smat%data_c, 1), 0.0, 0.0, 1, ne, mout, eigTemp, work, lwork, rwork, lrwork, iwork, liwork, error)
136 :
137 : print *, "Magma2"
138 :
139 : if (error /= 0) then
140 : write (*, *) 'magmaf_zhegvdx error code: ', error
141 : call juDFT_error("Magma failed to diagonalize Hamiltonian (2)", calledby="magma.F90")
142 : end if
143 : end if
144 :
145 : allocate (t_mat::zmat)
146 : call zmat%alloc(hmat%l_real, hmat%matsize1, ne)
147 : do i = 1, ne
148 : eig(i) = eigTemp(i)
149 : if (hmat%l_real) then
150 : zmat%data_r(:, i) = hmat%data_r(:, i)
151 : else
152 : zmat%data_c(:, i) = hmat%data_c(:, i)
153 : end if
154 : end do
155 : #endif
156 0 : end subroutine magma_gev
157 :
158 :
159 0 : subroutine magma_diag(self, hmat, ne, eig, zmat)
160 : !Simple driver to solve Generalized Eigenvalue Problem using magma routine
161 : implicit none
162 : class(t_solver_magma) :: self
163 : class(t_mat), intent(INOUT) :: hmat
164 : integer, intent(INOUT) :: ne
165 : class(t_mat), allocatable, intent(OUT) :: zmat
166 : real, intent(OUT) :: eig(:)
167 :
168 : integer :: info, m, n
169 : real :: abstol
170 : real, external :: dlamch
171 : real :: eigTemp(hmat%matsize1)
172 : #ifdef CPP_MAGMA
173 : call init()
174 : n = hmat%matsize1
175 : if (n /= hmat%matsize2) call judft_error("Non-square matrix in magma_diag")
176 : allocate (t_mat::zmat)
177 : call zmat%alloc(hmat%l_real, n, ne)
178 : abstol = 2*dlamch('S')
179 : if (hmat%l_real) then
180 : block !workspace locally
181 : integer:: isuppz(2*n), lrwork, liwork(1)
182 : real :: rwork_dum(1)
183 : real, allocatable :: rwork(:)
184 : integer, allocatable :: iwork(:)
185 : ! Workspace query
186 : call magmaf_dsyevr('V', 'I', 'U', n, hmat%data_r, size(hmat%data_r,1),&
187 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_r, size(zmat%data_r,1), &
188 : isuppz, rwork_dum, -1, liwork, -1, info)
189 : lrwork = rwork_dum(1)
190 : allocate (rwork(lrwork), iwork(liwork(1)))
191 : call magmaf_dsyevr('V', 'I', 'U', n, hmat%data_r, size(hmat%data_r,1), &
192 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_r, size(zmat%data_r,1), &
193 : isuppz, rwork, lrwork, iwork, liwork(1), info)
194 : end block
195 : else
196 : block !workspace locally
197 : integer:: isuppz(2*n), lwork, lrwork, liwork(1)
198 : complex:: work_dum(1)
199 : real :: rwork_dum(1)
200 : complex, allocatable :: work(:)
201 : real, allocatable :: rwork(:)
202 : integer, allocatable :: iwork(:)
203 : ! Workspace query
204 : call magmaf_zheevr('V', 'I', 'U', n, hmat%data_c, size(hmat%data_c,1), &
205 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_c, size(zmat%data_c,1), isuppz, work_dum, &
206 : -1, rwork_dum, -1, liwork, -1, info)
207 : lwork = work_dum(1)
208 : lrwork = rwork_dum(1)
209 : allocate (work(lwork), rwork(lrwork), iwork(liwork(1)))
210 : call magmaf_zheevr('V', 'I', 'U', n, hmat%data_c, size(hmat%data_c,1), &
211 : 0.0, 0.0, 1, ne, abstol, m, eigTemp, zmat%data_c, size(zmat%data_c,1), isuppz, work, &
212 : lwork, rwork, lrwork, iwork, liwork(1), info)
213 : end block
214 : end if
215 : eig(:min(size(eig), size(eigTemp))) = eigTemp(:min(size(eig), size(eigTemp)))
216 : #endif
217 0 : end subroutine magma_diag
218 0 : subroutine magma_diag_sp(self, hmat, ne, eig, zmat)
219 : !Simple driver to solve Standard Eigenvalue Problem using magma routine
220 : implicit none
221 : class(t_solver_magma) :: self
222 : class(t_mat), intent(INOUT) :: hmat
223 : integer, intent(INOUT) :: ne
224 : class(t_mat), allocatable, intent(OUT) :: zmat
225 : real, intent(OUT) :: eig(:)
226 : #ifdef CPP_MAGMA
227 : integer, parameter:: sp = selected_real_kind(6)
228 : integer :: info, m, n ,lwork
229 : real(sp) :: eigval(hmat%matsize1)
230 : call init()
231 : n = hmat%matsize1
232 : if (n /= hmat%matsize2) call judft_error("Non-square matrix in magma_diag")
233 : allocate (t_mat::zmat)
234 : call zmat%alloc(hmat%l_real, n, ne)
235 :
236 : if (hmat%l_real) then
237 : BLOCK
238 : REAL(kind=sp),allocatable:: h(:,:),z(:,:),eigval(:),work(:)
239 : integer,allocatable :: iwork(:),ifail(:)
240 : Allocate(h(size(hmat%data_r,1),size(hmat%data_r,2)))
241 : Allocate(eigval(size(hmat%data_r,1)),ifail(size(hmat%data_r,1)))
242 : Allocate(z(size(hmat%data_r,1),ne))
243 : h=hmat%data_r
244 :
245 : allocate(work(1),iwork(5*size(h,1)))
246 : call magmaf_ssyevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,1.0E-8_sp,m,eigval,z,size(z,1),work,-1,iwork,ifail,info)
247 : lwork=work(1)
248 : deallocate(work)
249 : allocate(work(lwork))
250 :
251 : call magmaf_ssyevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,1.0E-8_sp,m,eigval,z,size(z,1),work,lwork,iwork,ifail,info)
252 :
253 : eig(:ne)=eigval(:ne)
254 : zmat%data_r=z(:,:ne)
255 : deallocate(h,z,eigval,work,iwork)
256 : END BLOCK
257 : else
258 : BLOCK
259 : COMPLEX(kind=sp),allocatable:: h(:,:),z(:,:),work(:)
260 : REAL(kind=sp),allocatable:: eigval(:),rwork(:)
261 : integer,allocatable :: iwork(:),ifail(:)
262 : Allocate(h(size(hmat%data_c,1),size(hmat%data_c,2)))
263 : Allocate(eigval(size(hmat%data_c,1)),ifail(size(hmat%data_c,1)))
264 : Allocate(z(size(hmat%data_c,1),ne),rwork(7*size(hmat%data_c,1)))
265 : h=hmat%data_c
266 :
267 : allocate(work(1),iwork(5*size(hmat%data_c,1)))
268 : call magmaf_cheevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,m,eigval,z,size(z,1),work,-1,rwork,iwork,ifail,info)
269 : lwork=work(1)
270 : deallocate(work)
271 : allocate(work(lwork))
272 :
273 : call magmaf_cheevx('V','I','U',size(h,1),h,size(h,1),0.0,0.0,1,ne,0.0,m,eigval,z,size(z,1),work,lwork,rwork,iwork,ifail,info)
274 : eig=eigval(:ne)
275 : zmat%data_c=z(:,:ne)
276 : deallocate(h,z,eigval,work,rwork,iwork)
277 : END BLOCK
278 : end if
279 : #endif
280 0 : end subroutine magma_diag_sp
281 :
282 0 : subroutine magma_reduction(self, hmat, smat)
283 : !Simple driver to solve Generalized Eigenvalue Problem using magma routine
284 : class(t_solver_magma) :: self
285 : class(t_mat), intent(INOUT) :: hmat, smat
286 :
287 : integer :: info, n
288 : #ifdef CPP_MAGMA
289 : call init()
290 : n = smat%matsize1 !Matrix size
291 : if (n /= smat%matsize2 .or. n /= hmat%matsize1 .or. n /= hmat%matsize2) &
292 : call judft_error("Matices not square in magma_reduction")
293 : if (smat%l_real) then
294 : ! Perform Cholesky decomposition of B to obtain L (B = L * L^T)
295 : call magmaf_dpotrf('U', n, smat%data_r, size(smat%data_r,1), info)
296 : if (info /= 0) call juDFT_error("Error in Cholesky decomposition of B")
297 :
298 : ! Transform A to A' = L^-1 * A * L^-T using chegst
299 : call magmaf_dsygst(1, "U", n, hmat%data_r, size(hmat%data_r,1), smat%data_r, size(smat%data_r,1), info)
300 : if (info /= 0) call juDFT_error("Error in dsygst")
301 :
302 : else
303 : ! Perform Cholesky decomposition of B to obtain L (B = L * L^T)
304 : call magmaf_zpotrf('U', n, smat%data_c, size(smat%data_c,1), info)
305 : if (info /= 0) call juDFT_error("Error in Cholesky decomposition of B")
306 :
307 : ! Transform A to A' = L^-1 * A * L^-T using chegst
308 : call magmaf_zhegst(1, "U", n, hmat%data_c, size(hmat%data_c,1), smat%data_c, size(smat%data_c,1), info)
309 : if (info /= 0) call juDFT_error("Error in zhegst")
310 : end if
311 : #endif
312 0 : end subroutine magma_reduction
313 :
314 0 : subroutine magma_recover(self, smat, zmat)
315 : class(t_solver_magma) :: self
316 : class(t_mat), intent(INOUT) :: zmat, smat
317 : integer :: m, n, info
318 : #ifdef CPP_MAGMA
319 : call init()
320 : n = smat%matsize1
321 : m = zmat%matsize2
322 : if (n /= smat%matsize2 .or. n /= zmat%matsize1) call judft_error("Invalid matix sizes in reduction_magma")
323 : if (smat%l_real) then
324 : ! recover the generalized eigenvectors z by solving z' = l^t * z
325 : call magmaf_dtrtrs('U', 'N', 'N', n, m, smat%data_r, n, zMat%data_r, n, info)
326 : if (info /= 0) call juDFT_error("Error in back transformation (dpotrs)")
327 : else
328 : ! --> recover the generalized eigenvectors z by solving z' = l^t * z
329 : call magmaf_ztrtrs('U', 'N', 'N', n, m, smat%data_c, n, zMat%data_c, n, info)
330 : if (info /= 0) call juDFT_error("Error in back transformation (zpotrs)")
331 : end if
332 : #endif
333 0 : end subroutine
334 91 : end module m_magma
|
