TABLE OF CONTENTS
ABINIT/m_xgScalapack [ 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_profiling_abi 26 use m_xmpi 27 use m_errors 28 use m_slk 29 use m_xg 30 use m_xomp 31 32 #ifdef HAVE_MPI2 33 use mpi 34 #endif 35 36 implicit none 37 38 #ifdef HAVE_MPI1 39 include 'mpif.h' 40 #endif 41 42 43 private 44 45 integer, parameter :: M__SLK = 1 46 integer, parameter :: M__ROW = 1 47 integer, parameter :: M__COL = 2 48 integer, parameter :: M__UNUSED = 4 49 integer, parameter :: M__WORLD = 5 50 integer, parameter :: M__NDATA = 5 51 integer, parameter :: M__tim_init = 1690 52 integer, parameter :: M__tim_free = 1691 53 integer, parameter :: M__tim_heev = 1692 54 integer, parameter :: M__tim_hegv = 1693 55 integer, parameter :: M__tim_scatter = 1694 56 57 integer, parameter, public :: SLK_AUTO = -1 58 integer, parameter, public :: SLK_FORCED = 1 59 integer, parameter, public :: SLK_DISABLED = 0 60 integer, save :: M__CONFIG = SLK_AUTO 61 integer, save :: M__MAXDIM = 1000 62 63 type, public :: xgScalapack_t 64 integer :: comms(M__NDATA) 65 integer :: rank(M__NDATA) 66 integer :: size(M__NDATA) 67 integer :: coords(2) 68 integer :: ngroup 69 integer :: verbosity 70 type(grid_scalapack) :: grid 71 end type xgScalapack_t 72 73 public :: xgScalapack_init 74 public :: xgScalapack_free 75 public :: xgScalapack_heev 76 public :: xgScalapack_hegv 77 public :: xgScalapack_config 78 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
100 subroutine xgScalapack_init(xgScalapack,comm,maxDim,verbosity,usable) 101 102 103 !This section has been created automatically by the script Abilint (TD). 104 !Do not modify the following lines by hand. 105 #undef ABI_FUNC 106 #define ABI_FUNC 'xgScalapack_init' 107 use interfaces_18_timing 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 166 usable = .true. 167 maxProc = 2*((maxProc+1)/2) ! Round to next even number 168 end if 169 170 if ( xgScalapack%verbosity > 0 ) then 171 write(std_out,*) " xgScalapack will use", maxProc, "/", nproc, "MPIs" 172 end if 173 174 ngroup = nproc/maxProc 175 xgScalapack%ngroup = ngroup 176 177 if ( maxProc < nproc ) then 178 if ( xgScalapack%rank(M__WORLD) < maxProc*ngroup ) then 179 subgroup = xgScalapack%rank(M__WORLD)/maxProc 180 mycomm(1) = M__SLK 181 mycomm(2) = M__UNUSED 182 else 183 subgroup = ngroup+1 184 mycomm(1) = M__UNUSED 185 mycomm(2) = M__SLK 186 end if 187 call MPI_Comm_split(comm, subgroup, xgScalapack%rank(M__WORLD), xgScalapack%comms(mycomm(1)),ierr) 188 if ( ierr /= 0 ) then 189 MSG_ERROR("Error splitting communicator") 190 end if 191 xgScalapack%comms(mycomm(2)) = xmpi_comm_null 192 xgScalapack%rank(mycomm(1)) = xmpi_comm_rank(xgScalapack%comms(mycomm(1))) 193 xgScalapack%rank(mycomm(2)) = xmpi_undefined_rank 194 xgScalapack%size(mycomm(1)) = xmpi_comm_size(xgScalapack%comms(mycomm(1))) 195 xgScalapack%size(mycomm(2)) = nproc - xgScalapack%size(mycomm(1)) 196 else 197 call MPI_Comm_dup(comm,xgScalapack%comms(M__SLK),ierr) 198 if ( ierr /= 0 ) then 199 MSG_ERROR("Error duplicating communicator") 200 end if 201 xgScalapack%rank(M__SLK) = xmpi_comm_rank(xgScalapack%comms(M__SLK)) 202 xgScalapack%size(M__SLK) = nproc 203 end if 204 205 if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then 206 call build_grid_scalapack(xgScalapack%grid, xgScalapack%size(M__SLK), xgScalapack%comms(M__SLK)) 207 call BLACS_GridInfo(xgScalapack%grid%ictxt, & 208 xgScalapack%grid%dims(M__ROW), xgScalapack%grid%dims(M__COL),& 209 xgScalapack%coords(M__ROW), xgScalapack%coords(M__COL)) 210 211 !These values are the same as those computed by BLACS_GRIDINFO 212 !except in the case where the myproc argument is not the local proc 213 test_row = INT((xgScalapack%rank(M__SLK)) / xgScalapack%grid%dims(2)) 214 test_col = MOD((xgScalapack%rank(M__SLK)), xgScalapack%grid%dims(2)) 215 if ( test_row /= xgScalapack%coords(M__ROW) ) then 216 MSG_WARNING("Row id mismatch") 217 end if 218 if ( test_col /= xgScalapack%coords(M__COL) ) then 219 MSG_WARNING("Col id mismatch") 220 end if 221 end if 222 223 #else 224 usable = .false. 225 #endif 226 227 call timab(M__tim_init,2,tsec) 228 229 end subroutine xgScalapack_init 230 231 subroutine xgScalapack_config(myconfig,maxDim) 232 233 234 !This section has been created automatically by the script Abilint (TD). 235 !Do not modify the following lines by hand. 236 #undef ABI_FUNC 237 #define ABI_FUNC 'xgScalapack_config' 238 !End of the abilint section 239 240 integer, intent(in) :: myconfig 241 integer, intent(in) :: maxDim 242 if ( myconfig == SLK_AUTO) then 243 M__CONFIG = myconfig 244 MSG_COMMENT("xgScalapack in auto mode") 245 else if ( myconfig == SLK_DISABLED) then 246 M__CONFIG = myconfig 247 MSG_COMMENT("xgScalapack disabled") 248 else if ( myconfig > 0) then 249 M__CONFIG = myconfig 250 MSG_COMMENT("xgScalapack enabled") 251 else 252 MSG_WARNING("Bad value for xgScalapack config -> autodetection") 253 M__CONFIG = SLK_AUTO 254 end if 255 if ( maxDim > 0 ) then 256 M__MAXDIM = maxDim 257 end if 258 259 end subroutine xgScalapack_config 260 261 function toProcessorScalapack(xgScalapack) result(processor) 262 263 264 !This section has been created automatically by the script Abilint (TD). 265 !Do not modify the following lines by hand. 266 #undef ABI_FUNC 267 #define ABI_FUNC 'toProcessorScalapack' 268 !End of the abilint section 269 270 type(xgScalapack_t), intent(in) :: xgScalapack 271 type(processor_scalapack) :: processor 272 273 processor%myproc = xgScalapack%rank(M__SLK) 274 processor%comm = xgScalapack%comms(M__SLK) 275 processor%coords = xgScalapack%coords 276 processor%grid = xgScalapack%grid 277 end function toProcessorScalapack 278 279 !This is for testing purpose. 280 !May not be optimal since I do not control old implementation but at least gives a reference. 281 subroutine xgScalapack_heev(xgScalapack,matrixA,eigenvalues) 282 use iso_c_binding 283 284 !This section has been created automatically by the script Abilint (TD). 285 !Do not modify the following lines by hand. 286 #undef ABI_FUNC 287 #define ABI_FUNC 'xgScalapack_heev' 288 use interfaces_18_timing 289 !End of the abilint section 290 291 type(xgScalapack_t), intent(inout) :: xgScalapack 292 type(xgBlock_t) , intent(inout) :: matrixA 293 type(xgBlock_t) , intent(inout) :: eigenvalues 294 #ifdef HAVE_LINALG_SCALAPACK 295 double precision, pointer :: matrix(:,:) !(cplex*nbli_global,nbco_global) 296 double precision, pointer :: eigenvalues_tmp(:,:) 297 double precision, pointer :: vector(:) 298 double precision :: tsec(2) 299 integer :: cplex 300 integer :: istwf_k 301 integer :: nbli_global, nbco_global 302 type(c_ptr) :: cptr 303 integer :: req(2), status(MPI_STATUS_SIZE,2), ierr 304 #endif 305 306 #ifdef HAVE_LINALG_SCALAPACK 307 call timab(M__tim_heev,1,tsec) 308 309 ! Keep only working processors 310 if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then 311 312 call xgBlock_getSize(eigenvalues,nbli_global,nbco_global) 313 if ( cols(matrixA) /= nbli_global ) then 314 MSG_ERROR("Number of eigen values differ from number of vectors") 315 end if 316 317 if ( space(matrixA) == SPACE_C ) then 318 cplex = 2 319 istwf_k = 1 320 else 321 cplex = 1 322 istwf_k = 2 323 endif 324 325 call xgBlock_getSize(matrixA,nbli_global,nbco_global) 326 327 call xgBlock_reverseMap(matrixA,matrix,nbli_global,nbco_global) 328 call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,nbco_global,1) 329 cptr = c_loc(eigenvalues_tmp) 330 call c_f_pointer(cptr,vector,(/ nbco_global /)) 331 332 call compute_eigen1(xgScalapack%comms(M__SLK), & 333 toProcessorScalapack(xgScalapack), & 334 cplex,nbli_global,nbco_global,matrix,vector,istwf_k) 335 336 end if 337 338 call timab(M__tim_heev,2,tsec) 339 340 req(1:2)=-1 341 call xgScalapack_scatter(xgScalapack,matrixA,req(1)) 342 call xgScalapack_scatter(xgScalapack,eigenvalues,req(2)) 343 #ifdef HAVE_MPI 344 if ( any(req/=-1) ) then 345 call MPI_WaitAll(2,req,status,ierr) 346 write(*,*) "I wait" 347 if ( ierr /= 0 ) then 348 MSG_ERROR("Error waiting data") 349 endif 350 end if 351 #endif 352 #else 353 MSG_ERROR("ScaLAPACK support not available") 354 ABI_UNUSED(xgScalapack%verbosity) 355 ABI_UNUSED(matrixA%normal) 356 ABI_UNUSED(eigenvalues%normal) 357 #endif 358 359 end subroutine xgScalapack_heev 360 361 !This is for testing purpose. 362 !May not be optimal since I do not control old implementation but at least gives a reference. 363 subroutine xgScalapack_hegv(xgScalapack,matrixA,matrixB,eigenvalues) 364 use iso_c_binding 365 366 !This section has been created automatically by the script Abilint (TD). 367 !Do not modify the following lines by hand. 368 #undef ABI_FUNC 369 #define ABI_FUNC 'xgScalapack_hegv' 370 use interfaces_18_timing 371 !End of the abilint section 372 373 type(xgScalapack_t), intent(inout) :: xgScalapack 374 type(xgBlock_t) , intent(inout) :: matrixA 375 type(xgBlock_t) , intent(inout) :: matrixB 376 type(xgBlock_t) , intent(inout) :: eigenvalues 377 #ifdef HAVE_LINALG_SCALAPACK 378 double precision, pointer :: matrix1(:,:) !(cplex*nbli_global,nbco_global) 379 double precision, pointer :: matrix2(:,:) !(cplex*nbli_global,nbco_global) 380 double precision, pointer :: eigenvalues_tmp(:,:) 381 double precision, pointer :: vector(:) 382 double precision :: tsec(2) 383 integer :: cplex 384 integer :: istwf_k 385 integer :: nbli_global, nbco_global 386 type(c_ptr) :: cptr 387 integer :: req(2), status(MPI_STATUS_SIZE,2),ierr 388 #endif 389 390 #ifdef HAVE_LINALG_SCALAPACK 391 call timab(M__tim_hegv,1,tsec) 392 393 ! Keep only working processors 394 if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then 395 396 call xgBlock_getSize(eigenvalues,nbli_global,nbco_global) 397 if ( cols(matrixA) /= cols(matrixB) ) then 398 MSG_ERROR("Matrix A and B don't have the same number of vectors") 399 end if 400 401 if ( cols(matrixA) /= nbli_global ) then 402 MSG_ERROR("Number of eigen values differ from number of vectors") 403 end if 404 405 if ( space(matrixA) == SPACE_C ) then 406 cplex = 2 407 istwf_k = 1 408 else 409 cplex = 1 410 istwf_k = 2 411 endif 412 413 call xgBlock_getSize(matrixA,nbli_global,nbco_global) 414 415 call xgBlock_reverseMap(matrixA,matrix1,nbli_global,nbco_global) 416 call xgBlock_reverseMap(matrixB,matrix2,nbli_global,nbco_global) 417 call xgBlock_reverseMap(eigenvalues,eigenvalues_tmp,nbco_global,1) 418 cptr = c_loc(eigenvalues_tmp) 419 call c_f_pointer(cptr,vector,(/ nbco_global /)) 420 421 call compute_eigen2(xgScalapack%comms(M__SLK), & 422 toProcessorScalapack(xgScalapack), & 423 cplex,nbli_global,nbco_global,matrix1,matrix2,vector,istwf_k) 424 end if 425 426 call timab(M__tim_hegv,2,tsec) 427 428 req(1:2)=-1 429 call xgScalapack_scatter(xgScalapack,matrixA,req(1)) 430 call xgScalapack_scatter(xgScalapack,eigenvalues,req(2)) 431 #ifdef HAVE_MPI 432 if ( any(req/=-1) ) then 433 call MPI_WaitAll(2,req,status,ierr) 434 if ( ierr /= 0 ) then 435 MSG_ERROR("Error waiting data") 436 endif 437 end if 438 #endif 439 #else 440 MSG_ERROR("ScaLAPACK support not available") 441 ABI_UNUSED(xgScalapack%verbosity) 442 ABI_UNUSED(matrixA%normal) 443 ABI_UNUSED(matrixB%normal) 444 ABI_UNUSED(eigenvalues%normal) 445 #endif 446 447 end subroutine xgScalapack_hegv 448 449 450 subroutine xgScalapack_scatter(xgScalapack,matrix,req) 451 452 453 !This section has been created automatically by the script Abilint (TD). 454 !Do not modify the following lines by hand. 455 #undef ABI_FUNC 456 #define ABI_FUNC 'xgScalapack_scatter' 457 use interfaces_18_timing 458 !End of the abilint section 459 460 type(xgScalapack_t), intent(in ) :: xgScalapack 461 type(xgBlock_t) , intent(inout) :: matrix 462 integer , intent( out) :: req 463 double precision, pointer :: tab(:,:) 464 double precision :: tsec(2) 465 integer :: cols, rows 466 integer :: ierr 467 integer :: sendto, receivefrom 468 integer :: lap 469 470 call timab(M__tim_scatter,1,tsec) 471 472 call xgBlock_getSize(matrix,rows,cols) 473 call xgBlock_reverseMap(matrix,tab,rows,cols) 474 475 ! If we did the he(e|g)v and we are the first group 476 if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null .and. xgScalapack%rank(M__WORLD)<xgScalapack%size(M__SLK) ) then 477 lap = xgScalapack%ngroup 478 sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK) 479 if ( sendto < xgScalapack%size(M__WORLD) ) then 480 !do while ( sendto < xgScalapack%size(M__WORLD) ) 481 !call xmpi_send(tab,sendto,sendto,xgScalapack%comms(M__WORLD),ierr) 482 call xmpi_isend(tab,sendto,sendto,xgScalapack%comms(M__WORLD),req,ierr) 483 !write(*,*) xgScalapack%rank(M__WORLD), "sends to", sendto 484 if ( ierr /= 0 ) then 485 MSG_ERROR("Error sending data") 486 end if 487 !lap = lap+1 488 !sendto = xgScalapack%rank(M__WORLD) + lap*xgScalapack%size(M__SLK) 489 !end do 490 end if 491 else if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then 492 receivefrom = MODULO(xgScalapack%rank(M__WORLD), xgScalapack%size(M__SLK)) 493 if ( receivefrom >= 0 ) then 494 !call xmpi_recv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),ierr) 495 call xmpi_irecv(tab,receivefrom,xgScalapack%rank(M__WORLD),xgScalapack%comms(M__WORLD),req,ierr) 496 !write(*,*) xgScalapack%rank(M__WORLD), "receive from", receivefrom 497 if ( ierr /= 0 ) then 498 MSG_ERROR("Error receiving data") 499 end if 500 end if 501 !else 502 !MSG_BUG("Error scattering data") 503 end if 504 505 call timab(M__tim_scatter,2,tsec) 506 507 end subroutine xgScalapack_scatter 508 509 510 subroutine xgScalapack_free(xgScalapack) 511 512 513 !This section has been created automatically by the script Abilint (TD). 514 !Do not modify the following lines by hand. 515 #undef ABI_FUNC 516 #define ABI_FUNC 'xgScalapack_free' 517 use interfaces_18_timing 518 !End of the abilint section 519 520 type(xgScalapack_t), intent(inout) :: xgScalapack 521 double precision :: tsec(2) 522 #ifdef HAVE_LINALG_SCALAPACK 523 integer :: ierr 524 #endif 525 526 call timab(M__tim_free,1,tsec) 527 #ifdef HAVE_LINALG_SCALAPACK 528 if ( xgScalapack%comms(M__SLK) /= xmpi_comm_null ) then 529 call BLACS_GridExit(xgScalapack%grid%ictxt) 530 call MPI_Comm_free(xgScalapack%comms(M__SLK),ierr) 531 end if 532 if ( xgScalapack%comms(M__UNUSED) /= xmpi_comm_null ) then 533 call MPI_Comm_free(xgScalapack%comms(M__UNUSED),ierr) 534 end if 535 #else 536 ABI_UNUSED(xgScalapack%verbosity) 537 #endif 538 call timab(M__tim_free,2,tsec) 539 540 end subroutine xgScalapack_free 541 542 end module m_xgScalapack