Line data Source code
1 : !!****m* ABINIT/m_xgScalapack
2 : !! NAME
3 : !! m_xgScalapack
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2017-2026 ABINIT group (J. Bieder)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_xgScalapack
22 :
23 : use defs_basis, only : std_err, std_out, dp, ABI_GPU_DISABLED, ABI_GPU_OPENMP
24 : use m_abicore
25 : USE_MPI
26 : use m_xmpi
27 : use m_errors
28 : use m_slk
29 : use m_xg
30 : use m_xomp
31 : use m_time, only: timab
32 :
33 : implicit none
34 :
35 : #ifdef HAVE_MPI1
36 : include 'mpif.h'
37 : #endif
38 :
39 :
40 : private
41 :
42 : integer, parameter :: M__SLK = 1
43 : integer, parameter :: M__ROW = 1
44 : integer, parameter :: M__COL = 2
45 : integer, parameter :: M__UNUSED = 4
46 : integer, parameter :: M__WORLD = 5
47 : integer, parameter :: M__NDATA = 5
48 : integer, parameter :: M__tim_init = 1690
49 : integer, parameter :: M__tim_free = 1691
50 : integer, parameter :: M__tim_heev = 1692
51 : integer, parameter :: M__tim_hegv = 1693
52 : integer, parameter :: M__tim_scatter = 1694
53 :
54 : integer, parameter, public :: SLK_AUTO = -1
55 : integer, parameter, public :: SLK_FORCED = 1
56 : integer, parameter, public :: SLK_DISABLED = 0
57 : integer, save :: M__CONFIG = SLK_AUTO
58 : integer, save :: M__MAXDIM = 1000
59 :
60 : type, public :: xgScalapack_t
61 : integer :: comms(M__NDATA)
62 : integer :: rank(M__NDATA)
63 : integer :: size(M__NDATA)
64 : integer :: coords(2)
65 : integer :: ngroup
66 : integer :: verbosity
67 : type(slk_grid_t) :: grid
68 : end type xgScalapack_t
69 :
70 : public :: xgScalapack_init
71 : public :: xgScalapack_free
72 : public :: xgScalapack_heev
73 : public :: xgScalapack_hegv
74 : public :: xgScalapack_config
75 : contains
76 : !!***
77 :
78 : !!****f* m_xgScalapack/xgScalapack_init
79 : !! NAME
80 : !! xgScalapack_init
81 : !!
82 : !! FUNCTION
83 : !! Init the scalapack communicator for next operations.
84 : !! If the comm has too many cpus, then take only a subgroup of this comm
85 : !!
86 : !! INPUTS
87 : !!
88 : !! OUTPUT
89 : !!
90 : !! SOURCE
91 0 : subroutine xgScalapack_init(xgScalapack,comm,maxDim,verbosity,gpu_option,usable)
92 :
93 : type(xgScalapack_t), intent(inout) :: xgScalapack
94 : integer , intent(in ) :: comm
95 : integer , intent(in ) :: maxDim
96 : integer , intent(in ) :: verbosity
97 : logical , intent(in ) :: gpu_option
98 : logical , intent( out) :: usable
99 : double precision :: tsec(2)
100 : #ifdef HAVE_LINALG_MKL_THREADS
101 : integer :: mkl_get_max_threads
102 : #endif
103 : #ifdef HAVE_LINALG_OPENBLAS_THREADS
104 : integer :: openblas_get_num_threads
105 : #endif
106 : #ifdef HAVE_LINALG_NVPL_THREADS
107 : integer :: nvpl_blas_get_max_threads
108 : #endif
109 : integer :: nthread
110 : #ifdef HAVE_LINALG_SCALAPACK
111 : integer :: maxProc
112 : integer :: nproc
113 : integer :: ngroup
114 : integer :: subgroup
115 : integer :: mycomm(2)
116 : integer :: ierr
117 : integer :: test_row
118 : integer :: test_col
119 : #else
120 : ABI_UNUSED(comm)
121 : ABI_UNUSED(maxDim)
122 : #endif
123 :
124 0 : call timab(M__tim_init,1,tsec)
125 :
126 0 : xgScalapack%comms = xmpi_comm_null
127 0 : xgScalapack%rank = xmpi_undefined_rank
128 0 : xgScalapack%verbosity = verbosity
129 :
130 0 : nthread = 1
131 : #ifdef HAVE_LINALG_MKL_THREADS
132 0 : nthread = mkl_get_max_threads()
133 : #elif HAVE_LINALG_OPENBLAS_THREADS
134 : nthread = openblas_get_num_threads()
135 : #elif HAVE_LINALG_NVPL_THREADS
136 : nthread = nvpl_blas_get_max_threads()
137 : #else
138 : nthread = xomp_get_num_threads(open_parallel=.true.)
139 : if ( nthread == 0 ) nthread = 1
140 : #endif
141 :
142 : #ifdef HAVE_LINALG_SCALAPACK
143 :
144 : nproc = xmpi_comm_size(comm)
145 : xgScalapack%comms(M__WORLD) = comm
146 : xgScalapack%rank(M__WORLD) = xmpi_comm_rank(comm)
147 : xgScalapack%size(M__WORLD) = nproc
148 :
149 : maxProc = (maxDim / (M__MAXDIM*nthread))+1 ! ( M__MAXDIM x M__MAXDIM matrice per MPI )
150 : if ( M__CONFIG > 0 .and. M__CONFIG <= nproc ) then
151 : maxProc = M__CONFIG
152 : else if ( maxProc > nproc ) then
153 : maxProc = nproc
154 : end if
155 :
156 : if ( maxProc == 1 .or. M__CONFIG == SLK_DISABLED) then
157 : usable = .false.
158 : return
159 : else if ( nthread > 1 ) then ! disable scalapack with threads since it is not threadsafe
160 : ! This should be check with new elpa version en MPI+OpenMP
161 : if ( M__CONFIG > 0 ) then
162 : ABI_WARNING("xgScalapack turned off because you have threads")
163 : end if
164 : usable = .false.
165 : return
166 : else
167 : usable = .true.
168 : maxProc = 2*((maxProc+1)/2) ! Round to next even number
169 : end if
170 :
171 : if ( xgScalapack%verbosity > 0 ) then
172 : write(std_out,*) " xgScalapack will use", maxProc, "/", nproc, "MPIs"
173 : end if
174 :
175 : ngroup = nproc/maxProc
176 : xgScalapack%ngroup = ngroup
177 :
178 : if ( maxProc < nproc ) then
179 : if ( xgScalapack%rank(M__WORLD) < maxProc*ngroup ) then
180 : subgroup = xgScalapack%rank(M__WORLD)/maxProc
181 : mycomm(1) = M__SLK
182 : mycomm(2) = M__UNUSED
183 : else
184 : subgroup = ngroup+1
185 : mycomm(1) = M__UNUSED
186 : mycomm(2) = M__SLK
187 : end if
188 : call MPI_Comm_split(comm, subgroup, xgScalapack%rank(M__WORLD), xgScalapack%comms(mycomm(1)),ierr)
189 : if ( ierr /= 0 ) then
190 : ABI_ERROR("Error splitting communicator")
191 : end if
192 : xgScalapack%comms(mycomm(2)) = xmpi_comm_null
193 : xgScalapack%rank(mycomm(1)) = xmpi_comm_rank(xgScalapack%comms(mycomm(1)))
194 : xgScalapack%rank(mycomm(2)) = xmpi_undefined_rank
195 : xgScalapack%size(mycomm(1)) = xmpi_comm_size(xgScalapack%comms(mycomm(1)))
196 : xgScalapack%size(mycomm(2)) = nproc - xgScalapack%size(mycomm(1))
197 : else
198 : call MPI_Comm_dup(comm,xgScalapack%comms(M__SLK),ierr)
199 : if ( ierr /= 0 ) then
200 : ABI_ERROR("Error duplicating communicator")
201 : end if
202 : xgScalapack%rank(M__SLK) = xmpi_comm_rank(xgScalapack%comms(M__SLK))
203 : xgScalapack%size(M__SLK) = nproc
204 : end if
205 :
206 : if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
207 : call xgScalapack%grid%init(xgScalapack%size(M__SLK), xgScalapack%comms(M__SLK), gpu_option)
208 : call BLACS_GridInfo(xgScalapack%grid%comm, &
209 : xgScalapack%grid%dims(M__ROW), xgScalapack%grid%dims(M__COL),&
210 : xgScalapack%coords(M__ROW), xgScalapack%coords(M__COL))
211 :
212 : !These values are the same as those computed by BLACS_GRIDINFO
213 : !except in the case where the my_rank argument is not the local proc
214 : test_row = INT((xgScalapack%rank(M__SLK)) / xgScalapack%grid%dims(2))
215 : test_col = MOD((xgScalapack%rank(M__SLK)), xgScalapack%grid%dims(2))
216 : if ( test_row /= xgScalapack%coords(M__ROW) ) then
217 : ABI_WARNING("Row id mismatch")
218 : end if
219 : if ( test_col /= xgScalapack%coords(M__COL) ) then
220 : ABI_WARNING("Col id mismatch")
221 : end if
222 : end if
223 :
224 : #else
225 0 : usable = .false.
226 : #endif
227 :
228 0 : call timab(M__tim_init,2,tsec)
229 :
230 0 : end subroutine xgScalapack_init
231 :
232 5281 : subroutine xgScalapack_config(myconfig,maxDim)
233 :
234 : integer, intent(in) :: myconfig
235 : integer, intent(in) :: maxDim
236 5281 : if ( myconfig == SLK_AUTO) then
237 5237 : M__CONFIG = myconfig
238 5237 : ABI_COMMENT("xgScalapack in auto mode")
239 44 : else if ( myconfig == SLK_DISABLED) then
240 0 : M__CONFIG = myconfig
241 0 : ABI_COMMENT("xgScalapack disabled")
242 44 : else if ( myconfig > 0) then
243 44 : M__CONFIG = myconfig
244 44 : ABI_COMMENT("xgScalapack enabled")
245 : else
246 0 : ABI_WARNING("Bad value for xgScalapack config -> autodetection")
247 0 : M__CONFIG = SLK_AUTO
248 : end if
249 5281 : if ( maxDim > 0 ) then
250 5281 : M__MAXDIM = maxDim
251 : end if
252 :
253 5281 : end subroutine xgScalapack_config
254 :
255 : function toProcessorScalapack(xgScalapack) result(processor)
256 :
257 : type(xgScalapack_t), intent(in) :: xgScalapack
258 : type(slk_processor_t) :: processor
259 :
260 : processor%my_rank = xgScalapack%rank(M__SLK)
261 : processor%comm = xgScalapack%comms(M__SLK)
262 : processor%coords = xgScalapack%coords
263 : processor%grid = xgScalapack%grid
264 : end function toProcessorScalapack
265 :
266 : !This is for testing purpose.
267 : !May not be optimal since I do not control old implementation but at least gives a reference.
268 0 : subroutine xgScalapack_heev(xgScalapack,matrixA,eigenvalues,gpu_option)
269 : use, intrinsic :: iso_c_binding
270 : type(xgScalapack_t), intent(inout) :: xgScalapack
271 : type(xgBlock_t) , intent(inout) :: matrixA
272 : type(xgBlock_t) , intent(inout) :: eigenvalues
273 : integer, optional , intent(in) :: gpu_option
274 : #ifdef HAVE_LINALG_SCALAPACK
275 : double precision, pointer :: matrix(:,:) !(cplex*nbli_global,nbco_global)
276 : double precision, pointer :: eigenvalues_tmp(:,:)
277 : double precision, pointer :: vector(:)
278 : double precision :: tsec(2)
279 : integer :: cplex
280 : integer :: istwf_k
281 : integer :: nbli_global, nbco_global
282 : type(c_ptr) :: cptr
283 : integer :: req(2), status(MPI_STATUS_SIZE,2), ierr
284 : integer :: l_gpu_option,l_use_gpu_elpa
285 : #endif
286 :
287 : #ifdef HAVE_LINALG_SCALAPACK
288 : call timab(M__tim_heev,1,tsec)
289 :
290 : l_gpu_option=ABI_GPU_DISABLED
291 : if (present(gpu_option)) then
292 : l_gpu_option = gpu_option
293 : end if
294 : l_use_gpu_elpa=0
295 : #ifdef HAVE_LINALG_ELPA
296 : if (l_gpu_option/=ABI_GPU_DISABLED) l_use_gpu_elpa=1
297 : #endif
298 :
299 : ! Keep only working processors
300 : if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
301 :
302 : call xgBlock_getSize(eigenvalues,nbli_global,nbco_global)
303 : if ( cols(matrixA) /= nbli_global ) then
304 : ABI_ERROR("Number of eigen values differ from number of vectors")
305 : end if
306 :
307 : if ( space(matrixA) == SPACE_C ) then
308 : cplex = 2
309 : istwf_k = 1
310 : else
311 : cplex = 1
312 : istwf_k = 2
313 : endif
314 :
315 : call xgBlock_getSize(matrixA,nbli_global,nbco_global)
316 :
317 : if(l_gpu_option==ABI_GPU_OPENMP) then
318 : call xgBlock_copy_from_gpu(matrixA)
319 : call xgBlock_copy_from_gpu(eigenvalues)
320 : end if
321 :
322 : call xgBlock_reverseMap(matrixA,matrix,rows=nbli_global,cols=nbco_global)
323 : call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,rows=nbco_global,cols=1)
324 : cptr = c_loc(eigenvalues_tmp)
325 : call c_f_pointer(cptr,vector,(/ nbco_global /))
326 :
327 : call compute_eigen1(xgScalapack%comms(M__SLK), &
328 : toProcessorScalapack(xgScalapack), &
329 : cplex,nbli_global,nbco_global,matrix,vector,istwf_k,&
330 : use_gpu_elpa=l_use_gpu_elpa)
331 :
332 : end if
333 :
334 : call timab(M__tim_heev,2,tsec)
335 :
336 : req(1:2)=-1
337 : call xgScalapack_scatter(xgScalapack,matrixA,req(1))
338 : call xgScalapack_scatter(xgScalapack,eigenvalues,req(2))
339 : #ifdef HAVE_MPI
340 : if ( any(req/=-1) ) then
341 : call MPI_WaitAll(2,req,status,ierr)
342 : if ( ierr /= 0 ) then
343 : ABI_ERROR("Error waiting data")
344 : endif
345 : end if
346 : #endif
347 :
348 : if(l_gpu_option==ABI_GPU_OPENMP) then
349 : call xgBlock_copy_to_gpu(matrixA)
350 : call xgBlock_copy_to_gpu(eigenvalues)
351 : end if
352 :
353 : #else
354 0 : ABI_ERROR("ScaLAPACK support not available")
355 : ABI_UNUSED(xgScalapack%verbosity)
356 : ABI_UNUSED(matrixA%normal)
357 : ABI_UNUSED(eigenvalues%normal)
358 : #endif
359 :
360 0 : end subroutine xgScalapack_heev
361 :
362 : !This is for testing purpose.
363 : !May not be optimal since I do not control old implementation but at least gives a reference.
364 0 : subroutine xgScalapack_hegv(xgScalapack,matrixA,matrixB,eigenvalues,gpu_option)
365 : use, intrinsic :: iso_c_binding
366 : type(xgScalapack_t), intent(inout) :: xgScalapack
367 : type(xgBlock_t) , intent(inout) :: matrixA
368 : type(xgBlock_t) , intent(inout) :: matrixB
369 : type(xgBlock_t) , intent(inout) :: eigenvalues
370 : integer, optional , intent(in) :: gpu_option
371 : #ifdef HAVE_LINALG_SCALAPACK
372 : double precision, pointer :: matrix1(:,:) !(cplex*nbli_global,nbco_global)
373 : double precision, pointer :: matrix2(:,:) !(cplex*nbli_global,nbco_global)
374 : double precision, pointer :: eigenvalues_tmp(:,:)
375 : double precision, pointer :: vector(:)
376 : double precision :: tsec(2)
377 : integer :: cplex
378 : integer :: istwf_k
379 : integer :: nbli_global, nbco_global
380 : type(c_ptr) :: cptr
381 : integer :: req(2), status(MPI_STATUS_SIZE,2),ierr
382 : integer :: l_gpu_option,l_use_gpu_elpa
383 : #endif
384 :
385 : #ifdef HAVE_LINALG_SCALAPACK
386 : call timab(M__tim_hegv,1,tsec)
387 :
388 : l_gpu_option=ABI_GPU_DISABLED
389 : if (present(gpu_option)) then
390 : l_gpu_option = gpu_option
391 : end if
392 : l_use_gpu_elpa=0
393 : #ifdef HAVE_LINALG_ELPA
394 : if (l_gpu_option/=ABI_GPU_DISABLED) l_use_gpu_elpa=1
395 : #endif
396 :
397 : ! Keep only working processors
398 : if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
399 :
400 : call xgBlock_getSize(eigenvalues,nbli_global,nbco_global)
401 : if ( cols(matrixA) /= cols(matrixB) ) then
402 : ABI_ERROR("Matrix A and B don't have the same number of vectors")
403 : end if
404 :
405 : if ( cols(matrixA) /= nbli_global ) then
406 : ABI_ERROR("Number of eigen values differ from number of vectors")
407 : end if
408 :
409 : if ( space(matrixA) == SPACE_C ) then
410 : cplex = 2
411 : istwf_k = 1
412 : else
413 : cplex = 1
414 : istwf_k = 2
415 : endif
416 :
417 : call xgBlock_getSize(matrixA,nbli_global,nbco_global)
418 :
419 : if(l_gpu_option==ABI_GPU_OPENMP) then
420 : call xgBlock_copy_from_gpu(matrixA)
421 : call xgBlock_copy_from_gpu(matrixB)
422 : call xgBlock_copy_from_gpu(eigenvalues)
423 : end if
424 :
425 : call xgBlock_reverseMap(matrixA,matrix1,rows=nbli_global,cols=nbco_global)
426 : call xgBlock_reverseMap(matrixB,matrix2,rows=nbli_global,cols=nbco_global)
427 : call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,rows=nbco_global,cols=1)
428 : cptr = c_loc(eigenvalues_tmp)
429 : call c_f_pointer(cptr,vector,(/ nbco_global /))
430 :
431 : call compute_eigen2(xgScalapack%comms(M__SLK), &
432 : toProcessorScalapack(xgScalapack), &
433 : cplex,nbli_global,nbco_global,matrix1,matrix2,vector,istwf_k,&
434 : use_gpu_elpa=l_use_gpu_elpa)
435 : end if
436 :
437 : call timab(M__tim_hegv,2,tsec)
438 :
439 : req(1:2)=-1
440 : call xgScalapack_scatter(xgScalapack,matrixA,req(1))
441 : call xgScalapack_scatter(xgScalapack,eigenvalues,req(2))
442 : #ifdef HAVE_MPI
443 : if ( any(req/=-1) ) then
444 : call MPI_WaitAll(2,req,status,ierr)
445 : if ( ierr /= 0 ) then
446 : ABI_ERROR("Error waiting data")
447 : endif
448 : end if
449 : #endif
450 :
451 : if(l_gpu_option==ABI_GPU_OPENMP) then
452 : call xgBlock_copy_to_gpu(matrixA)
453 : call xgBlock_copy_to_gpu(eigenvalues)
454 : end if
455 :
456 : #else
457 0 : ABI_ERROR("ScaLAPACK support not available")
458 : ABI_UNUSED(xgScalapack%verbosity)
459 : ABI_UNUSED(matrixA%normal)
460 : ABI_UNUSED(matrixB%normal)
461 : ABI_UNUSED(eigenvalues%normal)
462 : #endif
463 :
464 0 : end subroutine xgScalapack_hegv
465 :
466 :
467 : subroutine xgScalapack_scatter(xgScalapack,matrix,req)
468 :
469 : type(xgScalapack_t), intent(in ) :: xgScalapack
470 : type(xgBlock_t) , intent(inout) :: matrix
471 : integer , intent( out) :: req
472 : double precision, pointer :: tab(:,:)
473 : double precision :: tsec(2)
474 : !integer :: cols, rows
475 : integer :: ierr
476 : integer :: sendto, receivefrom
477 : integer :: lap
478 :
479 : call timab(M__tim_scatter,1,tsec)
480 :
481 : !call xgBlock_getSize(matrix,rows,cols)
482 : !call xgBlock_reverseMap(matrix,tab,rows,cols)
483 : call xgBlock_reverseMap(matrix,tab)
484 :
485 : ! If we did the he(e|g)v and we are the first group
486 : if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null .and. xgScalapack%rank(M__WORLD)<xgScalapack%size(M__SLK) ) then
487 : lap = xgScalapack%ngroup
488 : sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK)
489 : if ( sendto < xgScalapack%size(M__WORLD) ) then
490 : !do while ( sendto < xgScalapack%size(M__WORLD) )
491 : !call xmpi_send(tab,sendto,sendto,xgScalapack%comms(M__WORLD),ierr)
492 : call xmpi_isend(tab,sendto,sendto,xgScalapack%comms(M__WORLD),req,ierr)
493 : !write(*,*) xgScalapack%rank(M__WORLD), "sends to", sendto
494 : if ( ierr /= 0 ) then
495 : ABI_ERROR("Error sending data")
496 : end if
497 : !lap = lap+1
498 : !sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK)
499 : !end do
500 : end if
501 : else if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then
502 : receivefrom = MODULO(xgScalapack%rank(M__WORLD), xgScalapack%size(M__SLK))
503 : if ( receivefrom >= 0 ) then
504 : !call xmpi_recv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),ierr)
505 : call xmpi_irecv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),req,ierr)
506 : !write(*,*) xgScalapack%rank(M__WORLD), "receive from", receivefrom
507 : if ( ierr /= 0 ) then
508 : ABI_ERROR("Error receiving data")
509 : end if
510 : end if
511 : !else
512 : !ABI_BUG("Error scattering data")
513 : end if
514 :
515 : call timab(M__tim_scatter,2,tsec)
516 :
517 : end subroutine xgScalapack_scatter
518 :
519 :
520 0 : subroutine xgScalapack_free(xgScalapack)
521 :
522 : type(xgScalapack_t), intent(inout) :: xgScalapack
523 : double precision :: tsec(2)
524 : #ifdef HAVE_LINALG_SCALAPACK
525 : integer :: ierr
526 : #endif
527 :
528 0 : call timab(M__tim_free,1,tsec)
529 : #ifdef HAVE_LINALG_SCALAPACK
530 : if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
531 : call BLACS_GridExit(xgScalapack%grid%comm)
532 : call MPI_Comm_free(xgScalapack%comms(M__SLK),ierr)
533 : end if
534 : if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then
535 : call MPI_Comm_free(xgScalapack%comms(M__UNUSED),ierr)
536 : end if
537 : #else
538 : ABI_UNUSED(xgScalapack%verbosity)
539 : #endif
540 0 : call timab(M__tim_free,2,tsec)
541 :
542 0 : end subroutine xgScalapack_free
543 :
544 0 : end module m_xgScalapack
545 : !!***
|
