From 6c9ec7931a8cceb3b41ace899499f519702f0130 Mon Sep 17 00:00:00 2001 From: tradke Date: Tue, 4 Dec 2001 21:47:26 +0000 Subject: Added (still incomplete) implementation of the new Hyperslab API. Not really usable yet, only single-processor case works. More changes to follow. git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHSlab/trunk@64 10716dce-81a3-4424-a2c8-48026a0d3035 --- src/CollectData1D.c | 56 +---- src/GetHyperslab.c | 573 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/Mapping.c | 519 +++++++++++++++++++++++++++++++++++++++++++++++ src/PUGHSlab.h | 12 +- src/PUGHSlabi.h | 16 +- src/make.code.defn | 2 +- 6 files changed, 1120 insertions(+), 58 deletions(-) create mode 100644 src/GetHyperslab.c create mode 100644 src/Mapping.c diff --git a/src/CollectData1D.c b/src/CollectData1D.c index ea5a052..38b5f60 100644 --- a/src/CollectData1D.c +++ b/src/CollectData1D.c @@ -29,16 +29,16 @@ CCTK_FILEVERSION(CactusPUGH_PUGHSlab_CollectData1D_c) } /* local function prototypes */ -int GetLocalLine(cGH *GH, pGExtras *extras, +int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl, + const int *origin, const int *directions, + int downsample, int length, + void **hdata, int *hsize, int proc); +int GetLocalLine(const cGH *GH, pGExtras *extras, int vardim, const int *S_l, const int *u_l, int ds, int len_l, int *locstart, int *locend, int *nlocpoints); -#ifdef CCTK_MPI -static MPI_Datatype CCTKtoMPItype(pGH *pughGH, int cctktype); -#endif - -static int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, +static int CollectLocalData1D(const cGH *GH, int vindex, int vtimelvl, const int *origin, const int *directions, int downsample, int length, void **hdata, int *hsize) @@ -110,7 +110,7 @@ static int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, } -int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl, +int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl, const int *origin, const int *directions, int downsample, int length, void **hdata, int *hsize, int proc) @@ -126,7 +126,7 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl, #ifdef CCTK_MPI int iproc; - pGH *pughGH; + const pGH *pughGH; CCTK_INT tmp, *rem_dcount; MPI_Datatype vmpitype; int *recvcnts, *displs; @@ -157,7 +157,7 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl, rem_dcount = NULL; recvcnts = NULL; displs = NULL; - vmpitype= CCTKtoMPItype(pughGH, vtype); + vmpitype = PUGH_MPIDataType (pughGH, vtype); /* Send local number of data elements */ tmp = (CCTK_INT) loc_dcount; @@ -236,41 +236,5 @@ int Hyperslab_CollectData1D(cGH *GH, int vindex, int vtimelvl, *hsize = all_dcount; } -#ifdef DEBUGME - printf(" CollectData1D done \n"); -#endif - return(0); + return (0); } - - -#ifdef CCTK_MPI -static MPI_Datatype CCTKtoMPItype(pGH *pughGH, int cctktype) -{ - MPI_Datatype mpitype; - switch (cctktype) { - - case CCTK_VARIABLE_CHAR: - mpitype = PUGH_MPI_CHAR; - break; - - case CCTK_VARIABLE_INT: - mpitype = PUGH_MPI_INT; - break; - - case CCTK_VARIABLE_REAL: - mpitype = PUGH_MPI_REAL; - break; - - case CCTK_VARIABLE_COMPLEX: - mpitype = pughGH->PUGH_mpi_complex; - break; - - default: - CCTK_WARN (1, "Unsupported MPI variable type"); - mpitype = 0; - break; - } - - return(mpitype); -} -#endif diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c new file mode 100644 index 0000000..0c883cc --- /dev/null +++ b/src/GetHyperslab.c @@ -0,0 +1,573 @@ + /*@@ + @file GetHyperslab.c + @date Sun 2 Dec 2001 + @author Thomas Radke + @desc + Routines to extract hyperslabs from CCTK array variables + @enddesc + @version $Id$ + @@*/ + + +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "CactusPUGH/PUGH/src/include/pugh.h" +#include "PUGHSlab.h" +#include "PUGHSlabi.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusPUGH_PUGHSlab_GetHyperslab_c) + +/******************************************************************** + ******************** Macro Definitions ************************ + ********************************************************************/ +/* define this if you want debugging output */ +/* #define DEBUG 1 */ + + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +static int GetLocalHyperslab (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata); +static const char *checkParameters (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + 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); + + + +CCTK_INT Hyperslab_Get (const cGH *GH, + CCTK_INT mapping_handle, + CCTK_INT vindex, + CCTK_INT timelevel, + CCTK_INT hdatatype, + void *hdata) +{ + int retval; + hslab_mapping_t *mapping; + + + mapping = PUGHSlabi_GetMapping (mapping_handle); + if (mapping == NULL) + { + return (-1); + } + + /* check mapping consistency */ + /*** FIXME ***/ + + /* get the processor-local hyperslab */ + retval = GetLocalHyperslab (GH, mapping, vindex, timelevel, hdatatype, hdata); + + return (retval); +} + + +CCTK_INT Hyperslab_GetList (const cGH *GH, + CCTK_INT mapping_handle, + CCTK_INT num_arrays, + const CCTK_INT *vindices /* num_arrays */, + const CCTK_INT *timelevels /* num_arrays */, + const CCTK_INT *hdatatypes /* num_arrays */, + void *const *hdata /* num_arrays */) +{ + int i, retval; + + + retval = 0; + for (i = 0; i < num_arrays; i++) + { + if (Hyperslab_Get (GH, mapping_handle, vindices[i], + timelevels ? timelevels[i] : 0, + hdatatypes ? hdatatypes[i] : -1, + hdata[i]) == 0) + { + retval++; + } + } + + return (retval); +} + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +/*@@ + @routine GetLocalHyperslab + @date Fri May 12 2000 + @author Thomas Radke + @desc + Extract a hyperslab from the processor-local chunk + of a domain-decomposed Cactus array variable. + + This routine delivers the local hyperslab data + to be collected into a global hyperslab + by Hyperslab_GetHyperslab(). + IO methods can call this routine as well to collect the + local hyperslab data and output it in parallel. + @enddesc + + @calls PUGHSlab_GetDatatypeConversionFn + + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH * + @vio in + @endvar + @var vindex + @vdesc index of variable to get a hyperslab from + @vtype int + @vio in + @endvar + @var timelevel + @vdesc timelvl of variable to get a hyperslab from + @vtype int + @vio in + @endvar + @var hdim + @vdesc dimensionality of the requested hyperslab + @vtype int + @vio in + @endvar + @var hdatatype + @vdesc CCTK datatype of the requested hyperslab + @vtype int + @vio in + @endvar + @var conversion_fn + @vdesc pointer to a user-supplied data conversion function + @vtype t_hslabConversionFn + @vio in + @endvar + @var global_startpoint + @vdesc global coordinates of the hyperslab origin + @vtype const int[dimensions of vindex] + @vio in + @endvar + @var directions + @vdesc directions which span the hyperslab + @vtype const int[hdim times dimensions of vindex] + @vio in + @endvar + @var extents + @vdesc number of grid points to follow in each hyperslab direction + starting from origin + Negative values are taken as extents up to the grid boundaries. + @vtype const int[hdim] + @vio in + @endvar + @var downsample + @vdesc downsampling values for each hyperslab dimension + @vtype const int[hdim] + @vio in + @endvar + @var hdata + @vdesc pointer to store the address of the hyperslab data buffer + @vtype void ** + @vio out + @endvar + @var free_data + @vdesc address of flag which decides whether the returned data needs + to be freed or not + @vtype int * + @vio out + @endvar + @var hsize + @vdesc sizes of the (local) hyperslab data buffer in each dimension + @vtype int[hdim] + @vio out + @endvar + @var hsize_global + @vdesc sizes of the global hyperslab data buffer in each dimension + @vtype int[hdim] + @vio out + @endvar + @var hoffset_global + @vdesc if not NULL, array to save the offsets of the local hyperslab + into the global one for each dimension + @vtype int[hdim] + @vio out + @endvar +@@*/ +static int GetLocalHyperslab (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata) +{ + int *point; /* looper over hyperslab dimensions */ + int *startpoint, /* hyperslab's local start and endpoint */ + *endpoint; /* within the variable's grid dimensions */ + int *downsample; /* the downsample[] vector extended to vdim */ + int *points_per_dim; /* points per subvolume */ + int myproc; /* local processor ID */ + int i; /* general looper */ + int vdim; /* looper over all source dimensions */ + int vdata_size, /* size of one data point in bytes for */ + hdata_size; /* source and hyperslab data */ + int dim0_points; /* number of hyperslab points in dim 0 */ + int dim0_hsize; /* byte size of hyperslab points in dim 0 */ + const char *typed_vdata; /* byte pointers into source and */ + char *typed_hdata; /* hyperslab data arrays */ + const void *vdata; + int retval; /* the return value (0 for success) */ + cGroup vinfo; /* variable's group info */ + pGH *pughGH; /* pointer to the current pGH */ + pGA *GA; /* the variable's GA structure from PUGH */ + const char *errormsg; /* error message string */ + t_hslabConversionFn conversion_fn; + + + /* do some plausibility checks */ + errormsg = checkParameters (GH, mapping, vindex, timelevel, hdata); + + /* immediately return in case of errors */ + if (errormsg) + { + CCTK_WARN (1, errormsg); + return (-1); + } + + /* check if there's any data to extract */ + if (mapping->totals == 0) + { + return (0); + } + + /* get the info on the variable to extract a hyperslab from */ + CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo); + + /* 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. */ + 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); + } + + return (retval); + } + + /* if datatype conversion was requested + get the appropriate predefined datatype conversion routine + in case the user didn't supply one by her own */ + if (hdatatype < 0) + { + hdatatype = vinfo.vartype; + } + conversion_fn = mapping->conversion_fn; + if (vinfo.vartype != hdatatype) + { + if (conversion_fn == NULL) + { + conversion_fn = PUGHSlabi_GetDatatypeConversionFn (vinfo.vartype, + hdatatype); + if (! conversion_fn) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "No predefined PUGHSlab routine available to convert " + "'%s' into '%s'", CCTK_VarTypeName (vinfo.vartype), + CCTK_VarTypeName (hdatatype)); + return (-1); + } + } + } + else if (conversion_fn) + { + CCTK_WARN (8, "Datatype conversion routine supplied but no datatype " + "conversion requested. Ignoring conversion routine..."); + conversion_fn = NULL; + } + + /* allocate the temporary arrays */ + point = (int *) malloc (5 * vinfo.dim * sizeof (int)); + startpoint = point + 1*vinfo.dim; + endpoint = point + 2*vinfo.dim; + downsample = point + 3*vinfo.dim; + points_per_dim = point + 4*vinfo.dim; + + memcpy (startpoint, mapping->local_startpoint, vinfo.dim * sizeof (int)); + memcpy (endpoint, mapping->local_endpoint, vinfo.dim * sizeof (int)); + memcpy (downsample, mapping->downsample, vinfo.dim * sizeof (int)); + + /* get the pGH pointer and the variable's GA structure */ + pughGH = PUGH_pGH (GH); + GA = (pGA *) pughGH->variables[vindex][timelevel]; + + /* get the local processor ID */ + myproc = CCTK_MyProc (GH); + + /* nested loop over vinfo.dim dimensions */ + /* NOTE: the following code assumes startpoint[vdim] < endpoint[vdim] */ + vdata = CCTK_VarDataPtrI (GH, timelevel, vindex); + + if (mapping->full_hyperslab && conversion_fn == NULL) + { + memcpy (hdata, vdata, mapping->totals * CCTK_VarTypeSize (vinfo.vartype)); + } + else + { + /* get the byte size of a single data point + in the variable and hyperslab data array */ + vdata_size = CCTK_VarTypeSize (vinfo.vartype); + hdata_size = CCTK_VarTypeSize (hdatatype); + + typed_hdata = (char *) hdata; + + /* compute the points_per_dim[] vector */ + /* NOTE: this could be computed at startup and kept in a GH extension + once we have one for thorn PUGHSlab */ + points_per_dim[0] = 1; + for (vdim = 1; vdim < vinfo.dim; vdim++) + { + points_per_dim[vdim] = points_per_dim[vdim-1] * + GA->extras->lnsize[vdim-1]; + } + + /* get the number of hyperslab points in lowest dimension + and their size in bytes */ + dim0_points = (endpoint[0] - startpoint[0]) / downsample[0]; + if ((endpoint[0] - startpoint[0]) % downsample[0]) + { + dim0_points++; + } + dim0_hsize = dim0_points * hdata_size; + + /* transform the ranges into byte ranges */ + for (i = 0; i < vinfo.dim; i++) + { + startpoint[i] *= vdata_size; + endpoint[i] *= vdata_size; + downsample[i] *= vdata_size; + } + + /* initialize the index vector to the local startpoint */ + memcpy (point, startpoint, vinfo.dim * sizeof (point[0])); + + /* do the nested loops starting with the innermost */ + vdim = 1; + while (1) + { + /* check for end of current loop */ + if (vinfo.dim > 1 && point[vdim] >= endpoint[vdim]) + { + /* increment outermost loopers */ + for (vdim++; vdim < vinfo.dim; vdim++) + { + point[vdim] += downsample[vdim]; + if (point[vdim] < endpoint[vdim]) + { + break; + } + } + + /* done if beyond outermost loop */ + if (vdim >= vinfo.dim) + { + break; + } + + /* reset innermost loopers */ + for (vdim--; vdim > 0; vdim--) + { + point[vdim] = startpoint[vdim]; + } + vdim = 1; + } + + /* get the byte pointer into the source array */ + typed_vdata = (const char *) vdata + point[0]; +#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++) + { + typed_vdata += point[i] * points_per_dim[i]; + } + + /* copy the data in lowest dimension: if possible copy all data points + in a row otherwise do it one by one */ + if (downsample[0] == vdata_size) + { + if (mapping->conversion_fn) + { + mapping->conversion_fn (typed_vdata, typed_hdata, dim0_points, 1,1); + } + else + { +#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); + } + } + else + { + if (mapping->conversion_fn) + { + mapping->conversion_fn (typed_vdata, typed_hdata, dim0_points, + downsample[0], 1); + typed_vdata += downsample[0] * dim0_points; + } + else + { + for (i = 0; i < dim0_hsize; i += hdata_size) + { + memcpy (typed_hdata + i, typed_vdata, vdata_size); + typed_vdata += downsample[0]; + } + } + } + typed_hdata += dim0_hsize; + + if (vinfo.dim > 1) + { + /* increment current looper */ + point[vdim] += downsample[vdim]; + } + else + { + /* exit loop if hyperslab dim is only 1D */ + break; + } + + } /* end of nested loops over all dimensions */ + } /* end of branch extracting the hyperslab data */ + + /* free allocated temporary memory */ + free (point); + + return (0); +} + + +static const char *checkParameters (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + void *hdata) +{ + cGroup vinfo; /* variable's group info */ +#if 0 + int i, vdim; /* looper */ + int num_directions; /* number of non-zero directions */ +#endif + + + /* check the variable index and timelevel */ + if (vindex < 0 || vindex >= CCTK_NumVars ()) + { + return ("Invalid variable index"); + } + if (timelevel < 0 || timelevel >= CCTK_NumTimeLevelsFromVarI (vindex)) + { + return ("Invalid timelevel"); + } + +#if 0 + /* check the passed pointers */ + if (! mapping->global_origin || ! directions || ! extents || ! mapping->downsample || + ! hdata || ! hsize) + { + return ("NULL pointer(s) passed as parameters"); + } + + /* check the extent and downsample parameters */ + for (vdim = 0; vdim < hdim; vdim++) + { + if (extents[vdim] == 0) + { + return ("Invalid hyperslab extent parameters"); + } + if (mapping->downsample[vdim] <= 0) + { + return ( "Invalid hyperslab downsample parameters"); + } + } +#endif + + /* get the info on the variable to extract a hyperslab from */ + if (CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo) < 0) + { + return ("Couldn't get group info"); + } + + /* check the variable's grouptype */ + if (vinfo.grouptype != CCTK_GF && vinfo.grouptype != CCTK_ARRAY) + { + return ("Invalid variable group type"); + } + + /* check the hyperslab dimension */ + if (mapping->hdim <= 0 || mapping->hdim > vinfo.dim) + { + return ("Invalid hyperslab dimension"); + } + +#if 0 + /* check the direction(s) of the hyperslab */ + for (i = 0; i < mapping->hdim; i++) + { + for (vdim = 0, num_directions = 0; vdim < vinfo.dim; vdim++) + { + if (directions[i * vinfo.dim + vdim]) + { + num_directions++; + } + } + if (num_directions == 0) + { + return ("Given direction vector is a null vector"); + } + if (num_directions != 1) + { + return ("Given direction vector isn't orthogonal"); + } + } +#endif + + /* check if PUGH is active */ + if (! PUGH_pGH (GH)) + { + return ("No GH extension for PUGH found. Did you activate thorn PUGH ?"); + } + + return (NULL); +} diff --git a/src/Mapping.c b/src/Mapping.c new file mode 100644 index 0000000..b53b050 --- /dev/null +++ b/src/Mapping.c @@ -0,0 +1,519 @@ + /*@@ + @file Mapping.c + @date Sun 2 Dec 2001 + @author Thomas Radke + @desc + Routines to define hyperslab mappings + @enddesc + @version $Id$ + @@*/ + + +#include +#include + +#include "cctk.h" + +#include "CactusPUGH/PUGH/src/include/pugh.h" + +#include "PUGHSlab.h" +#include "PUGHSlabi.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusPUGH_PUGHSlab_Mapping_c) + +/******************************************************************** + ******************** Macro Definitions ************************ + ********************************************************************/ +/* define this if you want debugging output */ +/* #define DEBUG 1 */ + +/* macros for getting the minimum/maximum value */ +#ifdef MIN +#undef MIN +#endif +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#ifdef MAX +#undef MAX +#endif +#define MAX(x, y) ((x) > (y) ? (x) : (y)) + +/* shortcuts for the local/global start/endpoints on this processor */ +#define MY_LOCAL_SP(extras, istag, dim) (extras->ownership[istag][0][dim]) +#define MY_LOCAL_EP(extras, istag, dim) (extras->ownership[istag][1][dim]) +#define MY_GLOBAL_SP(extras, myproc, istag, dim) \ + (extras->lb[myproc][dim] + MY_LOCAL_SP (extras, istag, dim)) +#define MY_GLOBAL_EP(extras, myproc, istag, dim) \ + (extras->lb[myproc][dim] + MY_LOCAL_EP (extras, istag, dim)) + + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +static int nmapping_list = 0; +static hslab_mapping_t *mapping_list = NULL; + +static int checkFullHyperslab (const pGA *GA, + const CCTK_INT *origin, + const CCTK_INT *extent, + const hslab_mapping_t *mapping); + + +CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( + const cGH *GH, + CCTK_INT vindex, + CCTK_INT hdim, + const CCTK_INT *direction /* vdim*hdim */, + const CCTK_INT *origin /* vdim */, + const CCTK_INT *extent /* vdim */, + const CCTK_INT *downsample /* hdim */, + CCTK_INT table_handle, + CCTK_INT target_proc, + t_hslabConversionFn conversion_fn, + CCTK_INT *hsize /* hdim */) +{ + int i, j, vdim, num_dirs, stagger_type, retval; + int stagger_index; + int myproc; + hslab_mapping_t *mapping; + const char *error_msg; + pGH *pughGH; /* pointer to the current pGH */ + pGA *GA; /* the variable's GA structure from PUGH */ + + + /* PUGHSlab doesn't use table information */ + if (table_handle >= 0) + { + CCTK_WARN (1, "Hyperslab_DefineGlobalMappingByIndex: table information is " + "ignored"); + } + + /* check parameter consistency */ + retval = 0; + error_msg = NULL; + vdim = CCTK_GroupDimFromVarI (vindex); + if (vdim <= 0) + { + error_msg = "invalid variable index given"; + retval = -1; + } + else if (hdim < 0 || hdim > vdim) + { + error_msg = "invalid hyperslab dimension given"; + retval = -2; + } + else if (hsize == NULL) + { + error_msg = "NULL pointer passed for hyperslab size return parameter"; + retval = -3; + } + else if (target_proc >= CCTK_nProcs (GH)) + { + error_msg = "invalid target procesor ID given"; + retval = -4; + } + else + { + for (i = 0; i < vdim; i++) + { + retval |= origin[i] < 0 || extent[i] <= 0; + if (downsample && i < hdim) + { + retval |= downsample[i] <= 0; + } + } + if (retval) + { + error_msg = "invalid hyperslab origin/extent/downsample given"; + retval = -5; + } + } + if (! retval) + { + mapping = (hslab_mapping_t *) malloc (sizeof (hslab_mapping_t)); + if (mapping) + { + mapping->vectors = (int *) malloc ((6*vdim + (2+vdim)*hdim) * sizeof (int)); + } + if (mapping == NULL || mapping->vectors == NULL) + { + if (mapping) + { + free (mapping); + } + error_msg = "couldn't allocate hyperslab mapping structure"; + retval = -6; + } + } + + /* return in case of errors */ + if (retval) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Hyperslab_DefineGlobalMappingByIndex: %s", error_msg); + return (retval); + } + + /* assign memory for the other vectors */ + mapping->local_startpoint = mapping->vectors + 0*vdim; + mapping->local_endpoint = mapping->vectors + 1*vdim; + mapping->global_startpoint = mapping->vectors + 2*vdim; + mapping->global_endpoint = mapping->vectors + 3*vdim; + mapping->do_dir = mapping->vectors + 4*vdim; + mapping->downsample = mapping->vectors + 5*vdim; + mapping->local_hsize = mapping->vectors + 6*vdim + 0*hdim; + mapping->global_hsize = mapping->vectors + 6*vdim + 1*hdim; + mapping->directions = mapping->vectors + 6*vdim + 2*hdim; + + /* check direction vectors */ + for (j = 0; j < hdim; j++) + { + num_dirs = 0; + for (i = 0; i < vdim; i++) + { + if (direction[j*vdim + i]) + { + mapping->directions[i] = 1; + num_dirs++; + } + } + mapping->is_diagonal_in_3D = num_dirs == 3 && hdim == 1; + if (num_dirs != 1 && ! mapping->is_diagonal_in_3D) + { + free (mapping->vectors); + free (mapping); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Hyperslab_DefineGlobalMappingByIndex: %d-direction vector " + "isn't axis-orthogonal", j); + return (-7); + } + } + for (i = 0; i < vdim; i++) + { + mapping->do_dir[i] = 0; + for (j = 0; j < hdim; j++) + { + if (direction[j*vdim + i]) + { + mapping->do_dir[i]++; + } + } + if (mapping->do_dir[i] > 1) + { + free (mapping->vectors); + free (mapping); + CCTK_WARN (1, "Hyperslab_DefineGlobalMappingByIndex: duplicate direction " + "vectors given"); + return (-8); + } + } + + /* get the pGH pointer and the variable's GA structure */ + pughGH = PUGH_pGH (GH); + GA = (pGA *) pughGH->variables[vindex][0]; + myproc = CCTK_MyProc (GH); + stagger_type = CCTK_GroupStaggerIndexGI (CCTK_GroupIndexFromVarI (vindex)); + + /* check extent */ + for (i = j = 0; i < vdim; i++) + { + if (mapping->do_dir[i]) + { + if (origin[i] + extent[j] > GA->extras->nsize[i]) + { + free (mapping->vectors); + free (mapping); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Hyperslab_DefineGlobalMappingByIndex: extent in " + "%d-direction exceeds grid size", j); + return (-8); + } + j++; + } + } + + /* now fill out the hyperslab mapping structure */ + for (i = j = 0; i < vdim; i++) + { + mapping->downsample[i] = 1; + + if (mapping->do_dir[i]) + { + if (downsample) + { + mapping->downsample[i] = downsample[j]; + } + mapping->global_hsize[j] = extent[j] / mapping->downsample[i]; + if (extent[j] % mapping->downsample[i]) + { + mapping->global_hsize[j]++; + } + /* subtract ghostzones for periodic BC */ + if (GA->connectivity->perme[i]) + { + mapping->global_hsize[j] -= 2 * GA->extras->nghostzones[i]; + } + j++; + } + } + + /* check whether the full local data patch was requested as hyperslab */ + mapping->full_hyperslab = checkFullHyperslab (GA, origin, extent, mapping); + if (mapping->full_hyperslab) + { + memset (mapping->local_startpoint, 0, vdim * sizeof (int)); + memcpy (mapping->local_endpoint, GA->extras->lnsize, vdim * sizeof (int)); + mapping->totals = GA->extras->npoints; + } + else + { + /* compute the global endpoint */ + for (i = j = 0; i < vdim; i++) + { + mapping->global_endpoint[i] = origin[i] + + mapping->do_dir[i] ? extent[j++] : 1; + } + + /* compute this processor's global startpoint from the global ranges */ + for (i = 0; i < vdim; i++) + { + stagger_index = CCTK_StaggerDirIndex (i, stagger_type); + + if (origin[i] < MY_GLOBAL_EP (GA->extras, myproc, stagger_index, i)) + { + mapping->global_startpoint[i] = origin[i]; + if (origin[i] < MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i)) + { + int npoints; + + npoints = (MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i) + - origin[i]) / mapping->downsample[i]; + if ((MY_GLOBAL_SP (GA->extras, myproc, stagger_index, i) + - origin[i]) % mapping->downsample[i]) + { + npoints++; + } + mapping->global_startpoint[i] += npoints * mapping->downsample[i]; + } + } + else + { + mapping->global_startpoint[i] = -1; + } + } + + /* compute the local start- and endpoint from the global ranges */ + mapping->totals = 1; + for (i = j = 0; i < vdim; i++) + { + stagger_index = CCTK_StaggerDirIndex (i, stagger_type); + + if (mapping->global_startpoint[i] >= 0 && + mapping->global_startpoint[i] < MY_GLOBAL_EP (GA->extras, myproc, + stagger_index, i)) + { + mapping->local_startpoint[i] = mapping->global_startpoint[i] - + GA->extras->lb[myproc][i]; + } + else + { + mapping->local_startpoint[i] = -1; + } + + if (mapping->global_endpoint[i] > MY_GLOBAL_SP (GA->extras, myproc, + stagger_index, i)) + { + mapping->local_endpoint[i] = MIN (MY_LOCAL_EP (GA->extras, stagger_index, i), + mapping->global_endpoint[i] - GA->extras->lb[myproc][i]); + } + else + { + mapping->local_endpoint[i] = -1; + } + +#ifdef DEBUG + printf ("direction %d: local ranges[%d, %d)\n", + i, mapping->local_startpoint[i], mapping->local_endpoint[i]); +#endif + + if (mapping->local_endpoint[i] < 0 || mapping->local_startpoint[i] < 0) + { + mapping->totals = 0; + mapping->local_endpoint[i] = mapping->local_startpoint[i]; + } + + if (mapping->do_dir[i]) + { + /* compute the local size in each hyperslab dimension */ + mapping->local_hsize[j] = (mapping->local_endpoint[i] - + mapping->local_startpoint[i]) / + mapping->downsample[i]; + if ((mapping->local_endpoint[i] - mapping->local_startpoint[i]) % + mapping->downsample[i]) + { + mapping->local_hsize[j]++; + } + mapping->totals *= mapping->local_hsize[j]; + j++; + } + } + } /* end of else branch for 'if (mapping->full_hyperslab)' */ + + /* diagonals in 3D are special: global_hsize is minimum of grid extents */ + if (mapping->is_diagonal_in_3D) + { + mapping->global_hsize[0] = GH->cctk_gsh[0]; + if (mapping->global_hsize[0] > GH->cctk_gsh[1]) + { + mapping->global_hsize[0] = GH->cctk_gsh[1]; + } + if (mapping->global_hsize[0] > GH->cctk_gsh[2]) + { + mapping->global_hsize[0] = GH->cctk_gsh[2]; + } + mapping->totals = mapping->global_hsize[0]; + } + +#ifdef DEBUG + printf ("total number of hyperslab data points: %d\n", mapping->totals); +#endif + + mapping->hdim = hdim; + mapping->vdim = vdim; + mapping->target_proc = target_proc; + mapping->conversion_fn = conversion_fn; + + /* add this mapping to the mapping list */ + if (mapping_list) + { + mapping_list->prev = mapping; + } + mapping->prev = NULL; + mapping->next = mapping_list; + mapping_list = mapping; + + mapping->handle = nmapping_list++; + + /* set the global hsize in the return arguments */ + if (hsize) + { + for (j = 0; j < hdim; j++) + { + hsize[j] = mapping->global_hsize[j]; + } + } + +#if 0 + if (mapping->totals > 0) + { + /* if requested, compute the offsets into the global hyperslab */ + if (hoffset_global) + { + for (i = j = 0; i < vdim; i++) + { + if (mapping->do_dir[i]) + { + if (mapping->full_hyperslab) + { + hoffset_global[j] = GA->extras->lb[myproc][i]; + } + else + { + hoffset_global[j] = (mapping->global_startpoint[i] - + origin[i]) / mapping->downsample[i]; + if (GA->connectivity->perme[i]) + { + hoffset_global[j] -= GA->extras->nghostzones[i]; + } + } +#ifdef DEBUG + printf ("hoffset_global, hsize in direction %d: %d, %d\n", + j, hoffset_global[j], mapping->local_hsize[j]); +#endif + j++; + } + } + } + } +#endif + + return (mapping->handle); +} + + +int Hyperslab_FreeMapping (CCTK_INT mapping_handle) +{ + hslab_mapping_t *mapping; + + + mapping = mapping_list; + while (mapping) + { + if (mapping->handle == mapping_handle) + { + if (mapping == mapping_list) + { + mapping_list = mapping_list->next; + } + else + { + if (mapping->next) + { + mapping->next->prev = mapping->prev; + } + if (mapping->prev) + { + mapping->prev->next = mapping->next; + } + } + free (mapping->vectors); + free (mapping); + return (0); + } + mapping = mapping->next; + } + + return (-1); +} + + +hslab_mapping_t *PUGHSlabi_GetMapping (int mapping_handle) +{ + hslab_mapping_t *mapping; + + + mapping = mapping_list; + while (mapping && mapping->handle != mapping_handle) + { + mapping = mapping->next; + } + + return (mapping); +} + + +static int checkFullHyperslab (const pGA *GA, + const CCTK_INT *origin, + const CCTK_INT *extent, + const hslab_mapping_t *mapping) +{ + int i; + int is_full_hyperslab; + + + is_full_hyperslab = mapping->hdim == GA->extras->dim; + + if (is_full_hyperslab) + { + for (i = 0; i < GA->extras->dim; i++) + { + is_full_hyperslab &= (origin[i] == 0); + is_full_hyperslab &= (extent[i] <= 0); + is_full_hyperslab &= (mapping->downsample[i] <= 1); + is_full_hyperslab &= (GA->connectivity->perme[i] == 0); + } + } + + return (is_full_hyperslab); +} diff --git a/src/PUGHSlab.h b/src/PUGHSlab.h index c2ff990..718edc6 100644 --- a/src/PUGHSlab.h +++ b/src/PUGHSlab.h @@ -25,12 +25,12 @@ typedef void (*t_hslabConversionFn) (const void *src, CCTK_INT dst_stride); /* function prototypes */ -CCTK_INT Hyperslab_Get (const cGH *GH, - CCTK_INT mapping_handle, - const CCTK_INT vindex, - const CCTK_INT timelevel, - const CCTK_INT hdatatype, - void *hdata); +CCTK_INT Hyperslab_Get (const cGH *GH, + CCTK_INT mapping_handle, + CCTK_INT vindex, + CCTK_INT timelevel, + CCTK_INT hdatatype, + void *hdata); CCTK_INT Hyperslab_GetList (const cGH *GH, CCTK_INT mapping_handle, CCTK_INT num_arrays, diff --git a/src/PUGHSlabi.h b/src/PUGHSlabi.h index f251544..47f8907 100644 --- a/src/PUGHSlabi.h +++ b/src/PUGHSlabi.h @@ -26,15 +26,21 @@ typedef struct hslab_mapping_t int target_proc; int hdim; int vdim; - int *global_origin; /* vdim */ - int *directions; /* hdim * vdim */ - int *extents; /* hdim */ - int *downsample; /* hdim */ - int *global_hsize; /* hdim */ + int *vectors; + int *local_startpoint; /* vdim */ + int *local_endpoint; /* vdim */ + int *global_startpoint; /* vdim */ + int *global_endpoint; /* vdim */ int *do_dir; /* vdim */ + int *downsample; /* vdim */ + int *global_hsize; /* hdim */ + int *local_hsize; /* hdim */ + int *directions; /* hdim * vdim */ int is_diagonal_in_3D; + unsigned int totals; t_hslabConversionFn conversion_fn; struct hslab_mapping_t *prev, *next; + int full_hyperslab; } hslab_mapping_t; diff --git a/src/make.code.defn b/src/make.code.defn index 2463bb7..7d322e8 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,4 +2,4 @@ # $Header$ # Source files in this directory -SRCS = Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c +SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c -- cgit v1.2.3