1!>========================================================================= 2!! 3!! Module: 4!! 5!! epsread_hdf5_m Originally by JRD Last Modified 12/2014 (FHJ) 6!! 7!! Routines to read header info and matrices from epsmat files in 8!! HDF5 format. 9!! 10!!========================================================================= 11 12#include "f_defs.h" 13 14module epsread_hdf5_m 15#ifdef HDF5 16 17 use global_m 18 use hdf5 19 use hdf5_io_m 20 use scalapack_m 21 !use io_utils_m 22 use timing_m, only: timing => common_timing 23 implicit none 24 25 private 26 27 public :: & 28 read_eps_params_hdf5, & 29 read_eps_grid_sizes_hdf5, & 30 read_eps_subspace_info, & 31 read_eps_matrix_diagonal_hdf5, & 32 read_eps_qgrid_hdf5, & 33 read_eps_sub_neig_of_q, & 34 read_eps_freqgrid_hdf5, & 35 read_eps_old_gvecs_hdf5, & 36 read_eps_gvecsofq_hdf5, & 37 read_vcoul_hdf5, & 38 read_eps_matrix_col_f_hdf5, & 39 read_eps_matrix_col_hdf5, & 40 read_eps_matrix_ser_hdf5, & 41 read_eps_matrix_par_hdf5, & 42 read_eps_matrix_par_f_hdf5 43#if defined MPI && defined USESCALAPACK 44 public :: read_eps_matrix_par_hdf5_general, & 45 read_eps_matrix_par_hdf5_general_f 46#endif 47 48contains 49 50 51subroutine read_eps_params_hdf5(pol, name, nband) 52 type(polarizability), intent(inout) :: pol 53 character(len=*), intent(in) :: name 54 integer, intent(out), optional :: nband 55 56 integer(HID_T) :: file_id ! File identifier 57 integer :: error 58 59 PUSH_SUB(read_eps_params_hdf5) 60 61 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 62 call hdf5_read_int(file_id, 'eps_header/params/matrix_type', pol%matrix_type, error) 63 call hdf5_read_logical(file_id, 'eps_header/params/has_advanced', pol%has_advanced, error) 64 call hdf5_read_int(file_id, 'eps_header/params/nmatrix', pol%nmatrix, error) 65 call hdf5_read_int(file_id, 'eps_header/params/matrix_flavor', pol%matrix_flavor, error) 66 call hdf5_read_int(file_id, 'eps_header/params/icutv', pol%icutv, error) 67 call hdf5_read_double(file_id, 'eps_header/params/ecuts', pol%ecuts, error) 68 if (present(nband)) & 69 call hdf5_read_int(file_id, 'eps_header/params/nband', nband, error) 70 call hdf5_read_double(file_id, 'eps_header/params/efermi', pol%efermi, error) 71 pol%efermi = pol%efermi*ryd 72 call hdf5_read_int(file_id, 'eps_header/params/intraband_flag', pol%intraband_flag, error) 73 call hdf5_read_double(file_id, 'eps_header/params/intraband_overlap_min', pol%intraband_overlap_min, error) 74 call hdf5_read_logical(file_id, 'eps_header/params/subsample', pol%subsample, error) 75 call h5fclose_f(file_id, error) 76 77 POP_SUB(read_eps_params_hdf5) 78 79end subroutine read_eps_params_hdf5 80 81!=================================================================================== 82 83subroutine read_eps_grid_sizes_hdf5(ng, nq, ecuts, nfreq, nfreq_imag, nmtxmax, qgrid, freq_dep, name) 84 integer, intent(out) :: ng 85 integer, intent(out) :: nq 86 real(DP), intent(out) :: ecuts 87 integer, intent(out) :: nfreq 88 integer, intent(out) :: nfreq_imag 89 integer, intent(out) :: nmtxmax 90 integer, intent(out) :: qgrid(3) 91 integer, intent(out) :: freq_dep 92 character(len=*), intent(in) :: name 93 94 integer(HID_T) :: file_id ! File identifier 95 integer :: error, version, flavor 96 logical :: exists 97 98 PUSH_SUB(read_eps_grid_sizes_hdf5) 99 100 call open_file(99, trim(name), status='old') 101 call close_file(99) 102 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 103 call hdf5_require_version(file_id, 'eps_header/versionnumber', VER_EPS_HDF5, trim(name)) 104 call hdf5_require_version(file_id, 'mf_header/versionnumber', VER_WFN_HDF5, trim(name)) 105 call hdf5_require_flavor(file_id, 'eps_header/flavor', SCALARSIZE, trim(name)) 106 call hdf5_require_flavor(file_id, 'mf_header/flavor', SCALARSIZE, trim(name)) 107 108 call hdf5_read_int(file_id, 'eps_header/qpoints/nq', nq, error) 109 call hdf5_read_int(file_id, 'mf_header/gspace/ng', ng, error) 110 call hdf5_read_int(file_id, 'eps_header/freqs/nfreq', nfreq, error) 111 call h5lexists_f(file_id, 'eps_header/freqs/nfreq_imag', exists, error) 112 if (exists.and.error==0) then 113 call hdf5_read_int(file_id, 'eps_header/freqs/nfreq_imag', nfreq_imag, error) 114 else 115 nfreq_imag = 0 116 endif 117 call hdf5_read_int(file_id, 'eps_header/freqs/freq_dep', freq_dep, error) 118 call hdf5_read_int(file_id, 'eps_header/gspace/nmtx_max', nmtxmax, error) 119 call hdf5_read_int_array(file_id, 'eps_header/qpoints/qgrid', (/3/), qgrid, error) 120 call hdf5_read_double(file_id, 'eps_header/params/ecuts', ecuts, error) 121 call h5fclose_f(file_id,error) 122 123 POP_SUB(read_eps_grid_sizes_hdf5) 124 125end subroutine read_eps_grid_sizes_hdf5 126 127!==================================================================================== 128 129subroutine read_eps_subspace_info(is_subspace, matrix_in_subspace_basis, keep_full_eps_static, & 130 neig_sub, chi_eigenvalue_cutoff, name) 131 logical, intent(out) :: is_subspace, matrix_in_subspace_basis, keep_full_eps_static 132 integer, intent(out) :: neig_sub 133 real(dp) :: chi_eigenvalue_cutoff 134 character(len=*), intent(in) :: name 135 136 integer(HID_T) :: file_id ! File identifier 137 integer :: error 138 logical :: exists 139 140 PUSH_SUB(read_eps_subspace_info) 141 142 is_subspace = .false. 143 matrix_in_subspace_basis = .false. 144 keep_full_eps_static = .false. 145 neig_sub = 0 146 chi_eigenvalue_cutoff = 0.0d0 147 148 call open_file(99, trim(name), status='old') 149 call close_file(99) 150 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 151 152 ! check if the subspace flag exists 153 call h5lexists_f(file_id, 'eps_header/params/subspace', exists, error) 154 if (exists .and. error==0) then 155 call hdf5_read_logical(file_id, 'eps_header/params/subspace', is_subspace, error) 156 if( is_subspace ) then 157 call hdf5_read_logical(file_id, 'eps_header/subspace/matrix_in_subspace_basis', matrix_in_subspace_basis, error) 158 call hdf5_read_logical(file_id, 'eps_header/subspace/keep_full_eps_static', keep_full_eps_static, error) 159 call hdf5_read_double(file_id, 'eps_header/subspace/eps_eigenvalue_cutoff', chi_eigenvalue_cutoff, error) 160 call hdf5_read_int(file_id, 'eps_header/subspace/neig_max', neig_sub, error) 161 end if 162 end if 163 call h5fclose_f(file_id,error) 164 165 POP_SUB(read_eps_subspace_info) 166 167end subroutine 168 169!==================================================================================== 170 171subroutine read_eps_matrix_diagonal_hdf5(nmtx, iq, epsdiag, name) 172 integer, intent(in) :: nmtx 173 integer, intent(in) :: iq 174 real(DP), intent(out) :: epsdiag(:,:) !< (matrix_flavor,nmtx) 175 character(len=*), intent(in) :: name 176 177 integer(HID_T) :: file_id ! File identifier 178 integer :: error, matrix_flavor 179 180 PUSH_SUB(read_eps_matrix_diagonal_hdf5) 181 182 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 183 call hdf5_read_int(file_id, 'eps_header/params/matrix_flavor', matrix_flavor, error) 184 if (matrix_flavor/=size(epsdiag,1)) then 185 write(0,*) 'ERROR: Got size(epsdiag,1)=',size(epsdiag,1),', but matrix_flavor=',matrix_flavor 186 call die('Internal error in read_eps_matrix_diagonal_hdf5', & 187 only_root_writes=.true.) 188 endif 189 call hdf5_read_double_hyperslab(file_id, 'mats/matrix-diagonal', & 190 (/matrix_flavor,nmtx,1/), (/0,0,iq-1/), epsdiag, error) 191 call h5fclose_f(file_id,error) 192 193 POP_SUB(read_eps_matrix_diagonal_hdf5) 194 195end subroutine read_eps_matrix_diagonal_hdf5 196 197!==================================================================================== 198 199subroutine read_eps_qgrid_hdf5(nq, qpts, nmtx, name) 200 integer, intent(in) :: nq 201 real(DP), intent(inout) :: qpts(:,:) !< (3,nq) 202 integer, intent(out) :: nmtx(:) !< (nq) 203 character(len=*), intent(in) :: name 204 205 integer(HID_T) :: file_id ! File identifier 206 integer :: error 207 208 PUSH_SUB(read_eps_qgrid_hdf5) 209 210 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 211 call hdf5_read_double_array(file_id, 'eps_header/qpoints/qpts', (/3,nq/), qpts, error) 212 call hdf5_read_int_array(file_id, 'eps_header/gspace/nmtx', (/nq/), nmtx, error) 213 call h5fclose_f(file_id, error) 214 215 POP_SUB(read_eps_qgrid_hdf5) 216 217end subroutine read_eps_qgrid_hdf5 218 219!=================================================================================== 220 221subroutine read_eps_sub_neig_of_q(nq, neig, name) 222 223 integer, intent(in) :: nq 224 integer, intent(out) :: neig(:) !< (nq) 225 character(len=*), intent(in) :: name 226 227 integer(HID_T) :: file_id ! File identifier 228 integer :: error 229 230 PUSH_SUB(read_eps_sub_neig_of_q) 231 232 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 233 call hdf5_read_int_array(file_id, 'eps_header/subspace/neig', (/nq/), neig, error) 234 call h5fclose_f(file_id, error) 235 236 POP_SUB(read_eps_sub_neig_of_q) 237 238end subroutine read_eps_sub_neig_of_q 239 240!=================================================================================== 241 242subroutine read_eps_freqgrid_hdf5(nfreq, dFreqGrid, dFreqBrd, name) 243 integer, intent(in) :: nfreq 244 real(DP), intent(out) :: dFreqGrid(:) !< (nfreq) 245 complex(DPC), intent(out) :: dFreqBrd(:) !< (nfreq) 246 character(len=*), intent(in) :: name 247 248 real(DP) :: freqs_tmp(2,nfreq) 249 integer :: iw 250 integer(HID_T) :: file_id ! File identifier 251 integer :: error 252 253 PUSH_SUB(read_eps_freqgrid_hdf5) 254 255 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 256 call hdf5_read_double_array(file_id, 'eps_header/freqs/freqs', (/2,nfreq/), freqs_tmp, error) 257 do iw=1,nfreq 258 dFreqGrid(iw) = freqs_tmp(1,iw) 259 dFreqBrd(iw) = CMPLX(0,freqs_tmp(2,iw)) 260 enddo 261 call h5fclose_f(file_id,error) 262 263 POP_SUB(read_eps_freqgrid_hdf5) 264 265end subroutine read_eps_freqgrid_hdf5 266 267!================================================================================= 268 269subroutine read_eps_old_gvecs_hdf5(ng, gvecs, name) 270 integer, intent(in) :: ng 271 integer, intent(out) :: gvecs(:,:) !< (3,ng) 272 character(len=*), intent(in) :: name 273 274 integer(HID_T) :: file_id ! File identifier 275 integer :: error 276 277 PUSH_SUB(read_eps_old_gvecs_hdf5) 278 279 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 280 call hdf5_read_int_array(file_id, 'mf_header/gspace/components', (/3,ng/), gvecs, error) 281 call h5fclose_f(file_id,error) 282 283 POP_SUB(read_eps_old_gvecs_hdf5) 284 285end subroutine read_eps_old_gvecs_hdf5 286 287!==================================================================================== 288 289subroutine read_eps_gvecsofq_hdf5(ng, gind_eps2rho, gind_rho2eps, ekin, iq, name) 290 integer, intent(in) :: ng !< Number of G-vectors 291 integer, intent(out) :: gind_eps2rho(:) !< (ng) 292 integer, intent(out) :: gind_rho2eps(:) !< (ng) 293 real(DP), intent(out) :: ekin(:) !< (ng) 294 integer, intent(in) :: iq 295 character(len=*), intent(in) :: name 296 297 integer(HID_T) :: file_id 298 integer :: error 299 integer :: countf(2), offsetf(2) 300 301 PUSH_SUB(read_eps_gvecsofq_hdf5) 302 303 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 304 305 countf(:) = (/ng, 1/) 306 offsetf(:) = (/0, iq-1/) 307 call hdf5_read_int_hyperslab(file_id, 'eps_header/gspace/gind_eps2rho', & 308 countf, offsetf, gind_eps2rho, error) 309 call hdf5_read_int_hyperslab(file_id, 'eps_header/gspace/gind_rho2eps', & 310 countf, offsetf, gind_rho2eps, error) 311 call hdf5_read_double_hyperslab(file_id, 'eps_header/gspace/ekin', & 312 countf, offsetf, ekin, error) 313 314 call h5fclose_f(file_id, error) 315 316 POP_SUB(read_eps_gvecsofq_hdf5) 317 318end subroutine read_eps_gvecsofq_hdf5 319 320!==================================================================================== 321 322subroutine read_vcoul_hdf5(vcoul, iq, name) 323 real(DP), intent(out) :: vcoul(:) !< (nmtx(iq)) 324 integer, intent(in) :: iq 325 character(len=*), intent(in) :: name 326 327 integer(HID_T) :: file_id 328 integer :: error, nmtx(1) 329 integer :: countf(2), offsetf(2) 330 331 PUSH_SUB(read_vcoul_hdf5) 332 333 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 334 335 call hdf5_read_int_hyperslab(file_id, 'eps_header/gspace/nmtx', (/1/), (/iq-1/), nmtx, error) 336 if (size(vcoul)<nmtx(1)) & 337 call die('Internal error: array vcoul too small in read_vcoul_hdf.', only_root_writes=.true.) 338 countf(:) = (/nmtx(1), 1/) 339 offsetf(:) = (/0, iq-1/) 340 call hdf5_read_double_hyperslab(file_id, 'eps_header/gspace/vcoul', & 341 countf, offsetf, vcoul, error) 342 343 call h5fclose_f(file_id, error) 344 345 POP_SUB(read_vcoul_hdf5) 346 347end subroutine read_vcoul_hdf5 348 349!=========================================================================================== 350 351subroutine read_eps_matrix_col_f_hdf5(retarded, nFreq, igp, nmtx, iq, is, name, advanced) 352 integer, intent(in) :: iq 353 integer, intent(in) :: is 354 integer, intent(in) :: nFreq 355 integer, intent(in) :: nmtx 356 integer, intent(in) :: igp 357 complex(DPC), intent(out) :: retarded(:,:) !< (nmtx,nFreq) 358 character(len=*), intent(in) :: name 359 complex(DPC), optional, intent(out) :: advanced(:,:) !< (nmtx,nFreq) 360 361 integer(HID_T) :: file_id ! File identifier 362 integer(HID_T) :: data_id ! Property list identifier 363 integer(HID_T) :: dataspace ! Property list identifier 364 integer(HID_T) :: memspace ! Property list identifier 365 integer :: error 366 integer(HSIZE_T) :: count(6), offset(6) 367 368 real(DP), allocatable :: data(:,:,:,:,:,:) 369 integer :: nmatrix_per_spin, nspin, buf_sz, version 370 logical :: has_advanced 371 372 PUSH_SUB(read_eps_matrix_col_f_hdf5) 373 374 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 375 376 ! FHJ: the default is never to read the advanced matrix, unless the file 377 ! version is <3 (on which case we didn`t store the Coulomb interaction) 378 call hdf5_read_int(file_id, 'eps_header/versionnumber', version, error) 379 call hdf5_read_int(file_id, 'eps_header/params/nmatrix', nmatrix_per_spin, error) 380 call hdf5_read_int(file_id, 'mf_header/kpoints/nspin', nspin, error) 381 nmatrix_per_spin = nmatrix_per_spin / nspin 382 has_advanced = .false. 383 buf_sz = 1 384 if (version<3) then 385 call hdf5_read_logical(file_id, 'eps_header/params/has_advanced', has_advanced, error) 386 if (present(advanced) .and. .not.has_advanced) & 387 call die('Inconsistent epsmat file: version<3, but no advanced matrix', only_root_writes=.true.) 388 endif 389 if (has_advanced) buf_sz = 2 390 391 call h5dopen_f(file_id, 'mats/matrix', data_id, error) 392 call h5dget_space_f(data_id,dataspace,error) 393 394 count(:) = (/2, nmtx, 1, nFreq, buf_sz, 1/) 395 call h5screate_simple_f(6, count, memspace, error) 396 offset(:) = (/0, 0, igp-1, 0, nmatrix_per_spin*(is-1), iq-1/) 397 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 398 399! XXX - this is bad because we double memory usage. Oh well. 400 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 401 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace) 402 retarded(1:nmtx, 1:nFreq) = CMPLX(data(1,1:nmtx,1,1:nFreq,1,1),data(2,1:nmtx,1,1:nFreq,1,1)) 403 if (has_advanced .and. present(advanced)) then 404 advanced(1:nmtx,1:nFreq) = CMPLX(data(1,1:nmtx,1,1:nFreq,2,1),data(2,1:nmtx,1,1:nFreq,2,1)) 405 endif 406 SAFE_DEALLOCATE(data) 407 call h5sclose_f(memspace, error) 408 409 if (.not.has_advanced .and. present(advanced)) call get_advanced_from_retarded() 410 411 call h5sclose_f(dataspace, error) 412 call h5dclose_f(data_id,error) 413 call h5fclose_f(file_id,error) 414 415 POP_SUB(read_eps_matrix_col_f_hdf5) 416 417contains 418 419 !> FHJ: This routine does the magic of generating advanced epsilon/chi out of 420 !! the retarded epsilon/chi. 421 !! We use the fact that W_A = (W_R)^H to write 422 !! epsinv_A = (epsinv_R * v)^H * (1/v) 423 subroutine get_advanced_from_retarded() 424 real(DP) :: vcoul(nmtx) 425 integer :: matrix_type, ig 426 427 PUSH_SUB(read_eps_matrix_col_f_hdf5.get_advanced_from_retarded) 428 429 call hdf5_read_int(file_id, 'eps_header/params/matrix_type', matrix_type, error) 430 vcoul(:) = 1d0 431 ! FHJ: There is no Coulomb interaction in a chimat file. 432 if (matrix_type<2) then 433 call hdf5_read_double_hyperslab(file_id, 'eps_header/gspace/vcoul', & 434 (/nmtx,1/), (/0,iq-1/), vcoul, error) 435 endif 436 437 count(:) = (/2, 1, nmtx, nFreq, buf_sz, 1/) 438 call h5screate_simple_f(6, count, memspace, error) 439 offset(:) = (/0, igp-1, 0, 0, nmatrix_per_spin*(is-1), iq-1/) 440 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 441 442 ! XXX - this is bad because we double memory usage. Oh well. 443 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 444 445 ! JRD XXX These reads have horrendous locality 446 ! FHJ: you are right, they have terrible I/O locality. But this is only 447 ! for the serial code. 448 449 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace) 450 ! FHJ: Note the complex conjugation here 451 advanced(1:nmtx,1:nFreq) = CMPLX(data(1,1,1:nmtx,1:nFreq,1,1),-data(2,1,1:nmtx,1:nFreq,1,1)) 452 ! FHJ: Correct asymmetry from Coulomb interaction 453 do ig = 1, nmtx 454 advanced(ig,1:nFreq) = advanced(ig,1:nFreq) * (vcoul(ig)/vcoul(igp)) 455 enddo 456 SAFE_DEALLOCATE(data) 457 call h5sclose_f(memspace, error) 458 459 POP_SUB(read_eps_matrix_col_f_hdf5.get_advanced_from_retarded) 460 461 end subroutine get_advanced_from_retarded 462 463 464end subroutine read_eps_matrix_col_f_hdf5 465 466!=========================================================================================== 467 468subroutine read_eps_matrix_col_hdf5(eps,j,nmtx,iq,is,name) 469 integer, intent(in) :: iq 470 integer, intent(in) :: is 471 integer, intent(in) :: nmtx 472 integer, intent(in) :: j 473 SCALAR, intent(out) :: eps(:) !< (nmtx) 474 character(len=*), intent(in) :: name 475 476 integer(HID_T) :: file_id ! File identifier 477 integer(HID_T) :: data_id ! Property list identifier 478 integer(HID_T) :: dataspace ! Property list identifier 479 integer(HID_T) :: memspace ! Property list identifier 480 integer :: error, rank 481 integer(HSIZE_T) :: count(6), offset(6) 482 483 real(DP), allocatable :: data(:,:,:,:,:,:) 484 485 PUSH_SUB(read_eps_matrix_col_hdf5) 486 487 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 488 489 call h5dopen_f(file_id, 'mats/matrix', data_id, error) 490 call h5dget_space_f(data_id,dataspace,error) 491 492 rank = 6 493 count(1) = SCALARSIZE 494 count(2) = nmtx 495 count(3) = 1 496 count(4) = 1 497 count(5) = 1 !mat 498 count(6) = 1 !iq 499 500 call h5screate_simple_f(rank, count, memspace, error) 501 502 offset(:) = 0 503 offset(3) = j - 1 504 offset(5) = is - 1 505 offset(6) = iq - 1 506 507 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 508 509! XXX - this is bad because we double memory usage. Oh well. 510 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 511 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace) 512! write(6,*) 'About to die?', nmtx, data(1,1,:,1,1,1) 513 eps(1:nmtx) = SCALARIFY2(data(1,1:nmtx,1,1,1,1),data(2,1:nmtx,1,1,1,1)) 514 SAFE_DEALLOCATE(data) 515 516 call h5sclose_f(memspace, error) 517 call h5sclose_f(dataspace, error) 518 call h5dclose_f(data_id,error) 519 call h5fclose_f(file_id,error) 520 521 POP_SUB(read_eps_matrix_col_hdf5) 522 523end subroutine read_eps_matrix_col_hdf5 524 525!=========================================================================================== 526 527subroutine read_eps_matrix_ser_hdf5(eps,nmtx,iq,is,name) 528 integer, intent(in) :: iq 529 integer, intent(in) :: is 530 integer, intent(in) :: nmtx 531 SCALAR, intent(out) :: eps(:,:) !< (nmtx,nmtx) 532 character(len=*), intent(in) :: name 533 534 integer(HID_T) :: file_id ! File identifier 535 integer(HID_T) :: data_id ! Property list identifier 536 integer(HID_T) :: dataspace ! Property list identifier 537 integer(HID_T) :: memspace ! Property list identifier 538 integer :: error, rank 539 integer(HSIZE_T) :: count(6), offset(6) 540 541 real(DP), allocatable :: data(:,:,:,:,:,:) 542 543 PUSH_SUB(read_eps_matrix_ser_hdf5) 544 545 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 546 547 call h5dopen_f(file_id, 'mats/matrix', data_id, error) 548 call h5dget_space_f(data_id,dataspace,error) 549 550 rank = 6 551 count(1) = SCALARSIZE 552 count(2) = nmtx 553 count(3) = nmtx 554 count(4) = 1 555 count(5) = 1 !mat 556 count(6) = 1 !iq 557 558 call h5screate_simple_f(rank, count, memspace, error) 559 560 offset(:) = 0 561 offset(5) = is - 1 562 offset(6) = iq - 1 563 564 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 565 566! XXX - this is bad because we double memory usage. Oh well. 567 SAFE_ALLOCATE(data,(count(1),count(2),count(3),count(4),count(5),count(6))) 568 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace) 569 eps(1:nmtx,1:nmtx) = SCALARIFY2(data(1,1:nmtx,1:nmtx,1,1,1),data(2,1:nmtx,1:nmtx,1,1,1)) 570 SAFE_DEALLOCATE(data) 571 572 call h5sclose_f(memspace, error) 573 call h5sclose_f(dataspace, error) 574 call h5dclose_f(data_id,error) 575 call h5fclose_f(file_id,error) 576 577 POP_SUB(read_eps_matrix_ser_hdf5) 578 579end subroutine read_eps_matrix_ser_hdf5 580 581!=========================================================================================== 582 583subroutine read_eps_matrix_par_hdf5(eps, nb, rank, npes, nmtx, iq, is, name) 584 SCALAR, intent(inout) :: eps(:,:) !< (neps,ngpown) 585 integer, intent(in) :: nb !< block size 586 integer, intent(in) :: rank !< processor rank for column distribution 587 integer, intent(in) :: npes !< number of processors over which we distribute 588 integer, intent(in) :: nmtx 589 integer, intent(in) :: iq 590 integer, intent(in) :: is 591 character(len=*), intent(in) :: name 592 593 integer(HID_T) :: file_id ! File identifier 594 integer(HID_T) :: data_id ! Property list identifier 595 integer(HID_T) :: plist_id ! Property list identifier 596 integer(HID_T) :: dataspace ! Property list identifier 597 integer(HID_T) :: memspace ! Property list identifier 598 integer :: error 599 integer :: comm, info 600 integer(HSIZE_T) :: count(6), offset(6), countm(6) 601 real(DP), allocatable :: data(:,:,:,:,:,:) 602 integer :: ngpown_max, igp, igp_loc 603 604 PUSH_SUB(read_eps_matrix_par_hdf5) 605 606 ngpown_max = NUMROC(nmtx, nb, 0, 0, npes) 607 SAFE_ALLOCATE(data,(SCALARSIZE,nmtx,1,1,1,1)) 608 609#ifdef MPI 610 comm = MPI_COMM_WORLD 611 info = MPI_INFO_NULL 612 613 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 614 call h5pset_fapl_mpio_f(plist_id, comm, info, error) 615 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error, access_prp = plist_id) 616 call h5pclose_f(plist_id,error) 617#else 618 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 619#endif 620 621 622 do igp_loc = 1, ngpown_max 623 624 igp = INDXL2G(igp_loc, nb, rank, 0, npes) 625 call h5dopen_f(file_id, 'mats/matrix', data_id, error) 626 call h5dget_space_f(data_id,dataspace,error) 627 628! JRD: The commented code in this routine represents efforts to use a single HDF5 read call 629! with an appropriate block and stride for each proc. In general, it currently appears to 630! perform worse (though this may be related to size of matrix. So, it is commented until 631! further investigation. 632 633 countm(1) = SCALARSIZE 634 countm(2) = nmtx 635 countm(3) = 1 636 countm(4) = 1 637 countm(5) = 1 638 countm(6) = 1 639 640 call h5screate_simple_f(6, countm, memspace, error) 641 642 count(1) = SCALARSIZE 643 count(2) = nmtx 644 count(3) = 1 645 count(4) = 1 646 count(5) = 1 647 count(6) = 1 648 649! Construct data and offset 650 651 if (igp <= nmtx) then 652 !if (ngpown .ne. 0 .and. my_igp .le. nmtx) then 653 654 offset(1)=0 655 offset(2)=0 656 offset(3)=igp-1 657 offset(4)=0 658 offset(5)=is-1 659 offset(6)=iq-1 660 661! Select hyperslab 662 663 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 664 665 else 666 667 call H5sselect_none_f(memspace,error); 668 call H5sselect_none_f(dataspace,error); 669 670 endif 671 672! Create property list for collective dataset read 673 674#ifdef MPI 675 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 676 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 677 678! Collectively read the file 679 680 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, dataspace, & 681 xfer_prp = plist_id) 682# else 683 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, dataspace) 684#endif 685 686 if (igp <= nmtx) then 687! if (ngpown .ne. 0 .and. my_igp .le. nmtx) then 688!XXX PROBABLY NEED THREADED LOOP HERE 689 eps(1:nmtx,igp_loc) = SCALARIFY2(data(1,1:nmtx,1,1,1,1),data(2,1:nmtx,1,1,1,1)) 690 endif 691 692 call h5sclose_f(memspace, error) 693 call h5sclose_f(dataspace, error) 694 call h5dclose_f(data_id, error) 695 696 enddo 697 698 call h5fclose_f(file_id,error) 699 SAFE_DEALLOCATE(data) 700 701 POP_SUB(read_eps_matrix_par_hdf5) 702 703end subroutine read_eps_matrix_par_hdf5 704 705!========================================================================================== 706 707subroutine read_eps_matrix_par_f_hdf5(retarded, nb, pool_comm, my_pool, npools, & 708 nmtx, nFreq, iq, is, name, advanced) 709 complex(DPC), intent(inout) :: retarded(:,:,:) !< (neps,ngpown,nFreq) 710 integer, intent(in) :: nb !< block size 711 integer, intent(in) :: pool_comm !< MPI comm for each pool 712 integer, intent(in) :: my_pool !< my pool, starting from 0 (0=no pools) 713 integer, intent(in) :: npools !< number of pools (1=no pools). 714 integer, intent(in) :: nmtx 715 integer, intent(in) :: nFreq 716 integer, intent(in) :: iq 717 integer, intent(in) :: is 718 character(len=*), intent(in) :: name 719 complex(DPC), optional, intent(inout) :: advanced(:,:,:) !< (nFreq,neps,ngpown) 720 721 integer(HID_T) :: file_id ! File identifier 722 integer(HID_T) :: data_id ! Property list identifier 723 integer(HID_T) :: plist_id ! Property list identifier 724 integer(HID_T) :: dataspace ! Property list identifier 725 integer(HID_T) :: memspace ! Property list identifier 726 integer :: error 727 integer(HSIZE_T) :: count(6), offset(6) 728 real(DP), allocatable :: data(:,:,:,:,:,:) 729 730 integer :: pool_rank !< processor rank for column distribution 731 integer :: npes_pool !< number of processors over which we distribute 732 integer :: ngpown_max, igp, igp_loc, nmatrix_per_spin, nspin, buf_sz, version 733 logical :: want_advanced, read_advanced 734 !type(progress_info) :: prog_info !< a user-friendly progress report 735 736 integer :: iw 737 738 PUSH_SUB(read_eps_matrix_par_f_hdf5) 739 740#ifdef MPI 741 call MPI_Comm_rank(pool_comm, pool_rank, mpierr) 742 call MPI_Comm_size(pool_comm, npes_pool, mpierr) 743 ! FHJ: We need the following BCast for processors left over out of the pool 744 call MPI_Bcast(npes_pool, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 745 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 746 call h5pset_fapl_mpio_f(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL, error) 747 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error, access_prp = plist_id) 748 call h5pclose_f(plist_id,error) 749#else 750 pool_rank = 0 751 npes_pool = 1 752 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 753#endif 754 755 want_advanced = present(advanced) 756 ! FHJ: the default is never to read the advanced matrix, unless the file 757 ! version is <3 (on which case we didn`t store the Coulomb interaction) 758 call hdf5_read_int(file_id, 'eps_header/versionnumber', version, error) 759 call hdf5_read_int(file_id, 'eps_header/params/nmatrix', nmatrix_per_spin, error) 760 call hdf5_read_int(file_id, 'mf_header/kpoints/nspin', nspin, error) 761 nmatrix_per_spin = nmatrix_per_spin / nspin 762 read_advanced = .false. 763 buf_sz = 1 764 if (version<3) then 765 call hdf5_read_logical(file_id, 'eps_header/params/has_advanced', read_advanced, error) 766 if (SCALARSIZE==2 .and. .not.read_advanced) & 767 call die('Inconsistent epsmat file: version<3, but no advanced matrix', only_root_writes=.true.) 768 endif 769 if (read_advanced) buf_sz = 2 770 771 ngpown_max = NUMROC(nmtx, nb, 0, 0, npes_pool) 772 SAFE_ALLOCATE(data,(2,nmtx,1,nFreq,buf_sz,1)) 773 774 call logit('Reading HDF5 file') 775 !call progress_init(prog_info, 'reading '//trim(name), & 776 ! 'distributed column', ngpown_max) 777! JRD XXX, we should do this as a single read instead of loop!!! 778 do igp_loc = 1, ngpown_max 779 !call progress_step(prog_info) 780 if (my_pool>=0) then 781 igp = INDXL2G(igp_loc, nb, pool_rank, 0, npes_pool) 782 else 783 igp = nmtx + 1 784 endif 785 call h5dopen_f(file_id, 'mats/matrix', data_id, error) 786 call h5dget_space_f(data_id,dataspace,error) 787 788 count(:) = (/2, nmtx, 1, nFreq, buf_sz, 1/) 789 call h5screate_simple_f(6, count, memspace, error) 790 791 ! Construct data and offset 792 if (igp <= nmtx) then 793 offset = (/0, 0, igp-1, 0, nmatrix_per_spin*(is-1), iq-1/) 794 call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, offset, count, error) 795 else 796 call H5sselect_none_f(memspace,error); 797 call H5sselect_none_f(dataspace,error); 798 endif 799 800! Create property list for collective dataset read 801! Read is serial for now 802 803#ifdef MPI 804 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 805 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 806 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace, & 807 xfer_prp = plist_id) 808#else 809 call h5dread_f(data_id, H5T_NATIVE_DOUBLE, data, count, error, memspace, dataspace) 810#endif 811 812 if (igp <= nmtx) then 813 retarded(1:nmtx,igp_loc,1:nFreq) = CMPLX(data(1,1:nmtx,1,1:nFreq,1,1),data(2,1:nmtx,1,1:nFreq,1,1)) 814 if (want_advanced .and. read_advanced) then 815 advanced(1:nmtx,igp_loc,1:nFreq) = CMPLX(data(1,1:nmtx,1,1:nFreq,2,1),data(2,1:nmtx,1,1:nFreq,2,1)) 816 endif 817 endif 818 819 call h5sclose_f(memspace, error) 820 call h5sclose_f(dataspace, error) 821 call h5dclose_f(data_id, error) 822 enddo 823 824 !call progress_free(prog_info) 825 SAFE_DEALLOCATE(data) 826 827 if (want_advanced .and. .not.read_advanced) call get_advanced_from_retarded() 828 829 call h5fclose_f(file_id,error) 830 831 POP_SUB(read_eps_matrix_par_f_hdf5) 832 833contains 834 835 !> FHJ: This routine does the magic of generating advanced epsilon/chi out of 836 !! the retarded epsilon/chi. 837 !! We use the fact that W_A = (W_R)^H to write 838 !! epsinv_A = (epsinv_R * v)^H * (1/v) 839 subroutine get_advanced_from_retarded() 840 real(DP) :: vcoul(nmtx) 841 integer :: ig, igp, igp_loc, matrix_type 842#ifdef MPI 843 complex(DPC) :: send_buf(ngpown_max,ngpown_max,nFreq,npes_pool), recv_buf(ngpown_max,ngpown_max,nFreq,npes_pool) 844 integer :: req_send(npes_pool), req_recv(npes_pool) 845 integer :: ig_loc, ngown_ipe, ngpown, ipe 846#endif 847 848 PUSH_SUB(read_eps_matrix_par_f_hdf5.get_advanced_from_retarded) 849 850 call logit('Getting advanced matrix from retarded') 851 call hdf5_read_int(file_id, 'eps_header/params/matrix_type', matrix_type, error) 852 vcoul(:) = 1d0 853 ! FHJ: There is no Coulomb interaction in a chimat file. 854 if (matrix_type<2) then 855 call hdf5_read_double_hyperslab(file_id, 'eps_header/gspace/vcoul', & 856 (/nmtx,1/), (/0,iq-1/), vcoul, error) 857 endif 858 859 if (my_pool<0) then 860 POP_SUB(read_eps_matrix_par_f_hdf5.get_advanced_from_retarded) 861 return 862 endif 863 864#ifdef MPI 865 ngpown = NUMROC(nmtx, nb, pool_rank, 0, npes_pool) 866 867 call logit('Preparing send buffers') 868 ! FHJ: Prepare the send buffer. The send buffer will retarded^H, but with the 869 ! columns of retarded^H (rows of retarded) reordered from ig=1..nmtx to irow=(ig_loc,ipe). 870 ! We do this to distribute all the matrix elements with a single MPI_Scatter. 871 send_buf = (0d0, 0d0) 872 do ipe = 0, npes_pool-1 873 ngown_ipe = NUMROC(nmtx, nb, ipe, 0, npes_pool) 874 ! JRD XXX Loops not ordered well. FHJ this well be killer on systems with lots of G vecs 875 ! FHJ: We are doing two matrix transpositions. We will *never* get good locality. 876 ! The indices are optimized here for MPI communication. I actually timed this 877 ! and saw that, with the opposite index order, MPI communication kills. 878 ! For large ng, we are still paying the price of a large matrix transposition, 879 ! which is still cheaper than lots of unordered MPI_BCasts`s or MPI_Gather`s. 880 do ig_loc = 1, ngown_ipe 881 ig = INDXL2G(ig_loc, nb, ipe, 0, npes_pool) 882 send_buf(1:ngpown,ig_loc,1:nFreq,ipe+1) = CONJG(retarded(ig,1:ngpown,1:nFreq)) 883 enddo 884 enddo 885 886 call logit('Setting up non-blocking communication') 887 call timing%start(timing%epscopy_comm) 888 ! FHJ: Each processor scatters its retarded^H matrix, stored in send_buf, to all 889 ! all other processors. Each PE gets nFreq*ngpown_max**2 complex numbers, 890 ! stored as a (ngpown_max,ngpown_max) matrix of vectors of length nFreq. 891 ! FHJ: Open non-blocking receiving buffers 892 do ipe = 0, npes_pool-1 893 call MPI_Irecv(recv_buf(1,1,1,ipe+1), nFreq*ngpown_max**2, MPI_COMPLEX_DPC, ipe, & 894 pool_rank*npes_pool + ipe, pool_comm, req_recv(ipe+1), mpierr) 895 enddo 896 ! FHJ: Send data in non-blocking way 897 do ipe = 0, npes_pool-1 898 call MPI_Isend(send_buf(1,1,1,ipe+1), nFreq*ngpown_max**2, MPI_COMPLEX_DPC, ipe, & 899 ipe*npes_pool + pool_rank, pool_comm, req_send(ipe+1), mpierr) 900 enddo 901 902 call logit('Receiving data') 903 do while (.true.) 904 ! FHJ: Work on whichever buffer arrived first. 905 call MPI_Waitany(npes_pool, req_recv, ipe, MPI_STATUS_IGNORE, mpierr) 906 if (ipe==MPI_UNDEFINED) exit 907 ipe = ipe - 1 908 ! FHJ: Processor "ipe" was kind enough to reorder its rows of retarded such in 909 ! a way that it corresponds to all columns of advanced that we own, and in the 910 ! correct order. However, each "ipe" originally owned a subset of all 911 ! columns of retarded, so here we received a subset of the rows of advanced. 912 ngown_ipe = NUMROC(nmtx, nb, ipe, 0, npes_pool) 913 ! JRD: XXX this loop is horrifying, FHJ - what were you thinking?? 914 ! FHJ: Again, we are doing a **two matrix transpositions in parallel**. 915 ! We will *never* get good locality. 916 ! Seriously, before complaining so much about this line, did you see 917 ! how efficient the MPI communication is -- which has a much smaller bandwitdh 918 ! than this in-memory transposition?! And how we send a small number of 919 ! MPI messages?! And how much all this is *much faster* than reading 920 ! data from disk twice, which has an even smaller bandwidth?! 921 ! There is no free lunch, you either pay the price of MPI communication 922 ! + disk I/O or memory transposition. 923 ! Memory transposition is still much cheaper, and I timed it. 924 ! And BTW: this code is faster than BLACS for the matrix transposition. 925 do iw = 1, nFreq 926 do igp_loc = 1, ngpown 927 do ig_loc = 1, ngown_ipe 928 ig = INDXL2G(ig_loc, nb, ipe, 0, npes_pool) 929 advanced(ig,igp_loc,iw) = recv_buf(ig_loc,igp_loc,iw,ipe+1) 930 enddo 931 enddo 932 enddo 933 enddo 934 935 call logit('Releasing send buffer') 936 call MPI_Waitall(npes_pool, req_send, MPI_STATUSES_IGNORE, mpierr) 937 call timing%stop(timing%epscopy_comm) 938 939 call logit('Fixing asymmetry in the Coulomb interaction') 940 ! FHJ: Finally, fix the asymmetry in the Coulomb interaction. 941 ! JRD XXX Should OMP here 942 ! FHJ: ok, but note that nFreq might be small. 943 do iw = 1, nFreq 944 do igp_loc = 1, ngpown 945 igp = INDXL2G(igp_loc, nb, pool_rank, 0, npes_pool) 946 do ig = 1, nmtx 947 advanced(ig,igp_loc,iw) = advanced(ig,igp_loc,iw) * (vcoul(ig)/vcoul(igp)) 948 enddo 949 enddo 950 enddo 951#else 952 ! FHJ: serial is just slightly simpler 953 do igp = 1, nmtx 954 do ig = 1, nmtx 955 advanced(ig,igp,1:nFreq) = CONJG(retarded(igp,ig,1:nFreq)) * (vcoul(ig)/vcoul(igp)) 956 enddo 957 enddo 958#endif 959 960 POP_SUB(read_eps_matrix_par_f_hdf5.get_advanced_from_retarded) 961 962 end subroutine get_advanced_from_retarded 963 964end subroutine read_eps_matrix_par_f_hdf5 965 966!========================================================================================== 967! Routine used in SUBSPACE method, can be used for other application also 968#if defined MPI && defined USESCALAPACK 969 970subroutine read_eps_matrix_par_hdf5_general(scal, matrix, nmtx, nmtx_col, iq, is, name, set_name, & 971 comm) 972 type(scalapack), intent(in) :: scal 973 complex(DPC), intent(inout) :: matrix(:,:) !< (local_row,local_col) 974 integer, intent(in) :: nmtx, nmtx_col ! global row and col of the matrix to read in 975 integer, intent(in) :: iq 976 integer, intent(in) :: is 977 character(len=*), intent(in) :: name 978 character(len=*), intent(in) :: set_name ! where to read from 979 integer, intent(in) :: comm ! MPI communicator 980 981 integer :: ii, jj, error, rank 982 real(DP), allocatable :: data(:,:,:,:,:,:) 983 984 integer(HID_T) :: file_id ! File identifier 985 integer(HID_T) :: dset_id ! Dataset identifier 986 integer(HID_T) :: filespace ! Dataspace identifier in file 987 integer(HID_T) :: memspace ! Dataspace identifier in mem 988 integer(HID_T) :: plist_id ! Property list identifier for parallel IO 989 990 integer(HSIZE_T) :: count(6), countm(6), offset(6), offsetm(6), stride(6), block(6) 991 integer(HSIZE_T) :: countr(6), offsetr(6), strider(6), blockr(6) 992 993 integer :: info, rowremainder, colremainder 994 995 PUSH_SUB(read_eps_matrix_par_hdf5_general) 996 997 ! MPI_INFO_NULL can be used for info in place of of any actual information 998 ! (also when using subgroups, hopefully) 999 info = MPI_INFO_NULL 1000 1001 ! synchronize subgroup 1002 call MPI_barrier(comm,mpierr) 1003 1004 rank=6 1005 countm(1) = 2 1006 countm(2) = scal%npr 1007 countm(3) = scal%npc 1008 countm(4) = 1 1009 countm(5) = 1 1010 countm(6) = 1 1011 1012 offsetm(:) = 0 1013 1014 count(1) = 1 1015 count(2) = scal%npr/scal%nbl 1016 count(3) = scal%npc/scal%nbl 1017 count(4) = 1 1018 count(5) = 1 1019 count(6) = 1 1020 1021 block(1) = 2 1022 block(2) = scal%nbl 1023 block(3) = scal%nbl 1024 block(4) = 1 1025 block(5) = 1 1026 block(6) = 1 1027 1028 offset(1) = 0 1029 offset(2) = scal%myprow*scal%nbl 1030 offset(3) = scal%mypcol*scal%nbl 1031 offset(4) = 0 1032 offset(5) = is-1 1033 offset(6) = iq-1 1034 1035 stride(1) = 1 1036 stride(2) = scal%nprow*scal%nbl 1037 stride(3) = scal%npcol*scal%nbl 1038 stride(4) = 1 1039 stride(5) = 1 1040 stride(6) = 1 1041 1042 SAFE_ALLOCATE(data,(countm(1),countm(2),countm(3),countm(4),countm(5),countm(6))) 1043 data = 0.0d0 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 h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error, access_prp = plist_id) 1049 call h5pclose_f(plist_id,error) 1050 1051 call h5dopen_f(file_id, trim(set_name), dset_id, error) 1052 call h5dget_space_f(dset_id,filespace,error) 1053 1054 call h5screate_simple_f(rank, countm, memspace, error) 1055 !XXX call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, countm, error) 1056 1057 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error, stride, block) 1058 1059 ! Bottom Rows 1060 rowremainder = mod(scal%npr,scal%nbl) 1061 if (rowremainder .ne. 0) then 1062 offsetr=offset 1063 countr=count 1064 blockr=block 1065 strider=stride 1066 offsetr(2)=nmtx-rowremainder 1067 countr(2)=rowremainder 1068 blockr(2)=1 1069 strider(2)=1 1070 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1071 !write(6,*) peinf%inode, "I have the bottom row", rowremainder, scal%npc 1072 endif 1073 1074 ! Right Columns 1075 colremainder = mod(scal%npc,scal%nbl) 1076 if (colremainder .ne. 0) then 1077 offsetr=offset 1078 countr=count 1079 blockr=block 1080 strider=stride 1081 offsetr(3)=nmtx_col - colremainder 1082 countr(3)=colremainder 1083 blockr(3)=1 1084 strider(3)=1 1085 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1086 !write(6,*) peinf%inode, "I have the right column", colremainder, scal%npr 1087 ! Bottom Corner of Matrix 1088 if (rowremainder .ne. 0) then 1089 offsetr=offset 1090 countr=count 1091 blockr=block 1092 strider=stride 1093 offsetr(2)=nmtx-rowremainder 1094 countr(2)=rowremainder 1095 blockr(2)=1 1096 strider(2)=1 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 bottom both" 1103 endif 1104 endif 1105 1106 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 1107 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 1108 1109 call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, filespace, & 1110 xfer_prp = plist_id) 1111 1112 do jj = 1, scal%npc 1113 do ii = 1, scal%npr 1114 matrix(ii,jj) = CMPLX(data(1,ii,jj,1,1,1),data(2,ii,jj,1,1,1)) 1115 enddo 1116 enddo 1117 1118 SAFE_DEALLOCATE(data) 1119 call h5dclose_f(dset_id, error) 1120 call h5sclose_f(memspace, error) 1121 call h5sclose_f(filespace, error) 1122 call h5fclose_f(file_id, error) 1123 1124 POP_SUB(read_eps_matrix_par_hdf5_general) 1125 1126end subroutine read_eps_matrix_par_hdf5_general 1127 1128subroutine read_eps_matrix_par_hdf5_general_f(scal, matrix, nmtx, nmtx_col, Nfreq, iq, is, name, set_name, & 1129 comm) 1130 type(scalapack), intent(in) :: scal 1131 complex(DPC), intent(inout) :: matrix(:,:,:) !< (local_row,local_col,Nfreq) 1132 integer, intent(in) :: nmtx, nmtx_col ! global row and col of the matrix to read in 1133 integer, intent(in) :: nfreq 1134 integer, intent(in) :: iq 1135 integer, intent(in) :: is 1136 character(len=*), intent(in) :: name 1137 character(len=*), intent(in) :: set_name ! where to read from 1138 integer, intent(in) :: comm ! MPI communicator 1139 1140 integer :: ii, jj, kk, error, rank 1141 real(DP), allocatable :: data(:,:,:,:,:,:) 1142 1143 integer(HID_T) :: file_id ! File identifier 1144 integer(HID_T) :: dset_id ! Dataset identifier 1145 integer(HID_T) :: filespace ! Dataspace identifier in file 1146 integer(HID_T) :: memspace ! Dataspace identifier in mem 1147 integer(HID_T) :: plist_id ! Property list identifier for parallel IO 1148 1149 integer(HSIZE_T) :: count(6), countm(6), offset(6), offsetm(6), stride(6), block(6) 1150 integer(HSIZE_T) :: countr(6), offsetr(6), strider(6), blockr(6) 1151 1152 integer :: info, rowremainder, colremainder 1153 1154 PUSH_SUB(read_eps_matrix_par_hdf5_general_f) 1155 1156 ! even if possible lets keep it simple for now and assume nmtx = nmtx_col 1157 ! meaning that we consider only square subspace matrices 1158 1159 ! MPI_INFO_NULL can be used for info in place of of any actual information 1160 ! (also when using subgroups, hopefully) 1161 info = MPI_INFO_NULL 1162 1163 ! synchronize subgroup 1164 call MPI_barrier(comm,mpierr) 1165 1166 rank=6 1167 countm(1) = 2 1168 countm(2) = scal%npr 1169 countm(3) = scal%npc 1170 countm(4) = Nfreq 1171 countm(5) = 1 1172 countm(6) = 1 1173 1174 offsetm(:) = 0 1175 1176 count(1) = 1 1177 count(2) = scal%npr/scal%nbl 1178 count(3) = scal%npc/scal%nbl 1179 count(4) = 1 1180 count(5) = 1 1181 count(6) = 1 1182 1183 block(1) = 2 1184 block(2) = scal%nbl 1185 block(3) = scal%nbl 1186 block(4) = Nfreq 1187 block(5) = 1 1188 block(6) = 1 1189 1190 offset(1) = 0 1191 offset(2) = scal%myprow*scal%nbl 1192 offset(3) = scal%mypcol*scal%nbl 1193 offset(4) = 0 1194 offset(5) = is-1 1195 offset(6) = iq-1 1196 1197 stride(1) = 1 1198 stride(2) = scal%nprow*scal%nbl 1199 stride(3) = scal%npcol*scal%nbl 1200 stride(4) = 1 1201 stride(5) = 1 1202 stride(6) = 1 1203 1204 SAFE_ALLOCATE(data,(countm(1),countm(2),countm(3),countm(4),countm(5),countm(6))) 1205 data = 0.0d0 1206 1207 call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error) 1208 call h5pset_fapl_mpio_f(plist_id, comm, info, error) 1209 1210 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error, access_prp = plist_id) 1211 call h5pclose_f(plist_id,error) 1212 1213 call h5dopen_f(file_id, trim(set_name), dset_id, error) 1214 call h5dget_space_f(dset_id,filespace,error) 1215 1216 call h5screate_simple_f(rank, countm, memspace, error) 1217 !XXX call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offsetm, countm, error) 1218 1219 call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error, stride, block) 1220 1221 ! Bottom Rows 1222 rowremainder = mod(scal%npr,scal%nbl) 1223 if (rowremainder .ne. 0) then 1224 offsetr=offset 1225 countr=count 1226 blockr=block 1227 strider=stride 1228 offsetr(2)=nmtx-rowremainder 1229 countr(2)=rowremainder 1230 blockr(2)=1 1231 strider(2)=1 1232 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1233 !write(6,*) peinf%inode, "I have the bottom row", rowremainder, scal%npc 1234 endif 1235 1236 ! Right Columns 1237 colremainder = mod(scal%npc,scal%nbl) 1238 if (colremainder .ne. 0) then 1239 offsetr=offset 1240 countr=count 1241 blockr=block 1242 strider=stride 1243 offsetr(3)=nmtx-colremainder 1244 countr(3)=colremainder 1245 blockr(3)=1 1246 strider(3)=1 1247 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1248 !write(6,*) peinf%inode, "I have the right column", colremainder, scal%npr 1249 ! Bottom Corner of Matrix 1250 if (rowremainder .ne. 0) then 1251 offsetr=offset 1252 countr=count 1253 blockr=block 1254 strider=stride 1255 offsetr(2)=nmtx-rowremainder 1256 countr(2)=rowremainder 1257 blockr(2)=1 1258 strider(2)=1 1259 offsetr(3)=nmtx-colremainder 1260 countr(3)=colremainder 1261 blockr(3)=1 1262 strider(3)=1 1263 call h5sselect_hyperslab_f(filespace, H5S_SELECT_OR_F, offsetr, countr, error, strider, blockr) 1264 !write(6,*) peinf%inode, "I have bottom both" 1265 endif 1266 endif 1267 1268 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 1269 call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error) 1270 1271 call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, data, countm, error, memspace, filespace, & 1272 xfer_prp = plist_id) 1273 1274 do kk =1, Nfreq 1275 do jj = 1, scal%npc 1276 do ii = 1, scal%npr 1277 matrix(ii,jj,kk) = CMPLX(data(1,ii,jj,kk,1,1),data(2,ii,jj,kk,1,1)) 1278 enddo 1279 enddo 1280 end do 1281 1282 SAFE_DEALLOCATE(data) 1283 call h5dclose_f(dset_id, error) 1284 call h5sclose_f(memspace, error) 1285 call h5sclose_f(filespace, error) 1286 call h5fclose_f(file_id, error) 1287 1288 POP_SUB(read_eps_matrix_par_hdf5_general_f) 1289 1290end subroutine read_eps_matrix_par_hdf5_general_f 1291 1292#endif 1293! END SUBSPACE specific routine 1294!========================================================================================== 1295 1296#endif 1297end module epsread_hdf5_m 1298