diff options
Diffstat (limited to 'src/GetHyperslab.c')
| -rw-r--r-- | src/GetHyperslab.c | 198 |
1 files changed, 164 insertions, 34 deletions
diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c index b6e0af0..b1206e2 100644 --- a/src/GetHyperslab.c +++ b/src/GetHyperslab.c @@ -39,19 +39,12 @@ static int GetLocalHyperslab (const cGH *GH, int timelevel, int hdatatype, void *hdata); -/******************************************************************** - ******************** External Routines ************************ - ********************************************************************/ -/* Gerd's routine to get 1D lines from a 3D array */ -int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl, - const int *origin, - const int *directions, - int downsampling, - int length, - void **hdata, - int *hsize, - int proc); - +static int GetDiagonalFromFrom3D (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata); CCTK_INT Hyperslab_Get (const cGH *GH, @@ -274,27 +267,11 @@ static int GetLocalHyperslab (const cGH *GH, return (0); } - /* FIXME: hack for getting diagonals from 3D variables - This is calling Gerd's CollectData1D() routine which can - extract non-axis-parallel lines too but is fixed to 3D data. */ + /* diagonals from 3D variables are treated special */ if (mapping->is_diagonal_in_3D) { - const int origin[] = {0, 0, 0}; - const int directions[] = {1, 1, 1}; - int dummy_hsize; - void *dummy_hdata; - - - retval = Hyperslab_CollectData1D (GH, vindex, timelevel, origin, directions, - 1, -1, &dummy_hdata, &dummy_hsize, - mapping->target_proc); - if (! retval) - { - memcpy (hdata, dummy_hdata, - mapping->totals * CCTK_VarTypeSize (vinfo.vartype)); - free (dummy_hdata); - } - + retval = GetDiagonalFromFrom3D (GH, mapping, vindex, timelevel, hdatatype, + hdata); return (retval); } @@ -427,7 +404,7 @@ static int GetLocalHyperslab (const cGH *GH, /* get the byte pointer into the source array */ typed_vdata = (const char *) vdata + point[0]; -#if 1 +#if 0 fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_FullName (vindex)); #endif for (i = 1; i < vinfo.dim; i++) @@ -445,7 +422,7 @@ fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_F } else { -#if 1 +#if 0 fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vdata, typed_hdata); #endif memcpy (typed_hdata, typed_vdata, dim0_hsize); @@ -489,3 +466,156 @@ fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vda return (0); } + + +static int GetDiagonalFromFrom3D (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata) +{ + int i, j, k, myproc, nprocs, linear_idx; + CCTK_INT local_npoints, npoints; + int vdatatype, htypesize, vtypesize; + const char *vdata; + char *local_hdata; + const pGH *pughGH; + const pGExtras *extras; +#ifdef CCTK_MPI + CCTK_INT *global_npoints; + int *recvcnts, *displs; + MPI_Datatype mpidatatype; +#endif + + + + /* get the pGH pointer and the variable's GA structure */ + pughGH = PUGH_pGH (GH); + extras = ((const pGA *) pughGH->variables[vindex][timelevel])->extras; + + vdatatype = CCTK_VarTypeI (vindex); + if (hdatatype < 0) + { + hdatatype = vdatatype; + } + htypesize = CCTK_VarTypeSize (hdatatype); + vtypesize = CCTK_VarTypeSize (vdatatype); + vdata = (const char *) CCTK_VarDataPtrI (GH, timelevel, vindex); + + myproc = CCTK_MyProc (GH); + nprocs = CCTK_nProcs (GH); + if (nprocs == 1) + { + local_hdata = (char *) hdata; + } + else + { + local_hdata = (char *) malloc (mapping->global_hsize[0] * htypesize); + } + + i = mapping->global_startpoint[0] - extras->lb[myproc][0]; + j = mapping->global_startpoint[1] - extras->lb[myproc][1]; + k = mapping->global_startpoint[2] - extras->lb[myproc][2]; + local_npoints = 0; + for (npoints = 0; npoints < mapping->global_hsize[0]; npoints++) + { + if (i >= extras->ownership[0][0][0] && i < extras->ownership[0][1][0] && + j >= extras->ownership[0][0][1] && j < extras->ownership[0][1][1] && + k >= extras->ownership[0][0][2] && k < extras->ownership[0][1][2]) + { + linear_idx = i + j*extras->hyper_volume[1] + k*extras->hyper_volume[2]; + if (vdatatype != hdatatype) + { + mapping->conversion_fn (vdata + linear_idx*vtypesize, local_hdata, + 1, 1, 1); + } + else + { + memcpy (local_hdata, vdata + linear_idx*vtypesize, htypesize); + } + local_hdata += htypesize; + local_npoints++; + } + i += mapping->downsample[0]; + j += mapping->downsample[1]; + k += mapping->downsample[2]; + } + local_hdata -= local_npoints * htypesize; + +#ifdef CCTK_MPI + if (nprocs > 1) + { + /* allocate communication buffers */ + myproc = CCTK_MyProc (GH); + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + global_npoints = (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT)); + recvcnts = (int *) malloc (2 * nprocs * sizeof (int)); + displs = recvcnts + nprocs; + } + else + { + global_npoints = NULL; recvcnts = displs = NULL; hdata = NULL; + } + + /* gather the number of local points on each processor */ + if (mapping->target_proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allgather (&local_npoints, 1, PUGH_MPI_INT, + global_npoints, 1, PUGH_MPI_INT, + pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Gather (&local_npoints, 1, PUGH_MPI_INT, + global_npoints, 1, PUGH_MPI_INT, + mapping->target_proc, + pughGH->PUGH_COMM_WORLD)); + } + + /* compute the receive count and displacement vectors */ + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + for (i = 0; i < nprocs; i++) + { + recvcnts[i] = (int) global_npoints[i]; + displs[i] = i == 0 ? 0 : displs[i-1] + recvcnts[i-1]; + } + } + + /* get the MPI datatype for the hyperslab data */ + mpidatatype = PUGH_MPIDataType (pughGH, hdatatype); + + /* gather the local hyperslab data from each processor */ + if (mapping->target_proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allgatherv (local_hdata, local_npoints, + mpidatatype, hdata, recvcnts, displs, + mpidatatype, pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Gatherv (local_hdata, local_npoints, + mpidatatype, hdata, recvcnts, displs, + mpidatatype, mapping->target_proc, + pughGH->PUGH_COMM_WORLD)); + } + + /* free communication buffers */ + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + free (global_npoints); + free (recvcnts); + } + } +#endif /* CCTK_MPI */ + + /* free allocated resources */ + if (nprocs > 1) + { + free (local_hdata); + } + + return (0); +} |