TABLE OF CONTENTS


ABINIT/m_xgScalapack [ Modules ]

[ Top ] [ Modules ]

NAME

  m_xgScalapack

FUNCTION

COPYRIGHT

  Copyright (C) 2017 ABINIT group (J. Bieder)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_xgScalapack
23 
24   use defs_basis, only : std_err, std_out, dp
25   use m_abicore
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 #ifdef HAVE_MPI2
34  use mpi
35 #endif
36 
37   implicit none
38 
39 #ifdef HAVE_MPI1
40  include 'mpif.h'
41 #endif
42 
43 
44   private
45 
46   integer, parameter :: M__SLK = 1
47   integer, parameter :: M__ROW = 1
48   integer, parameter :: M__COL = 2
49   integer, parameter :: M__UNUSED = 4
50   integer, parameter :: M__WORLD = 5
51   integer, parameter :: M__NDATA = 5
52   integer, parameter :: M__tim_init    = 1690
53   integer, parameter :: M__tim_free    = 1691
54   integer, parameter :: M__tim_heev    = 1692
55   integer, parameter :: M__tim_hegv    = 1693
56   integer, parameter :: M__tim_scatter = 1694
57 
58   integer, parameter, public :: SLK_AUTO = -1
59   integer, parameter, public :: SLK_FORCED = 1
60   integer, parameter, public :: SLK_DISABLED = 0
61   integer, save :: M__CONFIG = SLK_AUTO
62   integer, save :: M__MAXDIM = 1000
63 
64   type, public :: xgScalapack_t
65     integer :: comms(M__NDATA)
66     integer :: rank(M__NDATA)
67     integer :: size(M__NDATA)
68     integer :: coords(2)
69     integer :: ngroup
70     integer :: verbosity
71     type(grid_scalapack) :: grid
72   end type xgScalapack_t
73 
74   public :: xgScalapack_init
75   public :: xgScalapack_free
76   public :: xgScalapack_heev
77   public :: xgScalapack_hegv
78   public :: xgScalapack_config
79   contains

m_xgScalapack/xgScalapack_init [ Functions ]

[ Top ] [ m_xgScalapack ] [ Functions ]

NAME

  xgScalapack_init

FUNCTION

  Init the scalapack communicator for next operations.
  If the comm has to many cpus, then take only a subgroup of this comm

INPUTS

OUTPUT

PARENTS

      m_lobpcg2

CHILDREN

      blacs_gridexit,mpi_comm_free,timab

SOURCE

101   subroutine  xgScalapack_init(xgScalapack,comm,maxDim,verbosity,usable)
102 
103 
104 !This section has been created automatically by the script Abilint (TD).
105 !Do not modify the following lines by hand.
106 #undef ABI_FUNC
107 #define ABI_FUNC 'xgScalapack_init'
108 !End of the abilint section
109 
110     type(xgScalapack_t), intent(inout) :: xgScalapack
111     integer            , intent(in   ) :: comm
112     integer            , intent(in   ) :: maxDim
113     integer            , intent(in   ) :: verbosity
114     logical            , intent(  out) :: usable
115     double precision :: tsec(2)
116 #ifdef HAVE_LINALG_MKL_THREADS
117     integer :: mkl_get_max_threads
118 #endif
119     integer :: nthread
120 #ifdef HAVE_LINALG_SCALAPACK
121     integer :: maxProc
122     integer :: nproc
123     integer :: ngroup
124     integer :: subgroup
125     integer :: mycomm(2)
126     integer :: ierr
127     integer :: test_row
128     integer :: test_col
129 #else
130     ABI_UNUSED(comm)
131     ABI_UNUSED(maxDim)
132 #endif
133 
134     call timab(M__tim_init,1,tsec)
135 
136     xgScalapack%comms = xmpi_comm_null
137     xgScalapack%rank = xmpi_undefined_rank
138     xgScalapack%verbosity = verbosity
139 
140     nthread = 1
141 #ifdef HAVE_LINALG_MKL_THREADS
142     nthread =  mkl_get_max_threads()
143 #else
144     nthread = xomp_get_num_threads(open_parallel=.true.)
145     if ( nthread == 0 ) nthread = 1
146 #endif
147 
148 #ifdef HAVE_LINALG_SCALAPACK
149 
150     nproc = xmpi_comm_size(comm)
151     xgScalapack%comms(M__WORLD) = comm
152     xgScalapack%rank(M__WORLD) = xmpi_comm_rank(comm)
153     xgScalapack%size(M__WORLD) = nproc
154 
155     maxProc = (maxDim / (M__MAXDIM*nthread))+1 ! ( M__MAXDIM x M__MAXDIM matrice per MPI )
156     if ( M__CONFIG > 0 .and. M__CONFIG <= nproc ) then
157       maxProc = M__CONFIG
158     else if ( maxProc > nproc ) then
159       maxProc = nproc
160     end if
161 
162     if ( maxProc == 1 .or. M__CONFIG == SLK_DISABLED) then
163       usable = .false.
164       return
165     else if ( nthread > 1 ) then ! disable scalapack with threads since it is not threadsafe
166       ! This should be check with new elpa version en MPI+OpenMP
167       if ( M__CONFIG > 0 ) then
168         MSG_WARNING("xgScalapack turned off because you have threads")
169       end if
170       usable = .false.
171       return
172     else
173       usable = .true.
174       maxProc = 2*((maxProc+1)/2) ! Round to next even number
175     end if
176 
177     if ( xgScalapack%verbosity > 0 ) then
178       write(std_out,*) " xgScalapack will use", maxProc, "/", nproc, "MPIs"
179     end if
180 
181     ngroup = nproc/maxProc
182     xgScalapack%ngroup = ngroup
183 
184     if ( maxProc < nproc ) then
185       if ( xgScalapack%rank(M__WORLD) < maxProc*ngroup ) then
186         subgroup = xgScalapack%rank(M__WORLD)/maxProc
187         mycomm(1) = M__SLK
188         mycomm(2) = M__UNUSED
189       else
190         subgroup = ngroup+1
191         mycomm(1) = M__UNUSED
192         mycomm(2) = M__SLK
193       end if
194        call MPI_Comm_split(comm, subgroup, xgScalapack%rank(M__WORLD), xgScalapack%comms(mycomm(1)),ierr)
195        if ( ierr /= 0 ) then
196          MSG_ERROR("Error splitting communicator")
197        end if
198        xgScalapack%comms(mycomm(2)) = xmpi_comm_null
199        xgScalapack%rank(mycomm(1)) = xmpi_comm_rank(xgScalapack%comms(mycomm(1)))
200        xgScalapack%rank(mycomm(2)) = xmpi_undefined_rank
201        xgScalapack%size(mycomm(1)) = xmpi_comm_size(xgScalapack%comms(mycomm(1)))
202        xgScalapack%size(mycomm(2)) = nproc - xgScalapack%size(mycomm(1))
203     else
204        call MPI_Comm_dup(comm,xgScalapack%comms(M__SLK),ierr)
205        if ( ierr /= 0 ) then
206          MSG_ERROR("Error duplicating communicator")
207        end if
208        xgScalapack%rank(M__SLK) = xmpi_comm_rank(xgScalapack%comms(M__SLK))
209        xgScalapack%size(M__SLK) = nproc
210     end if
211 
212     if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
213       call build_grid_scalapack(xgScalapack%grid, xgScalapack%size(M__SLK), xgScalapack%comms(M__SLK))
214       call BLACS_GridInfo(xgScalapack%grid%ictxt, &
215         xgScalapack%grid%dims(M__ROW), xgScalapack%grid%dims(M__COL),&
216         xgScalapack%coords(M__ROW), xgScalapack%coords(M__COL))
217 
218      !These values are the same as those computed by BLACS_GRIDINFO
219      !except in the case where the myproc argument is not the local proc
220       test_row = INT((xgScalapack%rank(M__SLK)) / xgScalapack%grid%dims(2))
221       test_col = MOD((xgScalapack%rank(M__SLK)), xgScalapack%grid%dims(2))
222       if ( test_row /= xgScalapack%coords(M__ROW) ) then
223         MSG_WARNING("Row id mismatch")
224       end if
225       if ( test_col /= xgScalapack%coords(M__COL) ) then
226         MSG_WARNING("Col id mismatch")
227       end if
228     end if
229 
230 #else
231     usable = .false.
232 #endif
233 
234     call timab(M__tim_init,2,tsec)
235 
236   end subroutine xgScalapack_init
237 
238   subroutine xgScalapack_config(myconfig,maxDim)
239 
240 
241 !This section has been created automatically by the script Abilint (TD).
242 !Do not modify the following lines by hand.
243 #undef ABI_FUNC
244 #define ABI_FUNC 'xgScalapack_config'
245 !End of the abilint section
246 
247     integer, intent(in) :: myconfig
248     integer, intent(in) :: maxDim
249     if ( myconfig == SLK_AUTO) then
250       M__CONFIG = myconfig
251       MSG_COMMENT("xgScalapack in auto mode")
252     else if ( myconfig == SLK_DISABLED) then
253       M__CONFIG = myconfig
254       MSG_COMMENT("xgScalapack disabled")
255     else if ( myconfig > 0) then
256       M__CONFIG = myconfig
257       MSG_COMMENT("xgScalapack enabled")
258     else
259       MSG_WARNING("Bad value for xgScalapack config -> autodetection")
260       M__CONFIG = SLK_AUTO
261     end if
262     if ( maxDim > 0 ) then
263       M__MAXDIM = maxDim
264     end if
265 
266   end subroutine xgScalapack_config
267 
268   function toProcessorScalapack(xgScalapack) result(processor)
269 
270 
271 !This section has been created automatically by the script Abilint (TD).
272 !Do not modify the following lines by hand.
273 #undef ABI_FUNC
274 #define ABI_FUNC 'toProcessorScalapack'
275 !End of the abilint section
276 
277     type(xgScalapack_t), intent(in) :: xgScalapack
278     type(processor_scalapack) :: processor
279 
280     processor%myproc = xgScalapack%rank(M__SLK)
281     processor%comm = xgScalapack%comms(M__SLK)
282     processor%coords = xgScalapack%coords
283     processor%grid = xgScalapack%grid
284   end function toProcessorScalapack
285 
286   !This is for testing purpose.
287   !May not be optimal since I do not control old implementation but at least gives a reference.
288   subroutine xgScalapack_heev(xgScalapack,matrixA,eigenvalues)
289     use iso_c_binding
290 
291 !This section has been created automatically by the script Abilint (TD).
292 !Do not modify the following lines by hand.
293 #undef ABI_FUNC
294 #define ABI_FUNC 'xgScalapack_heev'
295 !End of the abilint section
296 
297     type(xgScalapack_t), intent(inout) :: xgScalapack
298     type(xgBlock_t)    , intent(inout) :: matrixA
299     type(xgBlock_t)    , intent(inout) :: eigenvalues
300 #ifdef HAVE_LINALG_SCALAPACK
301     double precision, pointer :: matrix(:,:) !(cplex*nbli_global,nbco_global)
302     double precision, pointer :: eigenvalues_tmp(:,:)
303     double precision, pointer :: vector(:)
304     double precision :: tsec(2)
305     integer :: cplex
306     integer :: istwf_k
307     integer :: nbli_global, nbco_global
308     type(c_ptr) :: cptr
309     integer :: req(2), status(MPI_STATUS_SIZE,2), ierr
310 #endif
311 
312 #ifdef HAVE_LINALG_SCALAPACK
313     call timab(M__tim_heev,1,tsec)
314 
315     ! Keep only working processors
316     if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
317 
318       call xgBlock_getSize(eigenvalues,nbli_global,nbco_global)
319       if ( cols(matrixA) /= nbli_global ) then
320         MSG_ERROR("Number of eigen values differ from number of vectors")
321       end if
322 
323       if ( space(matrixA) == SPACE_C ) then
324         cplex = 2
325         istwf_k = 1
326       else
327         cplex = 1
328         istwf_k = 2
329       endif
330 
331       call xgBlock_getSize(matrixA,nbli_global,nbco_global)
332 
333       call xgBlock_reverseMap(matrixA,matrix,nbli_global,nbco_global)
334       call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,nbco_global,1)
335       cptr = c_loc(eigenvalues_tmp)
336       call c_f_pointer(cptr,vector,(/ nbco_global /))
337 
338       call compute_eigen1(xgScalapack%comms(M__SLK), &
339         toProcessorScalapack(xgScalapack), &
340         cplex,nbli_global,nbco_global,matrix,vector,istwf_k)
341 
342     end if
343 
344     call timab(M__tim_heev,2,tsec)
345 
346     req(1:2)=-1
347     call xgScalapack_scatter(xgScalapack,matrixA,req(1))
348     call xgScalapack_scatter(xgScalapack,eigenvalues,req(2))
349 #ifdef HAVE_MPI
350     if ( any(req/=-1)  ) then
351       call MPI_WaitAll(2,req,status,ierr)
352       write(*,*) "I wait"
353       if ( ierr /= 0 ) then
354           MSG_ERROR("Error waiting data")
355       endif
356     end if
357 #endif
358 #else
359    MSG_ERROR("ScaLAPACK support not available")
360    ABI_UNUSED(xgScalapack%verbosity)
361    ABI_UNUSED(matrixA%normal)
362    ABI_UNUSED(eigenvalues%normal)
363 #endif
364 
365   end subroutine xgScalapack_heev
366 
367   !This is for testing purpose.
368   !May not be optimal since I do not control old implementation but at least gives a reference.
369   subroutine xgScalapack_hegv(xgScalapack,matrixA,matrixB,eigenvalues)
370     use iso_c_binding
371 
372 !This section has been created automatically by the script Abilint (TD).
373 !Do not modify the following lines by hand.
374 #undef ABI_FUNC
375 #define ABI_FUNC 'xgScalapack_hegv'
376 !End of the abilint section
377 
378     type(xgScalapack_t), intent(inout) :: xgScalapack
379     type(xgBlock_t)    , intent(inout) :: matrixA
380     type(xgBlock_t)    , intent(inout) :: matrixB
381     type(xgBlock_t)    , intent(inout) :: eigenvalues
382 #ifdef HAVE_LINALG_SCALAPACK
383     double precision, pointer :: matrix1(:,:) !(cplex*nbli_global,nbco_global)
384     double precision, pointer :: matrix2(:,:) !(cplex*nbli_global,nbco_global)
385     double precision, pointer :: eigenvalues_tmp(:,:)
386     double precision, pointer :: vector(:)
387     double precision :: tsec(2)
388     integer :: cplex
389     integer :: istwf_k
390     integer :: nbli_global, nbco_global
391     type(c_ptr) :: cptr
392     integer :: req(2), status(MPI_STATUS_SIZE,2),ierr
393 #endif
394 
395 #ifdef HAVE_LINALG_SCALAPACK
396     call timab(M__tim_hegv,1,tsec)
397 
398     ! Keep only working processors
399     if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
400 
401       call xgBlock_getSize(eigenvalues,nbli_global,nbco_global)
402       if ( cols(matrixA) /= cols(matrixB) ) then
403         MSG_ERROR("Matrix A and B don't have the same number of vectors")
404       end if
405 
406       if ( cols(matrixA) /= nbli_global ) then
407         MSG_ERROR("Number of eigen values differ from number of vectors")
408       end if
409 
410       if ( space(matrixA) == SPACE_C ) then
411         cplex = 2
412         istwf_k = 1
413       else
414         cplex = 1
415         istwf_k = 2
416       endif
417 
418       call xgBlock_getSize(matrixA,nbli_global,nbco_global)
419 
420       call xgBlock_reverseMap(matrixA,matrix1,nbli_global,nbco_global)
421       call xgBlock_reverseMap(matrixB,matrix2,nbli_global,nbco_global)
422       call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,nbco_global,1)
423       cptr = c_loc(eigenvalues_tmp)
424       call c_f_pointer(cptr,vector,(/ nbco_global /))
425 
426       call compute_eigen2(xgScalapack%comms(M__SLK), &
427         toProcessorScalapack(xgScalapack), &
428         cplex,nbli_global,nbco_global,matrix1,matrix2,vector,istwf_k)
429     end if
430 
431     call timab(M__tim_hegv,2,tsec)
432 
433     req(1:2)=-1
434     call xgScalapack_scatter(xgScalapack,matrixA,req(1))
435     call xgScalapack_scatter(xgScalapack,eigenvalues,req(2))
436 #ifdef HAVE_MPI
437     if ( any(req/=-1)   ) then
438       call MPI_WaitAll(2,req,status,ierr)
439       if ( ierr /= 0 ) then
440           MSG_ERROR("Error waiting data")
441       endif
442     end if
443 #endif
444 #else
445    MSG_ERROR("ScaLAPACK support not available")
446    ABI_UNUSED(xgScalapack%verbosity)
447    ABI_UNUSED(matrixA%normal)
448    ABI_UNUSED(matrixB%normal)
449    ABI_UNUSED(eigenvalues%normal)
450 #endif
451 
452   end subroutine xgScalapack_hegv
453 
454 
455   subroutine xgScalapack_scatter(xgScalapack,matrix,req)
456 
457 
458 !This section has been created automatically by the script Abilint (TD).
459 !Do not modify the following lines by hand.
460 #undef ABI_FUNC
461 #define ABI_FUNC 'xgScalapack_scatter'
462 !End of the abilint section
463 
464     type(xgScalapack_t), intent(in   ) :: xgScalapack
465     type(xgBlock_t)    , intent(inout) :: matrix
466     integer            , intent(  out) :: req
467     double precision, pointer :: tab(:,:)
468     double precision :: tsec(2)
469     integer :: cols, rows
470     integer :: ierr
471     integer :: sendto, receivefrom
472     integer :: lap
473 
474     call timab(M__tim_scatter,1,tsec)
475 
476     call xgBlock_getSize(matrix,rows,cols)
477     call xgBlock_reverseMap(matrix,tab,rows,cols)
478 
479     ! If we did the he(e|g)v and we are the first group
480     if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null .and. xgScalapack%rank(M__WORLD)<xgScalapack%size(M__SLK) ) then
481       lap = xgScalapack%ngroup
482       sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK)
483       if ( sendto < xgScalapack%size(M__WORLD) ) then
484       !do while ( sendto < xgScalapack%size(M__WORLD) )
485         !call xmpi_send(tab,sendto,sendto,xgScalapack%comms(M__WORLD),ierr)
486         call xmpi_isend(tab,sendto,sendto,xgScalapack%comms(M__WORLD),req,ierr)
487         !write(*,*) xgScalapack%rank(M__WORLD), "sends to", sendto
488         if ( ierr /= 0 ) then
489           MSG_ERROR("Error sending data")
490         end if
491         !lap = lap+1
492         !sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK)
493       !end do
494       end if
495     else if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then
496       receivefrom = MODULO(xgScalapack%rank(M__WORLD), xgScalapack%size(M__SLK))
497       if ( receivefrom >= 0 ) then
498         !call xmpi_recv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),ierr)
499         call xmpi_irecv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),req,ierr)
500         !write(*,*) xgScalapack%rank(M__WORLD), "receive from", receivefrom
501         if ( ierr /= 0 ) then
502           MSG_ERROR("Error receiving data")
503         end if
504       end if
505     !else
506       !MSG_BUG("Error scattering data")
507     end if
508 
509     call timab(M__tim_scatter,2,tsec)
510 
511   end subroutine xgScalapack_scatter
512 
513 
514   subroutine  xgScalapack_free(xgScalapack)
515 
516 
517 !This section has been created automatically by the script Abilint (TD).
518 !Do not modify the following lines by hand.
519 #undef ABI_FUNC
520 #define ABI_FUNC 'xgScalapack_free'
521 !End of the abilint section
522 
523     type(xgScalapack_t), intent(inout) :: xgScalapack
524     double precision :: tsec(2)
525 #ifdef HAVE_LINALG_SCALAPACK
526     integer :: ierr
527 #endif
528 
529     call timab(M__tim_free,1,tsec)
530 #ifdef HAVE_LINALG_SCALAPACK
531     if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then
532       call BLACS_GridExit(xgScalapack%grid%ictxt)
533       call MPI_Comm_free(xgScalapack%comms(M__SLK),ierr)
534     end if
535     if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then
536       call MPI_Comm_free(xgScalapack%comms(M__UNUSED),ierr)
537     end if
538 #else
539     ABI_UNUSED(xgScalapack%verbosity)
540 #endif
541     call timab(M__tim_free,2,tsec)
542 
543   end subroutine xgScalapack_free
544 
545 end module m_xgScalapack