rendered paste body//$ mpicc mpitest.c -o mpitest -std=c99 -lhdf5 -DMPI_DEBUG//$ NDOMAINS=2 mpirun -np 4 ./mpitest// http://www.hdfgroup.org/ftp/HDF5/examples/examples-by-api/hdf5-examples/1_8/C/H5D/#include <stdlib.h>#include <stdbool.h>#include <math.h>#include <assert.h>#include <unistd.h>#include <stdio.h>#include <mpi.h>#include <hdf5.h>typedef struct{ int rank; int size; int *members; MPI_Comm comm; MPI_Group group; void *parent;} MPI_Object;typedef struct{ int start; int stop; int step;} slice;int take(slice s, int idx){ int count = s.stop-s.start; assert(0 <= idx && idx < count); return s.start + idx*s.step;}// ------------------------------------------------------------------const char datfile[] = "parallelc.hdf5";const char logfile[] = "parallelc.log";static int devnull;#define H5RUN(x,y) if((x=y) < 0) { fprintf(stderr, "H5 %s:%d returned %d.\n", __FILE__, __LINE__, x); return -1; }#define MPIRUN(x,y) if((x=y) != MPI_SUCCESS) { fprintf(stderr, "MPI %s:%d returned %d.\n", __FILE__, __LINE__, x); return -1; }#define MPI_WAIT_ORDERED 10000#define MPI_PRE_ORDERED(comm_obj) \ for(int i=0; i<(comm_obj).rank; i++) { \ MPI_Barrier((comm_obj).comm); \ usleep(MPI_WAIT_ORDERED); \ }#define MPI_POST_ORDERED(comm_obj) \ for(int i=(comm_obj).rank; i<(comm_obj).size; i++) { \ usleep(MPI_WAIT_ORDERED); \ MPI_Barrier((comm_obj).comm); \ }// ------------------------------------------------------------------typedef struct{ int nbands; int mynbands; MPI_Object comm;} BandDescriptor;BandDescriptor new_band_descriptor(MPI_Object band_comm){ BandDescriptor bd; bd.nbands = 24; // Determine band descriptors (from band_descriptor.py) bd.mynbands = bd.nbands / band_comm.size; assert(bd.mynbands * band_comm.size == bd.nbands); bd.comm = band_comm; return bd;}// ------------------------------------------------------------------typedef struct{ int N_c[3]; // global number of grid points per axis incl. padding bool pbc_c[3]; int parsize_c[3]; int parpos_c[3]; double h_c[3]; int beg_c[3]; int end_c[3]; int n_c[3]; // local number of grid points per axis int gpts_c[3]; // global number of grid points per axis MPI_Object comm;} GridDescriptor;GridDescriptor new_grid_descriptor(MPI_Object domain_comm){ GridDescriptor gd; gd.N_c[0] = gd.N_c[1] = gd.N_c[2] = 12; gd.pbc_c[0] = gd.pbc_c[1] = gd.pbc_c[2] = false; // Determine domain decomposition (from domain.py) gd.parsize_c[0] = domain_comm.size; gd.parsize_c[1] = gd.parsize_c[2] = 1; assert(gd.parsize_c[0] * gd.parsize_c[1] * gd.parsize_c[2] == domain_comm.size); gd.parpos_c[0] = domain_comm.rank / (gd.parsize_c[1] * gd.parsize_c[2]); gd.parpos_c[1] = (domain_comm.rank % (gd.parsize_c[1] * gd.parsize_c[2])) / gd.parsize_c[2]; gd.parpos_c[2] = domain_comm.rank % gd.parsize_c[2]; assert(gd.parpos_c[0] * gd.parsize_c[1] * gd.parsize_c[2] \ + gd.parpos_c[1] * gd.parsize_c[2] \ + gd.parpos_c[2] == domain_comm.rank); // Determine grid descriptors (from grid_descriptor.py) gd.h_c[0] = gd.h_c[1] = gd.h_c[2] = 0.2; for(int c=0; c<3; c++) { double tmp; int n = gd.parsize_c[c]+1, *p = malloc(n*sizeof(int)); for(int i=0; i<n; i++) { tmp = i * gd.N_c[c]; tmp /= gd.parsize_c[c]; tmp += 0.4999; if(i==0 && !gd.pbc_c[c]) p[0] = 1; else p[i] = tmp; //printf("c=%d, i=%d, p=%d\n", c, i, p[i]); } gd.beg_c[c] = p[gd.parpos_c[c]]; gd.end_c[c] = p[gd.parpos_c[c] + 1]; gd.n_c[c] = gd.end_c[c] - gd.beg_c[c]; gd.gpts_c[c] = gd.N_c[c] - 1 + gd.pbc_c[c]; free(p); }#ifdef VERBOSE // Print ordered grid information MPI_Object world = *((MPI_Object*) domain_comm.parent); MPI_PRE_ORDERED(world); int ndigits = 1 + floor(log10(world.size)); printf("GRID[%0*d]: gd: %d/%d, pbc=%d,%d,%d; parsize=%d,%d,%d; parpos=%d,%d,%d; beg=%2d,%2d,%2d; end=%2d,%2d,%2d; n=%2d,%2d,%2d\n", ndigits, world.rank, domain_comm.rank, domain_comm.size, gd.pbc_c[0], gd.pbc_c[1], gd.pbc_c[2], gd.parsize_c[0], gd.parsize_c[1], gd.parsize_c[2], gd.parpos_c[0], gd.parpos_c[1], gd.parpos_c[2], gd.beg_c[0], gd.beg_c[1], gd.beg_c[2], gd.end_c[0], gd.end_c[1], gd.end_c[2], gd.n_c[0], gd.n_c[1], gd.n_c[2]); fflush(stdout); MPI_POST_ORDERED(world);#endif gd.comm = domain_comm; return gd;}// ------------------------------------------------------------------int check_fspace(MPI_Object world, hid_t fspace_id, BandDescriptor bd, GridDescriptor gd, int ndims, const hsize_t* refsize){ htri_t valid, simple; H5RUN(valid, H5Sselect_valid(fspace_id)); H5RUN(simple, H5Sis_simple(fspace_id)); assert((int) H5Sget_simple_extent_ndims(fspace_id) == ndims); hsize_t *count = (hsize_t*) malloc(ndims*sizeof(hsize_t)); H5RUN(devnull, H5Sget_simple_extent_dims(fspace_id, (hsize_t*) count, (hsize_t*) NULL)); int npts = H5Sget_select_npoints(fspace_id);#ifdef VERBOSE char tbuf[32]; H5S_sel_type ftype; H5RUN(ftype, H5Sget_select_type(fspace_id)); switch(ftype) { case H5S_SEL_ERROR: strncpy(tbuf, "ERROR", 32); break; case H5S_SEL_NONE: strncpy(tbuf, "NONE", 32); break; case H5S_SEL_POINTS: strncpy(tbuf, "POINTS", 32); break; case H5S_SEL_HYPERSLABS: strncpy(tbuf, "HYPERSLABS", 32); break; case H5S_SEL_ALL: strncpy(tbuf, "ALL", 32); break; case H5S_SEL_N: strncpy(tbuf, "N", 32); break; default: strncpy(tbuf, "UNKNOWN", 32); } char cbuf[32]; switch(ndims) { case 3: snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]); break; case 4: snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]); break; default: strncpy(cbuf, "UNKNOWN", 32); } MPI_PRE_ORDERED(world); int ndigits = 1 + floor(log10(world.size)); printf("FILE[%0*d]: valid: %d, simple: %d, npts: %d, H5S: %s, count: %s\n", ndigits, world.rank, valid, simple, npts, tbuf, cbuf); MPI_POST_ORDERED(world);#endif assert((int) valid>0); assert((int) simple>0); int refnpts = 1; for(int i=0; i<ndims; i++) { assert(count[i] == refsize[i]); refnpts *= refsize[i]; } //assert(npts == refnpts); //XXX BAD XXX get selected extent instead!!! free(count); return 0;}// ------------------------------------------------------------------int check_mspace(MPI_Object world, hid_t mspace_id, BandDescriptor bd, GridDescriptor gd, int ndims, const hsize_t* refsize){ htri_t valid, simple; H5RUN(valid, H5Sselect_valid(mspace_id)); H5RUN(simple, H5Sis_simple(mspace_id)); assert((int) H5Sget_simple_extent_ndims(mspace_id) == ndims); hsize_t *count = (hsize_t*) malloc(ndims*sizeof(hsize_t)); H5RUN(devnull, H5Sget_simple_extent_dims(mspace_id, (hsize_t*) count, (hsize_t*) NULL)); int npts = H5Sget_select_npoints(mspace_id);#ifdef VERBOSE char tbuf[32]; H5S_sel_type mtype; H5RUN(mtype, H5Sget_select_type(mspace_id)); switch(mtype) { case H5S_SEL_ERROR: strncpy(tbuf, "ERROR", 32); break; case H5S_SEL_NONE: strncpy(tbuf, "NONE", 32); break; case H5S_SEL_POINTS: strncpy(tbuf, "POINTS", 32); break; case H5S_SEL_HYPERSLABS: strncpy(tbuf, "HYPERSLABS", 32); break; case H5S_SEL_ALL: strncpy(tbuf, "ALL", 32); break; case H5S_SEL_N: strncpy(tbuf, "N", 32); break; default: strncpy(tbuf, "UNKNOWN", 32); } char cbuf[32]; switch(ndims) { case 3: snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]); break; case 4: snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]); break; default: strncpy(cbuf, "UNKNOWN", 32); } MPI_PRE_ORDERED(world); int ndigits = 1 + floor(log10(world.size)); printf("MEMO[%0*d]: valid: %d, simple: %d, npts: %d, H5S: %s, count: %s\n", ndigits, world.rank, valid, simple, npts, tbuf, cbuf); fflush(stdout); MPI_POST_ORDERED(world);#endif assert((int) valid>0); assert((int) simple>0); int refnpts = 1; for(int i=0; i<ndims; i++) { assert(count[i] == refsize[i]); refnpts *= refsize[i]; } //assert(npts == refnpts); //XXX BAD XXX can be zero if none selected!!! free(count); return 0;}// ------------------------------------------------------------------int check_dset(MPI_Object world, hid_t dset_id, BandDescriptor bd, GridDescriptor gd, int ndims, const hsize_t* count, int typesize){ int nbytes = H5Dget_storage_size(dset_id), refnbytes = typesize; for(int i=0; i<ndims; i++) { refnbytes *= count[i]; }#ifdef VERBOSE hid_t dcpl_id; H5RUN(dcpl_id, H5Dget_create_plist(dset_id)); char abuf[32]; H5D_alloc_time_t alloctime; H5RUN(devnull, H5Pget_alloc_time(dcpl_id, &alloctime)); switch(alloctime) { case H5D_ALLOC_TIME_ERROR: strncpy(abuf, "ERROR", 32); break; case H5D_ALLOC_TIME_DEFAULT: strncpy(abuf, "DEFAULT", 32); break; case H5D_ALLOC_TIME_EARLY: strncpy(abuf, "EARLY", 32); break; case H5D_ALLOC_TIME_LATE: strncpy(abuf, "LATE", 32); break; case H5D_ALLOC_TIME_INCR: strncpy(abuf, "INCR", 32); break; default: strncpy(abuf, "UNKNOWN", 32); } char fbuf[32]; H5D_fill_time_t filltime; H5RUN(devnull, H5Pget_fill_time(dcpl_id, &filltime)); switch(filltime) { case H5D_FILL_TIME_ERROR: strncpy(fbuf, "ERROR", 32); break; case H5D_FILL_TIME_ALLOC: strncpy(fbuf, "ALLOC", 32); break; case H5D_FILL_TIME_NEVER: strncpy(fbuf, "NEVER", 32); break; case H5D_FILL_TIME_IFSET: strncpy(fbuf, "IFSET", 32); break; default: strncpy(fbuf, "UNKNOWN", 32); } char lbuf[32]; H5D_layout_t layout; H5RUN(layout, H5Pget_layout(dcpl_id)); switch(layout) { case H5D_COMPACT: strncpy(lbuf, "COMPACT", 32); break; case H5D_CONTIGUOUS: strncpy(lbuf, "CONTIGUOUS", 32); break; case H5D_CHUNKED: strncpy(lbuf, "CHUNKED", 32); break; default: strncpy(lbuf, "UNKNOWN", 32); } H5Pclose(dcpl_id); char sbuf[32]; H5D_space_status_t stat; H5RUN(devnull, H5Dget_space_status(dset_id, &stat)); switch(stat) { case H5D_SPACE_STATUS_ERROR: strncpy(sbuf, "ERROR", 32); break; case H5D_SPACE_STATUS_NOT_ALLOCATED: strncpy(sbuf, "NOT_ALLOCATED", 32); break; case H5D_SPACE_STATUS_PART_ALLOCATED: strncpy(sbuf, "PART_ALLOCATED", 32); break; case H5D_SPACE_STATUS_ALLOCATED: strncpy(sbuf, "ALLOCATED", 32); break; default: strncpy(sbuf, "UNKNOWN", 32); } char cbuf[32]; switch(ndims) { case 3: snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]); break; case 4: snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]); break; default: strncpy(cbuf, "UNKNOWN", 32); } MPI_PRE_ORDERED(world); int ndigits = 1 + floor(log10(world.size)); printf("DSET[%0*d]: H5D: A=%s, F=%s, L=%s, S=%s, size: %d (%s*%d=%d)\n", ndigits, world.rank, abuf, fbuf, lbuf, sbuf, nbytes, cbuf, typesize, refnbytes); MPI_POST_ORDERED(world);#endif assert(nbytes == refnbytes); return 0;}// ------------------------------------------------------------------int run_write(MPI_Object world, MPI_Object domain_comm, MPI_Object band_comm){ // Print ordered processor distribution MPI_PRE_ORDERED(world); printf("WRITE: %d/%d, gd: %d/%d, bd: %d/%d\n", world.rank, world.size, domain_comm.rank, domain_comm.size, band_comm.rank, band_comm.size); fflush(stdout); MPI_POST_ORDERED(world); GridDescriptor gd = new_grid_descriptor(domain_comm); BandDescriptor bd = new_band_descriptor(band_comm); // Set up slices slice *indices = malloc(4*sizeof(slice)); indices[0].start = band_comm.rank * bd.mynbands; indices[0].stop = indices[0].start + bd.mynbands; indices[0].step = 1; for(int c=0; c<3; c++) { indices[c+1].start = gd.beg_c[c] - 1 + gd.pbc_c[c]; indices[c+1].stop = gd.end_c[c] - 1 + gd.pbc_c[c]; indices[c+1].step = 1; } // Fill in bogus but reproducible data in local domain int myngpts = gd.n_c[0] * gd.n_c[1] * gd.n_c[2]; double *psit_nG = malloc(bd.mynbands * myngpts * sizeof(double)); if(psit_nG) { double tmp, tmp2, r0_v[3]; for(int myn=0; myn<bd.mynbands; myn++) { int n = take(indices[0], myn), myg = 0; for(int c=0; c<3; c++) r0_v[c] = gd.pbc_c[c]?n*gd.h_c[c]*gd.N_c[c]:0; for(int gx=gd.beg_c[0]; gx<gd.end_c[0]; gx++) { tmp = pow(gd.h_c[0] * gx - r0_v[0], 2); for(int gy=gd.beg_c[1]; gy<gd.end_c[1]; gy++) { tmp2 = tmp + pow(gd.h_c[1] * gy - r0_v[1], 2); for(int gz=gd.beg_c[2]; gz<gd.end_c[2]; gz++) { psit_nG[myn*myngpts+myg] = sqrt(tmp2 + pow(gd.h_c[2] * gz - r0_v[2], 2)); myg++; } } } assert(myg==myngpts); } } // Fill in bogus but reproducible data on domain masters int ngpts = gd.gpts_c[0] * gd.gpts_c[1] * gd.gpts_c[2]; double *big_rhot_G = malloc(ngpts * sizeof(double)); if(big_rhot_G && domain_comm.rank == 0) { double tmp, tmp2, r0_v[3]; int g = 0; //XXX lazy global indexing should be offset by 1-pbc_c for(int c=0; c<3; c++) r0_v[c] = 0.5*gd.h_c[c]*gd.N_c[c]; for(int gx=0; gx<gd.gpts_c[0]; gx++) { tmp = exp(-pow(gd.h_c[0] * gx - r0_v[0], 2)); for(int gy=0; gy<gd.gpts_c[1]; gy++) { tmp2 = tmp * exp(-pow(gd.h_c[1] * gy - r0_v[1], 2)); for(int gz=0; gz<gd.gpts_c[2]; gz++) { big_rhot_G[g] = tmp2 * exp(-pow(gd.h_c[2] * gz - r0_v[2], 2)); g++; } } } assert(g==ngpts); } hid_t fapl_id = 0, f_id = 0, dxpl_id = 0, dcpl_id = 0; // File access property list (from io/hdf5.py) MPI_Info info; H5RUN(fapl_id, H5Pcreate(H5P_FILE_ACCESS));#ifndef MPI_DEBUG H5RUN(devnull, H5Pset_fclose_degree(fapl_id, H5F_CLOSE_STRONG));#endif if(world.size > 1) { // Copy MPI information to the file access property list (driver) char buf[20]; snprintf(buf, 20, "%d", world.size); MPIRUN(devnull, MPI_Info_create(&info)); MPIRUN(devnull, MPI_Info_set(info, "cb_nodes", buf)); //H5RUN(devnull, H5FD_mpio_init()); H5RUN(devnull, H5Pset_fapl_mpio(fapl_id, world.comm, info)); } else { H5RUN(devnull, H5Pset_fapl_stdio(fapl_id)); //H5RUN(devnull, H5Pset_fapl_log(fapl_id, logfile, H5FD_LOG_ALL, 256*1024)); } // Create file in write-mode using the file access property list from above H5RUN(f_id, H5Fcreate(datfile, H5F_ACC_TRUNC /*'w'*/, H5P_DEFAULT /*fcpl*/, fapl_id /*fapl*/)); // Speficy collective I/O for data transfer property list H5RUN(dxpl_id, H5Pcreate(H5P_DATASET_XFER)); if(world.size > 1) { H5RUN(devnull, H5Pset_dxpl_mpio(dxpl_id, H5FD_MPIO_COLLECTIVE)); //H5RUN(devnull, H5Pset_dxpl_mpio_chunk_opt(dxpl_id, H5FD_MPIO_CHUNK_MULTI_IO)); //new in 1.8x? //H5RUN(devnull, H5Pset_buffer(dxpl_id, 32*1024*1024)); //set transfer buffer to 32mb } // Specify early allocation for dataset creation property list if serial H5RUN(dcpl_id, H5Pcreate(H5P_DATASET_CREATE)); H5RUN(devnull, H5Pset_fill_time(dcpl_id, H5D_FILL_TIME_NEVER)); if(world.size == 1) { H5RUN(devnull, H5Pset_alloc_time(dcpl_id, H5D_ALLOC_TIME_EARLY)); } hid_t dset_id = 0, dtype_id = 0, fspace_id = 0, mspace_id = 0; if(psit_nG) { if(world.rank == 0) { printf("\npsit_nG\n"); fflush(stdout); } hsize_t count[] = {bd.mynbands, gd.n_c[0], gd.n_c[1], gd.n_c[2]}, start[] = {indices[0].start, indices[1].start, indices[2].start, indices[3].start}, fshape[] = {bd.nbands, gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}; // Create reduced-size memory dataspace H5RUN(mspace_id, H5Screate_simple(4, (const hsize_t*) count, NULL)); check_mspace(world, mspace_id, bd, gd, 4, (const hsize_t*) count); // Create fixed-size file dataspace H5RUN(fspace_id, H5Screate_simple(4, (const hsize_t*) fshape, (const hsize_t*) NULL)); // max_shape??? // Select the hyperslab that memory dataspace represents in file dataspace H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/)); //H5RUN(devnull, H5Sselect_none(fspace_id)); //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros)); check_fspace(world, fspace_id, bd, gd, 4, (const hsize_t*) fshape); // Create datatype and verify size and order H5RUN(dtype_id, H5Tcopy(H5T_NATIVE_DOUBLE)); H5RUN(devnull, H5Tset_order(dtype_id, H5T_ORDER_LE)); assert(H5Tget_size(dtype_id) == sizeof(double)); // Create psit_nG dataset H5RUN(dset_id, H5Dcreate(f_id, "PseudoWaveFunctions", dtype_id /*ftype*/, fspace_id /*fspace*/, dcpl_id /*dcpl*/)); check_dset(world, dset_id, bd, gd, 4, (const hsize_t*) fshape, sizeof(double)); H5RUN(devnull, H5Dwrite(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, psit_nG)); // Close all open ids H5RUN(devnull, H5Dclose(dset_id)); H5RUN(devnull, H5Tclose(dtype_id)); H5RUN(devnull, H5Sclose(fspace_id)); H5RUN(devnull, H5Sclose(mspace_id)); } if(big_rhot_G) { if(world.rank == 0) { printf("\nbig_rhot_G\n"); fflush(stdout); } hsize_t count[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}, start[] = {0, 0, 0}, fshape[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}; // Create reduced-size memory dataspace H5RUN(mspace_id, H5Screate_simple(3, (const hsize_t*) count, NULL)); if(domain_comm.rank != 0) { H5RUN(devnull, H5Sselect_none(mspace_id)); } check_mspace(world, mspace_id, bd, gd, 3, (const hsize_t*) count); // Create fixed-size file dataspace H5RUN(fspace_id, H5Screate_simple(3, (const hsize_t*) fshape, (const hsize_t*) NULL)); // max_shape??? // Select the hyperslab that memory dataspace represents in file dataspace if(domain_comm.rank == 0) { H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/)); } else { H5RUN(devnull, H5Sselect_none(fspace_id)); } //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros)); check_fspace(world, fspace_id, bd, gd, 3, (const hsize_t*) fshape); // Create datatype and verify size and order H5RUN(dtype_id, H5Tcopy(H5T_NATIVE_DOUBLE)); H5RUN(devnull, H5Tset_order(dtype_id, H5T_ORDER_LE)); assert(H5Tget_size(dtype_id) == sizeof(double)); // Create psit_nG dataset H5RUN(dset_id, H5Dcreate(f_id, "PseudoChargeDensity", dtype_id /*ftype*/, fspace_id /*fspace*/, dcpl_id /*dcpl*/)); check_dset(world, dset_id, bd, gd, 3, (const hsize_t*) fshape, sizeof(double)); H5RUN(devnull, H5Dwrite(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, big_rhot_G)); // Close all open ids H5RUN(devnull, H5Dclose(dset_id)); H5RUN(devnull, H5Tclose(dtype_id)); H5RUN(devnull, H5Sclose(fspace_id)); H5RUN(devnull, H5Sclose(mspace_id)); } // Close remaining ids H5RUN(devnull, H5Pclose(dcpl_id)); H5RUN(devnull, H5Pclose(dxpl_id)); H5RUN(devnull, H5Pclose(fapl_id)); if(world.size > 1) MPIRUN(devnull, MPI_Info_free(&info)); H5RUN(devnull, H5Fclose(f_id)); free(psit_nG); free(big_rhot_G); fflush(stdout); MPI_Barrier(world.comm); return 0;}// ------------------------------------------------------------------int run_read(MPI_Object world, MPI_Object domain_comm, MPI_Object band_comm){ // Print ordered processor distribution MPI_PRE_ORDERED(world); printf("READ: %d/%d, gd: %d/%d, bd: %d/%d\n", world.rank, world.size, domain_comm.rank, domain_comm.size, band_comm.rank, band_comm.size); fflush(stdout); MPI_POST_ORDERED(world); GridDescriptor gd = new_grid_descriptor(domain_comm); BandDescriptor bd = new_band_descriptor(band_comm); // Set up slices slice *indices = malloc(4*sizeof(slice)); indices[0].start = band_comm.rank * bd.mynbands; indices[0].stop = indices[0].start + bd.mynbands; indices[0].step = 1; for(int c=0; c<3; c++) { indices[c+1].start = gd.beg_c[c] - 1 + gd.pbc_c[c]; indices[c+1].stop = gd.end_c[c] - 1 + gd.pbc_c[c]; indices[c+1].step = 1; } // Allocate buffer for bogus but reproducible data in local domain int myngpts = gd.n_c[0] * gd.n_c[1] * gd.n_c[2]; double *psit_nG = malloc(bd.mynbands * myngpts * sizeof(double)); if(!psit_nG) return -1; // Allocate buffer for bogus but reproducible data on domain masters int ngpts = gd.gpts_c[0] * gd.gpts_c[1] * gd.gpts_c[2]; double *big_rhot_G = malloc(ngpts * sizeof(double)); if(!big_rhot_G) return -1; hid_t fapl_id = 0, f_id = 0, dxpl_id = 0; // File access property list (from io/hdf5.py) MPI_Info info; H5RUN(fapl_id, H5Pcreate(H5P_FILE_ACCESS));#ifndef MPI_DEBUG H5RUN(devnull, H5Pset_fclose_degree(fapl_id, H5F_CLOSE_STRONG));#endif if(world.size > 1) { // Copy MPI information to the file access property list (driver) char buf[20]; snprintf(buf, 20, "%d", world.size); MPIRUN(devnull, MPI_Info_create(&info)); MPIRUN(devnull, MPI_Info_set(info, "cb_nodes", buf)); //H5RUN(devnull, H5FD_mpio_init()); H5RUN(devnull, H5Pset_fapl_mpio(fapl_id, world.comm, info)); } else H5RUN(devnull, H5Pset_fapl_stdio(fapl_id)); //H5RUN(devnull, H5Pset_fapl_log(fapl_id, logfile, H5FD_LOG_ALL, 256*1024)); // Create file in write-mode using the file access property list from above H5RUN(f_id, H5Fopen(datfile, H5F_ACC_RDONLY /*'r'*/, fapl_id /*fapl*/)); // Speficy collective I/O for data transfer property list H5RUN(dxpl_id, H5Pcreate(H5P_DATASET_XFER)); if(world.size > 1) { H5RUN(devnull, H5Pset_dxpl_mpio(dxpl_id, H5FD_MPIO_COLLECTIVE)); //H5RUN(devnull, H5Pset_dxpl_mpio_chunk_opt(dxpl_id, H5FD_MPIO_CHUNK_MULTI_IO)); //new in 1.8x? //H5RUN(devnull, H5Pset_buffer(dxpl_id, 32*1024*1024)); //set transfer buffer to 32mb } hid_t dset_id = 0, dtype_id = 0, fspace_id = 0, mspace_id = 0; if(psit_nG) { if(world.rank == 0) { printf("\npsit_nG\n"); fflush(stdout); } hsize_t refsize[] = {bd.nbands, gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}, count[] = {bd.mynbands, gd.n_c[0], gd.n_c[1], gd.n_c[2]}, start[] = {indices[0].start, indices[1].start, indices[2].start, indices[3].start}; // Open psit_nG dataset H5RUN(dset_id, H5Dopen(f_id, "PseudoWaveFunctions")); check_dset(world, dset_id, bd, gd, 4, (const hsize_t*) refsize, sizeof(double)); // Open datatype and verify size and order H5RUN(dtype_id, H5Dget_type(dset_id)); assert((int) H5Tequal(dtype_id, H5T_NATIVE_DOUBLE) > 0); assert(H5Tget_size(dtype_id) == sizeof(double)); assert(H5Tget_order(dtype_id) == H5T_ORDER_LE); // Open file dataspace and verify fixed size H5RUN(fspace_id, H5Dget_space(dset_id)); assert(H5Sget_simple_extent_ndims(fspace_id) == 4); hsize_t fshape[4]; // Select the hyperslab that memory dataspace represents in file dataspace H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/)); //H5RUN(devnull, H5Sselect_none(fspace_id)); //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros)); check_fspace(world, fspace_id, bd, gd, 4, (const hsize_t*) refsize); // Create reduced-size memory dataspace H5RUN(mspace_id, H5Screate_simple(4, (const hsize_t*) count, NULL)); check_mspace(world, mspace_id, bd, gd, 4, (const hsize_t*) count); H5RUN(devnull, H5Dread(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, psit_nG)); // Close all open ids H5RUN(devnull, H5Dclose(dset_id)); H5RUN(devnull, H5Tclose(dtype_id)); H5RUN(devnull, H5Sclose(fspace_id)); H5RUN(devnull, H5Sclose(mspace_id)); } if(big_rhot_G) { if(world.rank == 0) { printf("\nbig_rhot_G\n"); fflush(stdout); } hsize_t refsize[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}, count[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]}, start[] = {0, 0, 0}; // Open psit_nG dataset H5RUN(dset_id, H5Dopen(f_id, "PseudoChargeDensity")); check_dset(world, dset_id, bd, gd, 3, (const hsize_t*) refsize, sizeof(double)); // Open datatype and verify size and order H5RUN(dtype_id, H5Dget_type(dset_id)); assert((int) H5Tequal(dtype_id, H5T_NATIVE_DOUBLE) > 0); assert(H5Tget_size(dtype_id) == sizeof(double)); assert(H5Tget_order(dtype_id) == H5T_ORDER_LE); // Open file dataspace and verify fixed size H5RUN(fspace_id, H5Dget_space(dset_id)); assert(H5Sget_simple_extent_ndims(fspace_id) == 3); hsize_t fshape[3]; // Select the hyperslab that memory dataspace represents in file dataspace if(domain_comm.rank == 0) { H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/)); } else { H5RUN(devnull, H5Sselect_none(fspace_id)); } //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros)); check_fspace(world, fspace_id, bd, gd, 3, (const hsize_t*) refsize); // Create reduced-size memory dataspace H5RUN(mspace_id, H5Screate_simple(3, (const hsize_t*) count, NULL)); if(domain_comm.rank != 0) { H5RUN(devnull, H5Sselect_none(mspace_id)); } check_mspace(world, mspace_id, bd, gd, 3, (const hsize_t*) count); H5RUN(devnull, H5Dread(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, big_rhot_G)); // Close all open ids H5RUN(devnull, H5Dclose(dset_id)); H5RUN(devnull, H5Tclose(dtype_id)); H5RUN(devnull, H5Sclose(fspace_id)); H5RUN(devnull, H5Sclose(mspace_id)); } H5RUN(devnull, H5Pclose(dxpl_id)); H5RUN(devnull, H5Pclose(fapl_id)); if(world.size > 1) MPIRUN(devnull, MPI_Info_free(&info)); H5RUN(devnull, H5Fclose(f_id)); if(psit_nG) { // Verify read data against bogus but reproducible equivalent double tmp, tmp2, r0_v[3]; for(int myn=0; myn<bd.mynbands; myn++) { int n = take(indices[0], myn), myg = 0; for(int c=0; c<3; c++) r0_v[c] = gd.pbc_c[c]?n*gd.h_c[c]*gd.N_c[c]:0; for(int gx=gd.beg_c[0]; gx<gd.end_c[0]; gx++) { tmp = pow(gd.h_c[0] * gx - r0_v[0], 2); for(int gy=gd.beg_c[1]; gy<gd.end_c[1]; gy++) { tmp2 = tmp + pow(gd.h_c[1] * gy - r0_v[1], 2); for(int gz=gd.beg_c[2]; gz<gd.end_c[2]; gz++) { assert(abs(psit_nG[myn*myngpts+myg] - sqrt(tmp2 + pow(gd.h_c[2] * gz - r0_v[2], 2))) < 1e-9); myg++; } } } assert(myg==myngpts); } free(psit_nG); } if(big_rhot_G) { MPIRUN(devnull, MPI_Bcast(big_rhot_G, ngpts, MPI_DOUBLE, 0, domain_comm.comm)); // Verify read data against bogus but reproducible equivalent double tmp, tmp2, r0_v[3]; int g = 0; //XXX lazy global indexing should be offset by 1-pbc_c for(int c=0; c<3; c++) r0_v[c] = 0.5*gd.h_c[c]*gd.N_c[c]; for(int gx=0; gx<gd.gpts_c[0]; gx++) { tmp = exp(-pow(gd.h_c[0] * gx - r0_v[0], 2)); for(int gy=0; gy<gd.gpts_c[1]; gy++) { tmp2 = tmp * exp(-pow(gd.h_c[1] * gy - r0_v[1], 2)); for(int gz=0; gz<gd.gpts_c[2]; gz++) { assert(abs(big_rhot_G[g] - tmp2 * exp(-pow(gd.h_c[2] * gz - r0_v[2], 2))) < 1e-9); g++; } } } assert(g==ngpts); free(big_rhot_G); } fflush(stdout); MPI_Barrier(world.comm); return 0;}// ------------------------------------------------------------------int main(int argc, char **argv){ int status; MPI_Init(&argc, &argv);#ifdef MPI_DEBUG // Default Errhandler is MPI_ERRORS_ARE_FATAL MPI_Errhandler_set(MPI_COMM_WORLD, MPI_ERRORS_RETURN);#endif MPI_Barrier(MPI_COMM_WORLD); // Basic info on world communicator int rank, nprocs, root; MPI_Object world; world.comm = MPI_COMM_WORLD; MPI_Comm_group(world.comm, &(world.group)); MPI_Comm_rank(world.comm, &(world.rank)); MPI_Comm_size(world.comm, &(world.size)); world.members = NULL; //XXX not needed, right? world.parent = NULL; //XXX not needed, right? // Determine communicator sizes rank = world.rank; nprocs = world.size; int ndomains = world.size, ngroups = 1; char *p = getenv("NDOMAINS"); if(p) { ndomains = atoi(p); if(ndomains <= 0 || nprocs < ndomains || nprocs % ndomains != 0) { if(rank == 0) fprintf(stderr, "Invalid number of domains (%d).\n", ndomains); MPI_Finalize(); return -1; } ngroups = nprocs / ndomains; } assert(ndomains * ngroups == nprocs); // Create domain communicator MPI_Object domain_comm; domain_comm.members = (int*) malloc(ndomains*sizeof(int)); root = rank - rank % ndomains; for(int i=0; i<ndomains; i++) domain_comm.members[i] = root + i; MPI_Group_incl(world.group, ndomains, domain_comm.members, &(domain_comm.group)); MPI_Comm_create(world.comm, domain_comm.group, &(domain_comm.comm)); // has a memory leak! MPI_Comm_rank(domain_comm.comm, &(domain_comm.rank)); MPI_Comm_size(domain_comm.comm, &(domain_comm.size)); domain_comm.parent = (void*) &world; // Create band communicator MPI_Object band_comm; band_comm.members = (int*) malloc(ngroups*sizeof(int)); root = rank % ndomains; for(int i=0; i<ngroups; i++) band_comm.members[i] = root + i * ndomains; MPI_Group_incl(world.group, ngroups, band_comm.members, &(band_comm.group)); MPI_Comm_create(world.comm, band_comm.group, &(band_comm.comm)); // has a memory leak! MPI_Comm_rank(band_comm.comm, &(band_comm.rank)); MPI_Comm_size(band_comm.comm, &(band_comm.size)); band_comm.parent = (void*) &world; bool write = true, read = true; p = getenv("ACTION"); if(p) { assert(strlen(p) == 1); switch(p[0]) { case 'w': case 'W': write = true; read = false; break; case 'r': case 'R': write = false; read = true; break; case 'b': case 'B': write = true; read = true; break; default: if(world.rank == 0) { fprintf(stderr, "\nWARNING: Action '%c' not understood.\n\n", p[0]); fflush(stderr); } } } if(write) status = run_write(world, domain_comm, band_comm); if(write && read && world.rank == 0) { printf("-----------------------------------------------\n"); fflush(stdout); } if(read) status = run_read(world, domain_comm, band_comm); free(domain_comm.members); free(band_comm.members); MPI_Finalize(); return status;}