Line data Source code
1 : ! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
2 : ! This file is part of FLEUR and available as free software under the conditions
3 : ! of the MIT license as expressed in the LICENSE file in more detail.
4 : !
5 : ! @authors: Miriam Hinzen, Gregor Michalicek
6 : ! Added MPI implementation, DW 2018
7 : !--------------------------------------------------------------------------------
8 : module m_chase
9 :
10 : use m_judft
11 : use m_constants
12 : use m_types_mat
13 : use m_types_mpimat
14 : use m_types_solver
15 : use , INTRINSIC :: iso_c_binding,ONLY: c_char
16 : #ifdef CPP_MPI
17 : use mpi
18 : #endif
19 : #ifdef CPP_CHASE
20 : use chase_diag !this is the interface to chase provided in the chase library
21 : #endif
22 : implicit none
23 :
24 : private
25 : type, extends(t_solver)::t_solver_chase
26 : contains
27 : procedure :: solve_std_sp => chase_diag_sp
28 : procedure :: solve_std_dp => chase_diag_dp
29 : end type
30 : public :: get_solver_chase
31 :
32 : contains
33 :
34 91 : function get_solver_chase() result(solver)
35 : type(t_solver_chase), pointer::solver
36 91 : allocate (solver)
37 91 : solver%name = "chase"
38 : #ifdef CPP_CHASE
39 : solver%available = .true.
40 : #else
41 91 : solver%available = .false.
42 : #endif
43 91 : solver%parallel = .true.
44 91 : solver%serial = .true.
45 91 : solver%generalized = .false.
46 91 : solver%standard = .true.
47 91 : solver%single_precision = .true.
48 91 : solver%transform = .false.
49 91 : solver%GPU = .true.
50 91 : end function
51 :
52 0 : subroutine chase_diag_dp(self, hmat, ne, eig, zmat)
53 : implicit none
54 : class(t_solver_chase) :: self
55 : class(t_mat), intent(INOUT) :: hmat
56 : integer, intent(INOUT) :: ne
57 : class(t_mat), allocatable, intent(OUT) :: zmat
58 : real, intent(OUT) :: eig(:)
59 :
60 : select type (hmat)
61 : type is (t_mat)
62 0 : call chase_serial_dp(hmat, ne, eig, zmat)
63 : type is (t_mpimat)
64 0 : call chase_mpi_dp(hmat, ne, eig, zmat)
65 : end select
66 0 : end subroutine
67 :
68 0 : subroutine chase_diag_sp(self, hmat, ne, eig, zmat)
69 : implicit none
70 : class(t_solver_chase) :: self
71 : class(t_mat), intent(INOUT) :: hmat
72 : integer, intent(INOUT) :: ne
73 : class(t_mat), allocatable, intent(OUT) :: zmat
74 : real, intent(OUT) :: eig(:)
75 :
76 : select type (hmat)
77 : type is (t_mat)
78 0 : call chase_serial_sp(hmat, ne, eig, zmat)
79 : type is (t_mpimat)
80 0 : call judft_bug("Chase in SP for parallel matrices not implemented")
81 : !call chase_mpi_diag_sp(hmat,ne,eig,zmat)
82 : end select
83 0 : end subroutine
84 :
85 0 : subroutine chase_serial_dp(hmat, ne, eig, zmat)
86 : !Simple driver to solve Standard Eigenvalue Problem using ChASE routine
87 : implicit none
88 : type(t_mat), intent(INOUT),VOLATILE :: hmat
89 : integer, intent(INOUT) :: ne
90 : class(t_mat), allocatable, intent(OUT) :: zmat
91 : real, intent(OUT) :: eig(:)
92 :
93 : !These chase parameters might to be adjusted
94 : real, parameter :: tol = 1e-10
95 : character(kind=c_char), parameter :: mode = 'R', opt = 'S',qr='C'
96 : integer, parameter :: deg = 20
97 : #ifdef CPP_CHASE
98 : integer :: nex !extra search space
99 : integer :: init !status variable
100 : !chase will modify these variables in call to xchase even though these are not arguments!!
101 : real, allocatable, VOLATILE :: zr(:, :), eigval(:)
102 : complex, allocatable, VOLATILE :: zc(:, :)
103 : nex = 0.2*ne
104 : allocate (eigval(ne+nex))
105 : allocate (t_mat::zmat)
106 : call zmat%init(hmat%l_real, hmat%matsize1, ne)
107 : call hmat%u2l() !chase needs full matrix not only upper part!
108 : if (hmat%l_real) then
109 : allocate (zr(hmat%matsize1, ne + nex))
110 : ! Initialize of ChASE
111 :
112 : call dchase_init(size(hmat%data_r,1), ne, nex, hmat%data_r, size(hmat%data_r,1),zr, eigval, init)
113 : !Solve eigenvalue problem
114 : call dchase(deg, tol, mode, opt,qr)
115 : ! finalize and clean up
116 : call dchase_finalize(init)
117 :
118 : zmat%data_r = zr(:, :ne)
119 : else
120 : allocate (zc(hmat%matsize1, ne + nex))
121 : ! Initialize of ChASE
122 : call zchase_init(hmat%matsize1, ne, nex, hmat%data_c, size(hmat%data_c,1), zc, eigval, init)
123 : !Solve eigenvalue problem
124 : call zchase(deg, tol, mode, opt,qr)
125 : ! finalize and clean up
126 : call zchase_finalize(init)
127 : zmat%data_c = zc(:, :ne)
128 : end if
129 : eig(:ne) = eigval(:ne)
130 : #endif
131 : end subroutine chase_serial_dp
132 :
133 : subroutine chase_serial_sp(hmat, ne, eig, zmat)
134 : !Simple driver to solve Standard Eigenvalue Problem using ChASE routine in single precision
135 : implicit none
136 : type(t_mat), intent(INOUT) :: hmat
137 : integer, intent(INOUT) :: ne
138 : class(t_mat), allocatable, intent(OUT) :: zmat
139 : real, intent(OUT) :: eig(:)
140 :
141 : integer, parameter:: sp = selected_real_kind(6)
142 : !These chase parameters might to be adjusted
143 : real(sp), parameter :: tol = 1e-6
144 : character, parameter :: mode = 'R', opt = 'S', qr='N'
145 : integer, parameter :: deg = 20
146 :
147 : #ifdef CPP_CHASE
148 : integer :: nex !extra search space
149 : integer :: init !status variable
150 : !chase will modify these variables in call to xchase even though these are not arguments!!
151 : real(sp), allocatable, volatile :: zr(:, :), eigval(:)
152 : complex(sp), allocatable, volatile :: zc(:, :)
153 : real(sp), allocatable :: hr(:, :)
154 : complex(sp), allocatable :: hc(:, :)
155 : nex = 0.2*ne
156 : allocate (eigval(ne))
157 : allocate (t_mat::zmat)
158 : call zmat%init(hmat%l_real, hmat%matsize1, ne)
159 :
160 : call hmat%u2l() !chase needs full matrix not only upper part!
161 : if (hmat%l_real) then
162 : allocate (zr(hmat%matsize1, ne + nex))
163 : allocate (hr(hmat%matsize1, hmat%matsize2))
164 : hr = hmat%data_r !cast to sp
165 : ! Initialize of ChASE
166 : call schase_init(hmat%matsize1, ne, nex, hr,size(hr,1), zr, eigval, init)
167 : !Solve eigenvalue problem
168 : call schase(deg, tol, mode, opt,qr)
169 : ! finalize and clean up
170 : call schase_finalize(init)
171 : zmat%data_r = zr(:, :ne)
172 : else
173 : allocate (zc(hmat%matsize1, ne + nex))
174 : allocate (hc(hmat%matsize1, hmat%matsize2))
175 : hc = hmat%data_c
176 : ! Initialize of ChASE
177 : call cchase_init(hmat%matsize1, ne, nex, hc, size(hc,1),zc, eigval, init)
178 : !Solve eigenvalue problem
179 : call cchase(deg, tol, mode, opt,qr)
180 : ! finalize and clean up
181 : call cchase_finalize(init)
182 : zmat%data_c = zc(:, :ne)
183 : end if
184 : eig(:ne) = eigval(:ne)
185 : #endif
186 : end subroutine chase_serial_sp
187 :
188 : subroutine chase_mpi_dp(hmat, ne, eig, zmat)
189 : !Simple driver to solve Standard Eigenvalue Problem using ChASE routine
190 : implicit none
191 : type(t_mpimat), intent(INOUT) :: hmat
192 : integer, intent(INOUT) :: ne
193 : class(t_mat), allocatable, intent(OUT) :: zmat
194 : real, intent(OUT) :: eig(:)
195 :
196 : !These chase parameters might to be adjusted
197 : real, parameter :: tol = 1e-10
198 : character, parameter :: mode = 'R', opt = 'S', qr='N'
199 : character, parameter :: grid_major = "C" !major of 2D MPI grid. Row major: grid_major=’R’, column major: grid_major=’C’
200 : integer, parameter :: deg = 20
201 : #ifdef CPP_CHASE
202 : integer:: mbsize, nbsize, irsrc, icsrc, dim0, dim1, myprow, mypcol
203 : integer :: comm_1d, comm_2d, ierr
204 : integer :: nex !extra search space
205 : integer :: init !status variable
206 : !chase will modify these variables in call to xchase even though these are not arguments!!
207 : real, allocatable, volatile :: eigval(:)
208 : type(t_mpimat), volatile :: ztemp
209 : nex = 0.2*ne
210 : allocate (eigval(ne))
211 :
212 : !setup ChASE
213 : mbsize = hmat%blacsdata%blacs_desc(5) !block size for the block-cyclic distribution for the rows of global matrix
214 : nbsize = hmat%blacsdata%blacs_desc(6) !block size for the block-cyclic distribution for the cloumns of global matrix
215 : irsrc = hmat%blacsdata%blacs_desc(7) !process row over which the first row of the global matrix h is distributed
216 : icsrc = hmat%blacsdata%blacs_desc(8) !process column over which the first column of the global matrix h is distributed.
217 :
218 : !determine the processor grid
219 : !dim0 :row number of 2D MPI grid
220 : !dim1 :column number of 2D MPI grid
221 : call BLACS_GRIDINFO(hmat%blacsdata%blacs_desc(2), dim0, dim1, MYPROW, MYPCOL)
222 : call create_mpi_comms(hmat%blacsdata%mpi_com, hmat%blacsdata%blacs_desc(2), comm_2d, comm_1d)
223 :
224 : call ztemp%init(hmat%l_real, hmat%global_size1, ne + nex, comm_1d, MPIMAT_COLUMN_BLOCK_CYCLIC)
225 :
226 : call hmat%u2l() !chase needs full matrix not only upper part!
227 : if (hmat%l_real) then
228 : ! Initialize of ChASE
229 : call pdchase_init_blockcyclic(hmat%global_size1, ne, nex, mbsize, nbsize, hmat%data_r, hmat%matsize1, &
230 : ztemp%data_r, eigval, dim0, dim1, grid_major, irsrc, icsrc, comm_2d, init)
231 : !Solve eigenvalue problem
232 : call dchase(deg, tol, mode, opt,qr)
233 : ! finalize and clean up
234 : call dchase_finalize(init)
235 : else
236 : ! Initialize of ChASE
237 : call pzchase_init_blockcyclic(hmat%global_size1, ne, nex, mbsize, nbsize, hmat%data_c, hmat%matsize1, &
238 : ztemp%data_c, eigval, dim0, dim1, grid_major, irsrc, icsrc, comm_2d, init)
239 : !Solve eigenvalue problem
240 : call zchase(deg, tol, mode, opt,qr)
241 : ! finalize and clean up
242 : call zchase_finalize(init)
243 : end if
244 : !create zmat in correct distribution
245 : allocate (t_mpimat::zmat)
246 : call zmat%init(hmat%l_real, hmat%matsize1, ne + NEX, hmat%blacsdata%mpi_com, MPIMAT_ROWCYCLIC)
247 : call zmat%copy(ztemp, 1, 1)
248 : eig(:ne) = eigval(:ne)
249 :
250 : call MPI_COMM_FREE(comm_1d, ierr)
251 : call MPI_COMM_FREE(comm_2d, ierr)
252 : #endif
253 : end subroutine chase_mpi_dp
254 :
255 : subroutine create_mpi_comms(parent_comm, icontext, comm_2d, comm_1d)
256 : integer, intent(in) :: parent_comm, icontext
257 : integer, intent(out):: comm_2d, comm_1d
258 : #ifdef CPP_CHASE
259 :
260 : integer:: isize, ierr, i, col, row, no_col, no_row, myprow, mypcol
261 : integer, allocatable:: map1d(:), map2d(:)
262 : integer :: group, group_2d, group_1d
263 :
264 : call BLACS_GRIDINFO(icontext, no_row, no_col, MYPROW, MYPCOL)
265 : call MPI_COMM_SIZE(parent_comm, isize, ierr)
266 : allocate (map2d(isize))
267 : allocate (map1d(no_row))
268 :
269 : map1d = isize !largest value
270 : do i = 0, isize - 1
271 : call blacs_pcoord(icontext, i, ROW, COL)
272 : map2d(row + no_row*(col - 1)) = i
273 : map1d(row) = min(map1d(row), i)
274 : end do
275 :
276 : !create 2d communicator
277 : call MPI_COMM_group(parent_comm, group, ierr)
278 : call MPI_Group_incl(group, isize, map2d, group_2d, ierr)
279 : call MPI_COMM_create_group(parent_comm, 1,group_2d, comm_2d, ierr)
280 :
281 : !create 1d communicator
282 : call MPI_COMM_group(parent_comm, group, ierr)
283 : call MPI_Group_incl(group, size(map1d), map1d, group_1d, ierr)
284 : call MPI_COMM_create_group(parent_comm, 1,group_1d, comm_1d, ierr)
285 :
286 : call MPI_group_free(group_2d, ierr)
287 : call MPI_group_free(group_1d, ierr)
288 : call MPI_group_free(group, ierr)
289 : #endif
290 :
291 : end subroutine
292 :
293 0 : end module m_chase
|
