1!================================================================================= 2! 3! Module write_matrix_m 4! 5! (1) write_matrix_d() Originally by JRD Last Modified 5/1/2008 (JRD) 6! 7! This program writes a distributed matrix like chimat or epsmat to file. 8! 9! (2) write_matrix_f() Originally by JRD Last Modified 2/5/2009 (CHP) 10! 11! Modification of write_matrix_d for full-frequency. 12! 13!================================================================================= 14 15#include "f_defs.h" 16 17module write_matrix_m 18 19 use global_m 20#ifdef HDF5 21 use hdf5 22#endif 23 use hdf5_io_m 24 use scalapack_m 25 use io_utils_m 26 use timing_m, only: timing => common_timing 27 28 implicit none 29 30 private 31 32 public :: & 33 write_matrix_d, & 34 write_matrix_d_sub, & 35 write_matrix_f 36#ifdef HDF5 37 public :: & 38 write_matrix_ser_hdf, & 39 write_matrix_f_ser_hdf, & 40 write_matrix_d_hdf, & 41 write_matrix_f_hdf, & 42 write_matrix_diagonal_hdf, & 43 write_gvec_indices_hdf, & 44 write_vcoul_hdf 45#ifdef USESCALAPACK 46 public :: & 47 write_matrix_d_par_hdf, & 48 write_matrix_d_par_hdf_sub, & 49 write_matrix_f_par_hdf 50#endif 51#endif 52 53contains 54 55!=================================================================================== 56 57subroutine write_matrix_d(scal,matrix,nmtx,iunit) 58 type(scalapack), intent(in) :: scal 59 SCALAR, intent(in) :: matrix(:,:) !< (scal%npr,scal%npc) 60 integer, intent(in) :: nmtx 61 integer, intent(in) :: iunit 62 63 integer :: ii, jj 64#ifdef USESCALAPACK 65 SCALAR, allocatable :: tempcol(:),tempcol2(:) 66 integer :: irow, icol, irowm, icolm 67 integer :: icurr 68#endif 69 type(progress_info) :: prog_info !< a user-friendly progress report 70 71 PUSH_SUB(write_matrix_d) 72 73 if (peinf%verb_debug .and. peinf%inode==0) then 74 write(6,*) 'Writing matrix: ', nmtx, iunit 75 write(6,*) 76 endif 77 78#ifdef USESCALAPACK 79 SAFE_ALLOCATE(tempcol, (nmtx)) 80 SAFE_ALLOCATE(tempcol2, (nmtx)) 81 82 icurr=0 83 84 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 85 do jj = 1, nmtx 86 call progress_step(prog_info, jj) 87! if (peinf%inode .eq. 0) then 88! write(6,*) ' In loop: ', ii 89! endif 90 icol=MOD(INT(((jj-1)/scal%nbl)+TOL_SMALL),scal%npcol) 91 tempcol=0d0 92 if (icol .eq. scal%mypcol) then 93 do ii = 1, nmtx 94 irow=MOD(INT(((ii-1)/scal%nbl)+TOL_SMALL),scal%nprow) 95 if (irow .eq. scal%myprow) then 96 icurr=icurr+1 97 icolm=INT((icurr-1)/scal%npr+TOL_SMALL)+1 98 irowm=MOD((icurr-1),scal%npr)+1 99 tempcol(ii)=matrix(irowm,icolm) 100 101! if (icolm .gt. scal%npc .or. irowm.gt.scal%npr) then 102! write(6,*) 'Error: ', scal%npc,scal%npr,icolm,irowm 103! endif 104 105 endif 106 enddo 107 endif 108 if (peinf%inode .eq. 0) then 109 tempcol2=0d0 110 endif 111 call MPI_REDUCE(tempcol,tempcol2,nmtx,MPI_SCALAR,MPI_SUM,0, & 112 MPI_COMM_WORLD,mpierr) 113 if (peinf%inode .eq. 0) then 114 write(iunit) (tempcol2(ii),ii=1,nmtx) 115 endif 116 117 call MPI_barrier(MPI_COMM_WORLD,mpierr) 118 119 enddo 120 call progress_free(prog_info) 121 122 SAFE_DEALLOCATE(tempcol) 123 SAFE_DEALLOCATE(tempcol2) 124 125! if(peinf%inode .eq. 0) then 126! write(6,*) ' Done Writing chimat: ' 127! endif 128 129#else 130 131 if (peinf%inode .eq. 0) then 132 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 133 do jj = 1, nmtx 134 call progress_step(prog_info, jj) 135 write(iunit) (matrix(ii, jj), ii = 1, nmtx) 136 enddo 137 call progress_free(prog_info) 138 endif 139 140#endif 141 142 POP_SUB(write_matrix_d) 143 144 return 145end subroutine write_matrix_d 146 147subroutine write_matrix_d_sub(scal,matrix,nmtx,iunit,neig_sub) 148 type(scalapack), intent(in) :: scal 149 complex(DPC), intent(in) :: matrix(:,:) !< (scal%npr,scal%npc) 150 integer, intent(in) :: nmtx 151 integer, intent(in) :: iunit 152 integer, intent(in), optional :: neig_sub 153 154 integer :: ii, jj, nmtx_col 155#ifdef USESCALAPACK 156 complex(DPC), allocatable :: tempcol(:),tempcol2(:) 157 integer :: irow, icol, irowm, icolm 158 integer :: icurr 159#endif 160 type(progress_info) :: prog_info !< a user-friendly progress report 161 162 PUSH_SUB(write_matrix_d_sub) 163 164 if (peinf%verb_debug .and. peinf%inode==0) then 165 write(6,*) 'Writing matrix: ', nmtx, iunit 166 write(6,*) 167 endif 168 169 ! neig_sub allows to write only neig_sub columns of the actual matrix 170 nmtx_col = nmtx 171 IF(PRESENT(neig_sub)) nmtx_col = neig_sub 172 173#ifdef USESCALAPACK 174 SAFE_ALLOCATE(tempcol, (nmtx)) 175 SAFE_ALLOCATE(tempcol2, (nmtx)) 176 177 icurr=0 178 179 call progress_init(prog_info, 'writing matrix', 'column', nmtx_col) 180 do jj = 1, nmtx_col 181 call progress_step(prog_info, jj) 182! if (peinf%inode .eq. 0) then 183! write(6,*) ' In loop: ', ii 184! endif 185 icol=MOD(INT(((jj-1)/scal%nbl)+TOL_SMALL),scal%npcol) 186 tempcol=0d0 187 if (icol .eq. scal%mypcol) then 188 do ii = 1, nmtx 189 irow=MOD(INT(((ii-1)/scal%nbl)+TOL_SMALL),scal%nprow) 190 if (irow .eq. scal%myprow) then 191 icurr=icurr+1 192 icolm=INT((icurr-1)/scal%npr+TOL_SMALL)+1 193 irowm=MOD((icurr-1),scal%npr)+1 194 tempcol(ii)=matrix(irowm,icolm) 195 196! if (icolm .gt. scal%npc .or. irowm.gt.scal%npr) then 197! write(6,*) 'Error: ', scal%npc,scal%npr,icolm,irowm 198! endif 199 200 endif 201 enddo 202 endif 203 if (peinf%inode .eq. 0) then 204 tempcol2=0d0 205 endif 206 call MPI_REDUCE(tempcol,tempcol2,nmtx,MPI_COMPLEX_DPC,MPI_SUM,0, & 207 MPI_COMM_WORLD,mpierr) 208 if (peinf%inode .eq. 0) then 209 write(iunit) (tempcol2(ii),ii=1,nmtx) 210 endif 211 212 call MPI_barrier(MPI_COMM_WORLD,mpierr) 213 214 enddo 215 call progress_free(prog_info) 216 217 if (peinf%inode .eq. 0) then 218 ! write empty rows for the missin column so sigma will work also with previously 219 ! generated epsmat (this can go in the future) 220 do jj = nmtx_col + 1, nmtx 221 write(iunit) 222 end do 223 end if 224 225 SAFE_DEALLOCATE(tempcol) 226 SAFE_DEALLOCATE(tempcol2) 227 228! if(peinf%inode .eq. 0) then 229! write(6,*) ' Done Writing chimat: ' 230! endif 231 232#else 233 234 if (peinf%inode .eq. 0) then 235 call progress_init(prog_info, 'writing matrix', 'column', nmtx_col) 236 do jj = 1, nmtx_col 237 call progress_step(prog_info, jj) 238 write(iunit) (matrix(ii, jj), ii = 1, nmtx) 239 enddo 240 call progress_free(prog_info) 241 !XXXX 242 do jj = nmtx_col + 1, nmtx 243 write(iunit) 244 end do 245 !XXXX 246 endif 247 248#endif 249 250 POP_SUB(write_matrix_d_sub) 251 252 return 253end subroutine write_matrix_d_sub 254 255!================================================================================= 256 257subroutine write_matrix_f(scal,nfreq,retarded,nmtx,iunit,nfreq_group,advanced) 258 type(scalapack), intent(in) :: scal 259 integer, intent(in) :: nfreq 260 complex(DPC), intent(in) :: retarded(:,:,:) !< (scal%npr,scal%npc,nfreq_in_group) 261 integer, intent(in) :: nmtx 262 integer, intent(in) :: iunit 263 integer, intent(in) :: nfreq_group 264 complex(DPC), optional, intent(in) :: advanced(:,:,:) !< (scal%npr,scal%npc,nfreq_in_group) 265 266 integer :: ii, jj, ifreq 267#ifdef USESCALAPACK 268 complex(DPC), allocatable :: tempcolR(:,:),tempcolR2(:,:) 269 complex(DPC), allocatable :: tempcolA(:,:),tempcolA2(:,:) 270 integer :: irow, icol, irowm, icolm,freq_grp_ind,ifreq_para 271 integer :: icurr 272#endif 273 type(progress_info) :: prog_info !< a user-friendly progress report 274 logical :: has_advanced 275 276 PUSH_SUB(write_matrix_f) 277 278 if (peinf%verb_debug .and. peinf%inode==0) then 279 write(6,*) 'Writing matrix: ', nfreq, nmtx, iunit 280 write(6,*) 281 endif 282 283 has_advanced = present(advanced) 284#ifdef USESCALAPACK 285 SAFE_ALLOCATE(tempcolR, (nfreq,nmtx)) 286 SAFE_ALLOCATE(tempcolR2, (nfreq,nmtx)) 287#ifdef CPLX 288 SAFE_ALLOCATE(tempcolA, (nfreq,nmtx)) 289 SAFE_ALLOCATE(tempcolA2, (nfreq,nmtx)) 290#endif 291 292 icurr=0 293 294 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 295 do jj = 1, nmtx 296 call progress_step(prog_info, jj) 297! if (peinf%inode .eq. 0) then 298! write(6,*) ' In loop: ', ii 299! endif 300 icol=MOD(INT(((jj-1)/scal%nbl)+TOL_SMALL),scal%npcol) 301 tempcolR=0d0 302#ifdef CPLX 303 tempcolA=0d0 304#endif 305 if (icol .eq. scal%mypcol) then 306 do ii = 1, nmtx 307 irow=MOD(INT(((ii-1)/scal%nbl)+TOL_SMALL),scal%nprow) 308 if (irow .eq. scal%myprow) then 309 icurr=icurr+1 310 icolm=INT((icurr-1)/scal%npr+TOL_SMALL)+1 311 irowm=MOD((icurr-1),scal%npr)+1 312 do ifreq=1,nfreq 313 freq_grp_ind=mod(ifreq-1,nfreq_group) 314 ifreq_para=(ifreq+nfreq_group-1)/nfreq_group 315 if(freq_grp_ind .eq. peinf%rank_mtxel) then 316 tempcolR(ifreq,ii)=retarded(irowm,icolm,ifreq_para) 317#ifdef CPLX 318 if (has_advanced) then 319 tempcolA(ifreq,ii)=advanced(irowm,icolm,ifreq_para) 320 endif 321#endif 322 endif 323 enddo 324 endif 325 enddo 326 endif 327 if (peinf%inode .eq. 0) then 328 tempcolR2=0d0 329#ifdef CPLX 330 if (has_advanced) then 331 tempcolA2=0d0 332 endif 333#endif 334 endif 335 call MPI_REDUCE(tempcolR(1,1),tempcolR2(1,1),nfreq*nmtx, & 336 MPI_COMPLEX_DPC,MPI_SUM,0,MPI_COMM_WORLD,mpierr) 337#ifdef CPLX 338 if (has_advanced) then 339 call MPI_REDUCE(tempcolA(1,1),tempcolA2(1,1),nfreq*nmtx, & 340 MPI_COMPLEX_DPC,MPI_SUM,0,MPI_COMM_WORLD,mpierr) 341 endif 342#endif 343 if (peinf%inode .eq. 0) then 344 do ii = 1, nmtx 345 write(iunit) (tempcolR2(ifreq,ii),ifreq=1,nfreq) 346 enddo 347#ifdef CPLX 348 if (has_advanced) then 349 do ii = 1, nmtx 350 write(iunit) (tempcolA2(ifreq,ii),ifreq=1,nfreq) 351 enddo 352 else 353 do ii = 1, nmtx 354 write(iunit) 355 enddo 356 endif 357#endif 358 endif 359 360 call MPI_barrier(MPI_COMM_WORLD,mpierr) 361 362 enddo 363 call progress_free(prog_info) 364 365 SAFE_DEALLOCATE(tempcolR) 366 SAFE_DEALLOCATE(tempcolR2) 367#ifdef CPLX 368 SAFE_DEALLOCATE(tempcolA) 369 SAFE_DEALLOCATE(tempcolA2) 370#endif 371 372#else 373 374 if(peinf%inode .eq. 0) then 375 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 376 do jj = 1, nmtx 377 call progress_step(prog_info, jj) 378 do ii = 1, nmtx 379 write(iunit) (retarded(ii, jj, ifreq), ifreq= 1, nfreq) 380 enddo 381#ifdef CPLX 382 if (has_advanced) then 383 do ii = 1, nmtx 384 write(iunit) (advanced(ii, jj, ifreq),ifreq = 1, nfreq) 385 enddo 386 else 387 do ii = 1, nmtx 388 write(iunit) 389 enddo 390 endif 391#endif 392 enddo 393 call progress_free(prog_info) 394 endif 395 396#endif 397 398 POP_SUB(write_matrix_f) 399 400 return 401end subroutine write_matrix_f 402 403#ifdef HDF5 404 405!======================================================================== 406! JRD: The HDF5 Equivalents of the above routines. 407!======================================================================== 408 409!========================================================================================== 410 411subroutine write_matrix_diagonal_hdf(epsdiag,nmtx,iq,isize,name) 412 real(DP), intent(in) :: epsdiag(:,:,:) !< (isize,nmtx,1) 413 integer, intent(in) :: nmtx 414 integer, intent(in) :: iq 415 integer, intent(in) :: isize 416 character(len=*), intent(in) :: name 417 418 integer(HID_T) :: file_id ! File identifier 419 integer(HID_T) :: dset_id ! Dataset identifier 420 integer(HID_T) :: filespace ! Dataspace identifier in file 421 integer(HID_T) :: memspace ! Dataspace identifier in mem 422 423 integer(HSIZE_T) :: dims(3), offset(3), offsetm(3) 424 425 integer :: error, rank 426 427 PUSH_SUB(write_matrix_diagonal_hdf) 428 429 call open_file(99, trim(name), status='old') 430 call close_file(99) 431 432 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 433 434! Write Array 435 436 rank = 3 437 dims(1) = isize 438 dims(2) = nmtx 439 dims(3) = 1 440 offset(1) = 0 441 offset(2) = 0 442 offset(3) = iq - 1 443 offsetm(:) = 0 444 445 call h5dopen_f(file_id, 'mats/matrix-diagonal', dset_id, error) 446 call h5screate_simple_f(rank, dims, memspace, error) 447 call h5dget_space_f(dset_id,filespace,error) 448 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, dims, error) 449 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, dims, error) 450 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, epsdiag, dims, error, memspace, filespace) 451 call h5dclose_f(dset_id, error) 452 call h5sclose_f(memspace, error) 453 call h5sclose_f(filespace, error) 454 455 call h5fclose_f(file_id, error) 456 457 POP_SUB(write_matrix_diagonal_hdf) 458 459end subroutine write_matrix_diagonal_hdf 460 461!=================================================================================== 462 463subroutine write_matrix_ser_hdf(matrix,nmtx,iq,is,name) 464 SCALAR, intent(in) :: matrix(:,:) !< (nmtx,nmtx) 465 integer, intent(in) :: nmtx 466 integer, intent(in) :: iq 467 integer, intent(in) :: is 468 character(len=*), intent(in) :: name 469 470 integer :: error, rank, ii, jj 471 472 integer(HID_T) :: file_id ! File identifier 473 integer(HID_T) :: dset_id ! Dataset identifier 474 integer(HID_T) :: filespace ! Dataspace identifier in file 475 integer(HID_T) :: memspace ! Dataspace identifier in mem 476 477 integer(HSIZE_T) :: count(6), offset(6), offsetm(6) 478 479 real(DP), allocatable :: data(:,:,:,:,:,:) 480 481 PUSH_SUB(write_matrix_ser_hdf) 482 483 call open_file(99, trim(name), status='old') 484 call close_file(99) 485 486 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 487 488 rank=6 489 count(1) = SCALARSIZE 490 count(2) = nmtx 491 count(3) = nmtx 492 count(4) = 1 493 count(5) = 1 494 count(6) = 1 495 496 offset(:) = 0 497 offset(5) = is - 1 498 offset(6) = iq - 1 499 500 offsetm(:) = 0 501 502 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 503 504 do jj = 1, nmtx 505 do ii = 1, nmtx 506 data(1,ii,jj,1,1,1) = dble(matrix(ii,jj)) 507#ifdef CPLX 508 data(2,ii,jj,1,1,1) = IMAG(matrix(ii,jj)) 509#endif 510 enddo 511 enddo 512 513 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 514 call h5screate_simple_f(rank, count, memspace, error) 515 call h5dget_space_f(dset_id,filespace,error) 516 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, count, error) 517 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error) 518 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, filespace) 519 520 SAFE_DEALLOCATE(data) 521 522 call h5dclose_f(dset_id, error) 523 call h5sclose_f(memspace, error) 524 call h5sclose_f(filespace, error) 525 526 call h5fclose_f(file_id, error) 527 528 POP_SUB(write_matrix_ser_hdf) 529 530end subroutine write_matrix_ser_hdf 531 532!======================================================================== 533 534subroutine write_matrix_f_ser_hdf(nfreq, retarded, nmtx, iq, is, name) 535 integer, intent(in) :: nfreq 536 complex(DPC), intent(in) :: retarded(:,:,:) !< (nmtx,nmtx,nfreq) 537 integer, intent(in) :: nmtx 538 integer, intent(in) :: iq 539 integer, intent(in) :: is 540 character(len=*), intent(in) :: name 541 542 integer :: ii, jj, error, rank 543 real(DP), allocatable :: data(:,:,:,:,:,:) 544 type(progress_info) :: prog_info !< a user-friendly progress report 545 546 integer(HID_T) :: file_id ! File identifier 547 integer(HID_T) :: dset_id ! Dataset identifier 548 integer(HID_T) :: filespace ! Dataspace identifier in file 549 integer(HID_T) :: memspace ! Dataspace identifier in mem 550 551 integer(HSIZE_T) :: count(6), offset(6), offsetm(6) 552 553 PUSH_SUB(write_matrix_f_ser_hdf) 554 555! DVF: this routine was built off of write_matrix_f_hdf to do the serial 556! writing of an hdf format matrix. This is needed for epsmat_old2hdf5.f90 557 558 rank=6 559 count(1) = 2 560 count(2) = nmtx 561 count(3) = 1 562 count(4) = nfreq 563 count(5) = 1 564 count(6) = 1 565 566 SAFE_ALLOCATE(data, (count(1),count(2),count(3),count(4),count(5),count(6))) 567 568 call open_file(99, trim(name), status='old') 569 call close_file(99) 570 571 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 572 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 573 call h5screate_simple_f(rank, count, memspace, error) 574 call h5dget_space_f(dset_id,filespace,error) 575 576 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 577! JRD XXX the fact that jj is not outer loop presents a bit of challenge 578! but this serial routine is just for legacy support 579 do jj = 1, nmtx 580 call progress_step(prog_info, jj) 581 data(1,:,1,:,1,1)=dble(retarded(:,jj,:)) 582 data(2,:,1,:,1,1)=IMAG(retarded(:,jj,:)) 583 584 offsetm(:) = 0 585 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, count, error) 586 587 offset(1)=0 588 offset(2)=0 589 offset(3)=jj-1 590 offset(4)=0 591 offset(5)=is-1 592 offset(6)=iq-1 593 594 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error) 595 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, filespace) 596 enddo 597 call progress_free(prog_info) 598 599 SAFE_DEALLOCATE(data) 600 call h5dclose_f(dset_id, error) 601 call h5sclose_f(memspace, error) 602 call h5sclose_f(filespace, error) 603 call h5fclose_f(file_id, error) 604 605 POP_SUB(write_matrix_f_ser_hdf) 606 607 return 608 609end subroutine write_matrix_f_ser_hdf 610 611!======================================================================== 612 613subroutine write_matrix_d_hdf(scal,matrix,nmtx,iq,is,name) 614 type(scalapack), intent(in) :: scal 615 SCALAR, intent(in) :: matrix(:,:) !< (scal%npr,scal%npc) 616 integer, intent(in) :: nmtx 617 integer, intent(in) :: iq 618 integer, intent(in) :: is 619 character(len=*), intent(in) :: name 620 621 integer :: ii, jj, error, size, rank 622#ifdef USESCALAPACK 623 real(DP), allocatable :: datatmp(:,:,:,:,:,:) 624 integer :: irow, icol, irowm, icolm 625 integer :: icurr 626#endif 627 real(DP), allocatable :: data(:,:,:,:,:,:) 628 type(progress_info) :: prog_info !< a user-friendly progress report 629 630 integer(HID_T) :: file_id ! File identifier 631 integer(HID_T) :: dset_id ! Dataset identifier 632 integer(HID_T) :: filespace ! Dataspace identifier in file 633 integer(HID_T) :: memspace ! Dataspace identifier in mem 634#ifdef USESCALAPACK 635! integer(HID_T) :: plist_id ! Property list identifier for parallel IO 636! ! Not used yet... 637#endif 638 639 integer(HSIZE_T) :: count(6), offset(6), offsetm(6) 640 641 PUSH_SUB(write_matrix_d_hdf) 642 643 if (peinf%verb_debug .and. peinf%inode==0) then 644 write(6,*) 'Writing matrix: ', nmtx 645 write(6,*) 646 endif 647 648! XXX: For now, we will still have only proc 0 write... 649! We should changes this to parallel writes. But doing 650! this effectively from the scalapack, block cyclic layout 651! seems a bit tricky. So, ignoring for now... 652 653 rank=6 654 count(1) = SCALARSIZE 655 count(2) = nmtx 656 count(3) = 1 657 count(4) = 1 658 count(5) = 1 659 count(6) = 1 660 661 if (peinf%inode .eq. 0) then 662 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 663 664 call open_file(99, trim(name), status='old') 665 call close_file(99) 666 667 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 668 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 669 call h5screate_simple_f(rank, count, memspace, error) 670 call h5dget_space_f(dset_id,filespace,error) 671 endif 672 673#ifdef USESCALAPACK 674 SAFE_ALLOCATE(datatmp, (count(1),count(2),count(3),count(4),count(5),count(6))) 675 icurr=0 676#endif 677 678 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 679 do jj = 1, nmtx 680 681 call progress_step(prog_info, jj) 682#ifdef USESCALAPACK 683 684 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_comm) 685 686 icol=MOD(INT(((jj-1)/scal%nbl)+TOL_SMALL),scal%npcol) 687 datatmp=0d0 688 if (icol .eq. scal%mypcol) then 689 do ii = 1, nmtx 690 irow=MOD(INT(((ii-1)/scal%nbl)+TOL_SMALL),scal%nprow) 691 if (irow .eq. scal%myprow) then 692 icurr=icurr+1 693 icolm=INT((icurr-1)/scal%npr+TOL_SMALL)+1 694 irowm=MOD((icurr-1),scal%npr)+1 695 datatmp(1,ii,1,1,1,1)=dble(matrix(irowm,icolm)) 696#ifdef CPLX 697 datatmp(2,ii,1,1,1,1)=IMAG(matrix(irowm,icolm)) 698#endif 699 endif 700 enddo 701 endif 702 if (peinf%inode .eq. 0) then 703 data=0d0 704 endif 705 706! XXX This is a big waste of communication. Should be fixed when do 707! parallel IO. 708 709 size = nmtx * SCALARSIZE 710 711 call MPI_REDUCE(datatmp,data,size,MPI_DOUBLE_PRECISION,MPI_SUM,0, & 712 MPI_COMM_WORLD,mpierr) 713 714 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_comm) 715 716#else 717 718 if (peinf%inode .eq. 0) then 719 do ii = 1, nmtx 720 data(1,ii,1,1,1,1) = dble(matrix(ii,jj)) 721#ifdef CPLX 722 data(2,ii,1,1,1,1) = IMAG(matrix(ii,jj)) 723#endif 724 enddo 725 endif 726 727#endif 728 729 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_io) 730 731 if (peinf%inode .eq. 0) then 732 733 offsetm(:) = 0 734 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, count, error) 735 736 offset(1)=0 737 offset(2)=0 738 offset(3)=jj-1 739 offset(4)=0 740 offset(5)=is-1 741 offset(6)=iq-1 742 743 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error) 744 745 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, filespace) 746 747 endif 748 749 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_io) 750 751#ifdef USESCALAPACK 752 call MPI_barrier(MPI_COMM_WORLD,mpierr) 753#endif 754 755 enddo 756 call progress_free(prog_info) 757 758#ifdef USESCALAPACK 759 SAFE_DEALLOCATE(datatmp) 760#endif 761 762 if (peinf%inode .eq. 0) then 763 SAFE_DEALLOCATE(data) 764 call h5dclose_f(dset_id, error) 765 call h5sclose_f(memspace, error) 766 call h5sclose_f(filespace, error) 767 call h5fclose_f(file_id, error) 768 endif 769 770 POP_SUB(write_matrix_d_hdf) 771 772 return 773end subroutine write_matrix_d_hdf 774 775!======================================================================== 776 777#ifdef USESCALAPACK 778 779subroutine write_matrix_d_par_hdf(scal,matrix,nmtx,iq,is,name) 780 type(scalapack), intent(in) :: scal 781 SCALAR, intent(in) :: matrix(:,:) !< (scal%npr,scal%npc) 782 integer, intent(in) :: nmtx 783 integer, intent(in) :: iq 784 integer, intent(in) :: is 785 character(len=*), intent(in) :: name 786 787 integer :: ii, jj, error, rank 788 real(DP), allocatable :: data(:,:,:,:,:,:) 789 790 integer(HID_T) :: file_id ! File identifier 791 integer(HID_T) :: dset_id ! Dataset identifier 792 integer(HID_T) :: filespace ! Dataspace identifier in file 793 integer(HID_T) :: memspace ! Dataspace identifier in mem 794 integer(HID_T) :: plist_id ! Property list identifier for parallel IO 795 796 integer(HSIZE_T) :: count(6), countm(6), offset(6), offsetm(6), stride(6), block(6) 797 integer(HSIZE_T) :: countr(6), offsetr(6), strider(6), blockr(6) 798 799 integer :: comm, info, rowremainder, colremainder 800 801 PUSH_SUB(write_matrix_d_par_hdf) 802 803 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_comm) 804 805! JRD: We need a barrier here or else parallel file opening gets mixed up with 806! peinf%inode 0 opening the file to write the diagonal (which is called first). 807 call MPI_barrier(MPI_COMM_WORLD,mpierr) 808 809 comm = MPI_COMM_WORLD 810 info = MPI_INFO_NULL 811 812 if (peinf%verb_debug .and. peinf%inode==0) then 813 write(6,*) 'Writing matrix: ', nmtx 814 write(6,*) 815 endif 816 817! JRD Should be ok with npr and npc = 0 818 819 !if (scal%npr .eq. 0 .or. scal%npc .eq. 0) then 820 ! write(6,*) peinf%inode,"Zero npr or npc!!", scal%npr, scal%npc 821 !endif 822 823 rank=6 824 countm(1) = SCALARSIZE 825 countm(2) = scal%npr 826 countm(3) = scal%npc 827 countm(4) = 1 828 countm(5) = 1 829 countm(6) = 1 830 831 offsetm(:) = 0 832 833 count(1) = 1 834 count(2) = scal%npr/scal%nbl 835 count(3) = scal%npc/scal%nbl 836 count(4) = 1 837 count(5) = 1 838 count(6) = 1 839 840 block(1) = SCALARSIZE 841 block(2) = scal%nbl 842 block(3) = scal%nbl 843 block(4) = 1 844 block(5) = 1 845 block(6) = 1 846 847 offset(1) = 0 848 offset(2) = scal%myprow*scal%nbl 849 offset(3) = scal%mypcol*scal%nbl 850 offset(4) = 0 851 offset(5) = is-1 852 offset(6) = iq-1 853 854 stride(1) = 1 855 stride(2) = scal%nprow*scal%nbl 856 stride(3) = scal%npcol*scal%nbl 857 stride(4) = 1 858 stride(5) = 1 859 stride(6) = 1 860 861 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 862 call h5pset_fapl_mpio_f(plist_id, comm, info, error) 863 864 call open_file(99, trim(name), status='old') 865 call close_file(99) 866 867 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error, access_prp = plist_id) 868 call h5pclose_f(plist_id,error) 869 870 SAFE_ALLOCATE(data,(countm(1),countm(2),countm(3),countm(4),countm(5),countm(6))) 871!XXX create data can we avoid duplication? 872 do jj = 1, scal%npc 873 !do jj = 1, scal%npc - mod(scal%npc,scal%nbl) 874 do ii = 1, scal%npr 875 !do ii = 1, scal%npr - mod(scal%npr,scal%nbl) 876 data(1,ii,jj,1,1,1) = dble(matrix(ii,jj)) 877#ifdef CPLX 878 data(2,ii,jj,1,1,1) = IMAG(matrix(ii,jj)) 879#endif 880 enddo 881 enddo 882 883 call h5screate_simple_f(rank, countm, memspace, error) 884 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, countm, error) 885 886 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 887 call h5dget_space_f(dset_id,filespace,error) 888 889 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error, stride, block) 890 891! Add in remainders, in case scal%nbl doesnt perfectly divide nmtx 892 893! Bottom Rows 894 rowremainder = mod(scal%npr,scal%nbl) 895 if (rowremainder .ne. 0) then 896 offsetr=offset 897 countr=count 898 blockr=block 899 strider=stride 900 offsetr(2)=nmtx-rowremainder 901 countr(2)=rowremainder 902 blockr(2)=1 903 strider(2)=1 904 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 905 !write(6,*) peinf%inode, "I have the bottom row", rowremainder, scal%npc 906 endif 907 908! Right Columns 909 colremainder = mod(scal%npc,scal%nbl) 910 if (colremainder .ne. 0) then 911 offsetr=offset 912 countr=count 913 blockr=block 914 strider=stride 915 offsetr(3)=nmtx-colremainder 916 countr(3)=colremainder 917 blockr(3)=1 918 strider(3)=1 919 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 920 !write(6,*) peinf%inode, "I have the right column", colremainder, scal%npr 921! Bottom Corner of Matrix 922 if (rowremainder .ne. 0) then 923 offsetr=offset 924 countr=count 925 blockr=block 926 strider=stride 927 offsetr(2)=nmtx-rowremainder 928 countr(2)=rowremainder 929 blockr(2)=1 930 strider(2)=1 931 offsetr(3)=nmtx-colremainder 932 countr(3)=colremainder 933 blockr(3)=1 934 strider(3)=1 935 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 936 !write(6,*) peinf%inode, "I have bottom both" 937 endif 938 endif 939 940 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 941 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 942 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_comm) 943 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_io) 944 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, filespace, & 945 xfer_prp = plist_id) 946 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_io) 947 call h5pclose_f(plist_id, error) 948 949 SAFE_DEALLOCATE(data) 950 call h5dclose_f(dset_id, error) 951 call h5sclose_f(memspace, error) 952 call h5sclose_f(filespace, error) 953 call h5fclose_f(file_id, error) 954 955 POP_SUB(write_matrix_d_par_hdf) 956 957 return 958end subroutine write_matrix_d_par_hdf 959 960!======================================================================================== 961! subspace routine for eigenvectors and full epsilon zero 962subroutine write_matrix_d_par_hdf_sub(scal, matrix, nmtx, nmtx_col, iq, is, name, set_name) 963 type(scalapack), intent(in) :: scal 964 complex(DPC), intent(in) :: matrix(:,:) !< (scal%npr,scal%npc) 965 integer, intent(in) :: nmtx, nmtx_col 966 integer, intent(in) :: iq 967 integer, intent(in) :: is 968 character(len=*), intent(in) :: name 969 character(len=*), intent(in) :: set_name 970 971 integer :: ii, jj, error, rank 972 real(DP), allocatable :: data(:,:,:,:,:,:) 973 974 integer(HID_T) :: file_id ! File identifier 975 integer(HID_T) :: dset_id ! Dataset identifier 976 integer(HID_T) :: filespace ! Dataspace identifier in file 977 integer(HID_T) :: memspace ! Dataspace identifier in mem 978 integer(HID_T) :: plist_id ! Property list identifier for parallel IO 979 980 integer(HSIZE_T) :: count(6), countm(6), offset(6), offsetm(6), stride(6), block(6) 981 integer(HSIZE_T) :: countr(6), offsetr(6), strider(6), blockr(6) 982 983 integer :: comm, info, rowremainder, colremainder 984 985 PUSH_SUB(write_matrix_d_par_hdf_sub) 986 987 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_comm) 988 989! JRD: We need a barrier here or else parallel file opening gets mixed up with 990! peinf%inode 0 opening the file to write the diagonal (which is called first). 991 call MPI_barrier(MPI_COMM_WORLD,mpierr) 992 993 comm = MPI_COMM_WORLD 994 info = MPI_INFO_NULL 995 996 if (peinf%verb_debug .and. peinf%inode==0) then 997 write(6,*) 'Writing matrix: ', nmtx, 'x', nmtx_col 998 write(6,*) 999 endif 1000 1001! JRD Should be ok with npr and npc = 0 1002 1003 !if (scal%npr .eq. 0 .or. scal%npc .eq. 0) then 1004 ! write(6,*) peinf%inode,"Zero npr or npc!!", scal%npr, scal%npc 1005 !endif 1006 1007 rank=6 1008 countm(1) = 2 1009 countm(2) = scal%npr 1010 countm(3) = scal%npc 1011 countm(4) = 1 1012 countm(5) = 1 1013 countm(6) = 1 1014 1015 offsetm(:) = 0 1016 1017 count(1) = 1 1018 count(2) = scal%npr/scal%nbl 1019 count(3) = scal%npc/scal%nbl 1020 count(4) = 1 1021 count(5) = 1 1022 count(6) = 1 1023 1024 block(1) = 2 1025 block(2) = scal%nbl 1026 block(3) = scal%nbl 1027 block(4) = 1 1028 block(5) = 1 1029 block(6) = 1 1030 1031 offset(1) = 0 1032 offset(2) = scal%myprow*scal%nbl 1033 offset(3) = scal%mypcol*scal%nbl 1034 offset(4) = 0 1035 offset(5) = is-1 1036 offset(6) = iq-1 1037 1038 stride(1) = 1 1039 stride(2) = scal%nprow*scal%nbl 1040 stride(3) = scal%npcol*scal%nbl 1041 stride(4) = 1 1042 stride(5) = 1 1043 stride(6) = 1 1044 1045 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 1046 call h5pset_fapl_mpio_f(plist_id, comm, info, error) 1047 1048 call open_file(99, trim(name), status='old') 1049 call close_file(99) 1050 1051 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error, access_prp = plist_id) 1052 call h5pclose_f(plist_id,error) 1053 1054 SAFE_ALLOCATE(data,(countm(1),countm(2),countm(3),countm(4),countm(5),countm(6))) 1055!XXX create data can we avoid duplication? 1056 do jj = 1, scal%npc 1057 !do jj = 1, scal%npc - mod(scal%npc,scal%nbl) 1058 do ii = 1, scal%npr 1059 !do ii = 1, scal%npr - mod(scal%npr,scal%nbl) 1060 data(1,ii,jj,1,1,1) = dble(matrix(ii,jj)) 1061 data(2,ii,jj,1,1,1) = IMAG(matrix(ii,jj)) 1062 enddo 1063 enddo 1064 1065 call h5screate_simple_f(rank, countm, memspace, error) 1066 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, countm, error) 1067 1068 call h5dopen_f(file_id, trim(set_name), dset_id, error) 1069 call h5dget_space_f(dset_id,filespace,error) 1070 1071 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error, stride, block) 1072 1073! Add in remainders, in case scal%nbl doesnt perfectly divide nmtx 1074 1075! Bottom Rows 1076 rowremainder = mod(scal%npr,scal%nbl) 1077 if (rowremainder .ne. 0) then 1078 offsetr=offset 1079 countr=count 1080 blockr=block 1081 strider=stride 1082 offsetr(2)=nmtx-rowremainder 1083 countr(2)=rowremainder 1084 blockr(2)=1 1085 strider(2)=1 1086 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1087 !write(6,*) peinf%inode, "I have the bottom row", rowremainder, scal%npc 1088 endif 1089 1090! Right Columns 1091 colremainder = mod(scal%npc,scal%nbl) 1092 if (colremainder .ne. 0) then 1093 offsetr=offset 1094 countr=count 1095 blockr=block 1096 strider=stride 1097 offsetr(3)=nmtx_col - colremainder 1098 countr(3)=colremainder 1099 blockr(3)=1 1100 strider(3)=1 1101 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1102 !write(6,*) peinf%inode, "I have the right column", colremainder, scal%npr 1103! Bottom Corner of Matrix 1104 if (rowremainder .ne. 0) then 1105 offsetr=offset 1106 countr=count 1107 blockr=block 1108 strider=stride 1109 offsetr(2)=nmtx-rowremainder 1110 countr(2)=rowremainder 1111 blockr(2)=1 1112 strider(2)=1 1113 offsetr(3)=nmtx_col-colremainder 1114 countr(3)=colremainder 1115 blockr(3)=1 1116 strider(3)=1 1117 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1118 !write(6,*) peinf%inode, "I have bottom both" 1119 endif 1120 endif 1121 1122 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 1123 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 1124 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_comm) 1125 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_io) 1126 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, filespace, & 1127 xfer_prp = plist_id) 1128 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_io) 1129 call h5pclose_f(plist_id, error) 1130 1131 SAFE_DEALLOCATE(data) 1132 call h5dclose_f(dset_id, error) 1133 call h5sclose_f(memspace, error) 1134 call h5sclose_f(filespace, error) 1135 call h5fclose_f(file_id, error) 1136 1137 POP_SUB(write_matrix_d_par_hdf_sub) 1138 1139 return 1140end subroutine write_matrix_d_par_hdf_sub 1141 1142!======================================================================================== 1143 1144subroutine write_matrix_f_par_hdf(scal, nfreq_in_group, retarded, nmtx, iq, is, name, nfreq_group,& 1145 set_name) 1146 type(scalapack), intent(in) :: scal 1147 integer, intent(in) :: nfreq_in_group 1148 complex(DPC), intent(in) :: retarded(:,:,:) !< (scal%npr,scal%npc,nfreq_in_group) 1149 integer, intent(in) :: nmtx 1150 integer, intent(in) :: iq 1151 integer, intent(in) :: is 1152 character(len=*), intent(in) :: name 1153 integer, intent(in) :: nfreq_group 1154 character(len=*), intent(in), optional :: set_name 1155 1156 integer :: ii, jj, error, rank 1157 real(DP), allocatable :: data(:,:,:,:,:,:) 1158 1159 integer(HID_T) :: file_id ! File identifier 1160 integer(HID_T) :: dset_id ! Dataset identifier 1161 integer(HID_T) :: filespace ! Dataspace identifier in file 1162 integer(HID_T) :: memspace ! Dataspace identifier in mem 1163 integer(HID_T) :: plist_id ! Property list identifier for parallel IO 1164 1165 integer(HSIZE_T) :: count(6), countm(6), offset(6), offsetm(6), stride(6), block(6) 1166 integer(HSIZE_T) :: countr(6), offsetr(6), strider(6), blockr(6) 1167 1168 integer :: comm, info, rowremainder, colremainder 1169 1170 logical :: default_set_name 1171 1172 PUSH_SUB(write_matrix_f_par_hdf) 1173 1174 default_set_name = .true. 1175 if( present(set_name) ) then 1176 default_set_name = .false. 1177 end if 1178 1179 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_comm) 1180 1181! JRD: We need a barrier here or else parallel file opening gets mixed up with 1182! peinf%inode 0 opening the file to write the diagonal (which is called first). 1183 call MPI_barrier(MPI_COMM_WORLD,mpierr) 1184 1185 comm = MPI_COMM_WORLD 1186 info = MPI_INFO_NULL 1187 1188 if (peinf%verb_debug .and. peinf%inode==0) then 1189 write(6,*) 'Writing matrix: ', nmtx 1190 write(6,*) 1191 endif 1192 1193! JRD Should be ok with npr and npc = 0 1194 1195 rank=6 1196 countm(1) = 2 1197 countm(2)=max(1,scal%npr) 1198 countm(3)=max(1,scal%npc) 1199 countm(4)=max(1,nfreq_in_group) 1200 countm(5) = 1 1201 countm(6) = 1 1202 1203 offsetm(:) = 0 1204 1205 count(1) = 1 1206 count(2) = scal%npr/scal%nbl 1207 count(3) = scal%npc/scal%nbl 1208 if(nfreq_group .gt. 1) then 1209 count(4) = nfreq_in_group 1210 else 1211 count(4) = 1 1212 endif 1213 count(5) = 1 1214 count(6) = 1 1215 1216 block(1) = 2 1217 block(2) = scal%nbl 1218 block(3) = scal%nbl 1219 if(nfreq_group .gt. 1) then 1220 block(4) = 1 1221 else 1222 block(4) = nfreq_in_group 1223 endif 1224 block(5) = 1 1225 block(6) = 1 1226 1227 offset(1) = 0 1228 offset(2) = scal%myprow*scal%nbl 1229 offset(3) = scal%mypcol*scal%nbl 1230 offset(4) = peinf%igroup_f 1231 offset(5) = is-1 1232 offset(6) = iq-1 1233 1234 stride(1) = 1 1235 stride(2) = scal%nprow*scal%nbl 1236 stride(3) = scal%npcol*scal%nbl 1237 stride(4) = nfreq_group 1238 stride(5) = 1 1239 stride(6) = 1 1240 1241 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 1242 call h5pset_fapl_mpio_f(plist_id, comm, info, error) 1243 1244 call open_file(99, trim(name), status='old') 1245 call close_file(99) 1246 1247 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error, access_prp = plist_id) 1248 call h5pclose_f(plist_id,error) 1249 1250 SAFE_ALLOCATE(data,(countm(1),countm(2),countm(3),countm(4),countm(5),countm(6))) 1251!XXX create data can we avoid duplication? 1252!XXX THREAD? Yes we should to get better bandwidth 1253 if (scal%npr/=0 .and. scal%npc/=0 .and. nfreq_in_group/=0) then 1254 data(1,:,:,:,1,1) = dble(retarded(:,:,:)) 1255 data(2,:,:,:,1,1) = IMAG(retarded(:,:,:)) 1256 endif 1257 1258 call h5screate_simple_f(rank, countm, memspace, error) 1259 1260 if(scal%npr*scal%npc .ne. 0) then 1261 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, countm, error) 1262 else 1263 call H5sselect_none_f(memspace,error) 1264 endif 1265 1266 if( default_set_name ) then 1267 ! use default set location 1268 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 1269 else 1270 ! use set location from input 1271 call h5dopen_f(file_id, trim(set_name), dset_id, error) 1272 end if 1273 call h5dget_space_f(dset_id,filespace,error) 1274 1275 if(scal%npr*scal%npc .ne. 0) then 1276 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error, stride, block) 1277 else 1278 call H5sselect_none_f(filespace,error) 1279 endif 1280 1281! Add in remainders, in case scal%nbl doesnt perfectly divide nmtx 1282 1283! The condition that nfreq_in_group .ne. 0 is here because this 1284! condition will only be met by the processors that are excluded 1285! from the calculation when using parallel frequencies with a 1286! number of processors not divisible by the number of frequency 1287! groups (for the paranoid: this condition that some processors don`t own any frequencies 1288! could also be met if the number of frequency groups 1289! requested by the user is greater than the total number of frequencies calculated, 1290! but we don`t allow that, i.e. the code dies if the user makes such a request). 1291! They will have a row/colremainder = 0 because they own no part of the dielectric 1292! matrix, but we don`t want them to be involved with the hyperslab operations because 1293! they have no data and are spectators. 1294 1295! Bottom Rows 1296 rowremainder = mod(scal%npr,scal%nbl) 1297 if (rowremainder .ne. 0 .and. nfreq_in_group .ne. 0) then 1298 offsetr=offset 1299 countr=count 1300 blockr=block 1301 strider=stride 1302 offsetr(2)=nmtx-rowremainder 1303 countr(2)=rowremainder 1304 blockr(2)=1 1305 strider(2)=1 1306 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1307 !write(6,*) peinf%inode, "I have the bottom row", rowremainder, scal%npc 1308 endif 1309 1310! Right Columns 1311 colremainder = mod(scal%npc,scal%nbl) 1312 if (colremainder .ne. 0 .and. nfreq_in_group .ne. 0) then 1313 offsetr=offset 1314 countr=count 1315 blockr=block 1316 strider=stride 1317 offsetr(3)=nmtx-colremainder 1318 countr(3)=colremainder 1319 blockr(3)=1 1320 strider(3)=1 1321 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1322 !write(6,*) peinf%inode, "I have the right column", colremainder, scal%npr 1323! Bottom Corner of Matrix 1324 if (rowremainder .ne. 0) then 1325 offsetr=offset 1326 countr=count 1327 blockr=block 1328 strider=stride 1329 offsetr(2)=nmtx-rowremainder 1330 countr(2)=rowremainder 1331 blockr(2)=1 1332 strider(2)=1 1333 offsetr(3)=nmtx-colremainder 1334 countr(3)=colremainder 1335 blockr(3)=1 1336 strider(3)=1 1337 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1338 !write(6,*) peinf%inode, "I have bottom both" 1339 endif 1340 endif 1341 1342 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 1343 !call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, error) 1344 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 1345 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_comm) 1346 if ( peinf%inode == 0 ) call timing%start(timing%eps_i_o_io) 1347 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, filespace, & 1348 xfer_prp = plist_id) 1349 if ( peinf%inode == 0 ) call timing%stop(timing%eps_i_o_io) 1350 call h5pclose_f(plist_id, error) 1351 1352 SAFE_DEALLOCATE(data) 1353 call h5dclose_f(dset_id, error) 1354 call h5sclose_f(memspace, error) 1355 call h5sclose_f(filespace, error) 1356 call h5fclose_f(file_id, error) 1357 1358 POP_SUB(write_matrix_f_par_hdf) 1359 1360 return 1361 1362end subroutine write_matrix_f_par_hdf 1363 1364#endif 1365 1366 1367!========================================================================================== 1368 1369subroutine write_matrix_f_hdf(scal, nfreq, retarded, nmtx, iq, is, name) 1370 type(scalapack), intent(in) :: scal 1371 integer, intent(in) :: nfreq 1372 complex(DPC), intent(in) :: retarded(:,:,:) !< (scal%npr,scal%npc,nfreq) 1373 integer, intent(in) :: nmtx 1374 integer, intent(in) :: iq 1375 integer, intent(in) :: is 1376 character(len=*), intent(in) :: name 1377 1378 integer :: ii, jj, error, size, rank 1379#ifdef USESCALAPACK 1380 real(DP), allocatable :: datatmp(:,:,:,:,:,:) 1381 integer :: irow, icol, irowm, icolm 1382 integer :: icurr 1383#endif 1384 real(DP), allocatable :: data(:,:,:,:,:,:) 1385 type(progress_info) :: prog_info !< a user-friendly progress report 1386 1387 integer(HID_T) :: file_id ! File identifier 1388 integer(HID_T) :: dset_id ! Dataset identifier 1389 integer(HID_T) :: filespace ! Dataspace identifier in file 1390 integer(HID_T) :: memspace ! Dataspace identifier in mem 1391 1392 integer(HSIZE_T) :: count(6), offset(6), offsetm(6) 1393 1394 PUSH_SUB(write_matrix_f_hdf) 1395 1396 if (peinf%verb_debug .and. peinf%inode==0) then 1397 write(6,*) 'Writing matrix: ', nmtx, nfreq 1398 write(6,*) 1399 endif 1400 1401! XXX: For now, we will still have only proc 0 write... 1402! We should changes this to parallel writes. But doing 1403! this effectively from the scalapack, block cyclic layout 1404! seems a bit tricky. So, ignoring for now... 1405 1406 rank=6 1407 count(1) = 2 1408 count(2) = nmtx 1409 count(3) = 1 1410 count(4) = nfreq 1411 count(5) = 1 1412 count(6) = 1 1413 1414 if (peinf%inode .eq. 0) then 1415 SAFE_ALLOCATE(data, (count(1),count(2),count(3),count(4),count(5),count(6))) 1416 1417 call open_file(99, trim(name), status='old') 1418 call close_file(99) 1419 1420 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 1421 call h5dopen_f(file_id, 'mats/matrix', dset_id, error) 1422 call h5screate_simple_f(rank, count, memspace, error) 1423 call h5dget_space_f(dset_id,filespace,error) 1424 endif 1425 1426#ifdef USESCALAPACK 1427 SAFE_ALLOCATE(datatmp, (count(1),count(2),count(3),count(4),count(5),count(6))) 1428 icurr=0 1429#endif 1430 1431 call progress_init(prog_info, 'writing matrix', 'column', nmtx) 1432 do jj = 1, nmtx 1433 call progress_step(prog_info, jj) 1434 1435#ifdef USESCALAPACK 1436 1437 icol=MOD(INT(((jj-1)/scal%nbl)+TOL_SMALL),scal%npcol) 1438 datatmp=0d0 1439! JRD XXX The below is going to be incredibly slow. Hoping all over memory. 1440 if (icol .eq. scal%mypcol) then 1441 do ii = 1, nmtx 1442 irow=MOD(INT(((ii-1)/scal%nbl)+TOL_SMALL),scal%nprow) 1443 if (irow .eq. scal%myprow) then 1444 icurr=icurr+1 1445 icolm=INT((icurr-1)/scal%npr+TOL_SMALL)+1 1446 irowm=MOD((icurr-1),scal%npr)+1 1447 datatmp(1,ii,1,:,1,1)=dble(retarded(irowm,icolm,:)) 1448 datatmp(2,ii,1,:,1,1)=IMAG(retarded(irowm,icolm,:)) 1449 endif 1450 enddo 1451 endif 1452 if (peinf%inode .eq. 0) then 1453 data=0d0 1454 endif 1455! XXX This is a big waste of communication. Should be fixed when do 1456! parallel IO. 1457 1458 size = nmtx*nfreq*2 1459 1460 call MPI_REDUCE(datatmp,data,size,MPI_DOUBLE_PRECISION,MPI_SUM,0, & 1461 MPI_COMM_WORLD,mpierr) 1462 1463#else 1464! JRD XXX The below is horrendously slow. Jumping all over memory 1465 if (peinf%inode .eq. 0) then 1466 do ii = 1, nmtx 1467 data(1,ii,1,:,1,1)=dble(retarded(ii,jj,:)) 1468 data(2,ii,1,:,1,1)=IMAG(retarded(ii,jj,:)) 1469 enddo 1470 endif 1471 1472#endif 1473 1474 if (peinf%inode .eq. 0) then 1475 1476 offsetm(:) = 0 1477 call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, count, error) 1478 1479 offset(1)=0 1480 offset(2)=0 1481 offset(3)=jj-1 1482 offset(4)=0 1483 offset(5)=is-1 1484 offset(6)=iq-1 1485 1486 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error) 1487 1488 call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, filespace) 1489 1490 endif 1491 1492#ifdef USESCALAPACK 1493 call MPI_barrier(MPI_COMM_WORLD,mpierr) 1494#endif 1495 1496 1497 enddo 1498 call progress_free(prog_info) 1499 1500#ifdef USESCALAPACK 1501 SAFE_DEALLOCATE(datatmp) 1502#endif 1503 1504 if (peinf%inode .eq. 0) then 1505 SAFE_DEALLOCATE(data) 1506 endif 1507 1508 if (peinf%inode .eq. 0) then 1509 call h5dclose_f(dset_id, error) 1510 call h5sclose_f(memspace, error) 1511 call h5sclose_f(filespace, error) 1512 call h5fclose_f(file_id, error) 1513 endif 1514 1515 POP_SUB(write_matrix_f_hdf) 1516 1517 return 1518end subroutine write_matrix_f_hdf 1519 1520!=================================================================================== 1521 1522subroutine write_gvec_indices_hdf(ng, gind_eps2rho, gind_rho2eps, ekin, iq, name) 1523 integer, intent(in) :: ng 1524 integer, intent(in) :: gind_eps2rho(:) !< (ng) 1525 integer, intent(in) :: gind_rho2eps(:) !< (ng) 1526 real(DP), intent(in) :: ekin(:) !< (ng) 1527 integer, intent(in) :: iq 1528 character(len=*), intent(in) :: name 1529 1530 integer(HID_T) :: file_id 1531 integer :: error, countf(2), offsetf(2) 1532 1533 PUSH_SUB(write_gvec_indices_hdf) 1534 1535 call open_file(99, trim(name), status='old') 1536 call close_file(99) 1537 1538 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 1539 countf(:) = (/ng, 1/) 1540 offsetf(:) = (/0, iq-1/) 1541 call hdf5_write_int_hyperslab(file_id, 'eps_header/gspace/gind_eps2rho', & 1542 countf, offsetf, gind_eps2rho, error) 1543 call hdf5_write_int_hyperslab(file_id, 'eps_header/gspace/gind_rho2eps', & 1544 countf, offsetf, gind_rho2eps, error) 1545 call hdf5_write_double_hyperslab(file_id, 'eps_header/gspace/ekin', & 1546 countf, offsetf, ekin, error) 1547 call h5fclose_f(file_id, error) 1548 1549 POP_SUB(write_gvec_indices_hdf) 1550 1551end subroutine write_gvec_indices_hdf 1552 1553!=================================================================================== 1554 1555subroutine write_vcoul_hdf(vcoul, iq, name) 1556 real(DP), intent(in) :: vcoul(:) !< (nmtx(iq)) 1557 integer, intent(in) :: iq 1558 character(len=*), intent(in) :: name 1559 1560 integer(HID_T) :: file_id 1561 integer :: error, countf(2), offsetf(2), nmtx(1) 1562 1563 PUSH_SUB(write_vcoul_hdf) 1564 1565 call open_file(99, trim(name), status='old') 1566 call close_file(99) 1567 1568 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 1569 call hdf5_read_int_hyperslab(file_id, 'eps_header/gspace/nmtx', (/1/), (/iq-1/), nmtx, error) 1570 if (size(vcoul)<nmtx(1)) & 1571 call die('Internal error: array vcoul too small in write_vcoul_hdf.', only_root_writes=.true.) 1572 countf(:) = (/nmtx(1), 1/) 1573 offsetf(:) = (/0, iq-1/) 1574 call hdf5_write_double_hyperslab(file_id, 'eps_header/gspace/vcoul', & 1575 countf, offsetf, vcoul, error) 1576 call h5fclose_f(file_id, error) 1577 1578 POP_SUB(write_vcoul_hdf) 1579 1580end subroutine write_vcoul_hdf 1581 1582 1583#endif 1584 1585end module write_matrix_m 1586 1587