From 667cbe19ba0fd827de899e77113520f6751b64ef Mon Sep 17 00:00:00 2001 From: tradke Date: Thu, 12 Dec 2002 13:56:26 +0000 Subject: Preliminary support for both the old and the new global interpolator API. git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@33 1c20744c-e24a-42ec-9533-f5004cb800e5 --- src/Operator.c | 571 ++++++++++++++++++++++++----------------------------- src/Startup.c | 18 +- src/pughInterpGH.h | 39 ++-- 3 files changed, 295 insertions(+), 333 deletions(-) diff --git a/src/Operator.c b/src/Operator.c index fc6e2cf..8064fd3 100644 --- a/src/Operator.c +++ b/src/Operator.c @@ -2,9 +2,9 @@ @file Operator.c @date Tue Apr 15 18:22:45 1997 @author Paul Walker - @desc + @desc Definition of interpolation operators for regular uniform grids. - @enddesc + @enddesc @history @date Sun Jul 04 1999 @@ -18,11 +18,9 @@ @@*/ #include -#include #include #include "cctk.h" -#include "cctk_Parameters.h" #include "CactusPUGH/PUGH/src/include/pugh.h" #include "pughInterpGH.h" @@ -35,52 +33,50 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Operator_c) /* macro do sort interpolation results from a single communication buffer into their appropriate output arrays */ -#define SORT_TYPED_ARRAY(cctk_type, source, dest, points, index_array, offset)\ - { \ - int i; \ - cctk_type *src, *dst; \ +#define SORT_TYPED_ARRAY(cctk_type) \ + { \ + int _i; \ + cctk_type *_src, *_dst; \ \ \ - src = (cctk_type *) source; \ - dst = (cctk_type *) dest; \ - for (i = 0; i < points; i++) \ - { \ - dst[index_array[i + offset]] = *src++; \ - } \ - source = (char *) src; \ - } - + _src = (cctk_type *) this->buf; \ + _dst = (cctk_type *) output_arrays[array]; \ + for (_i = 0; _i < myGH->N_points_from[proc]; _i++) \ + { \ + _dst[myGH->indices[_i + offset]] = *_src++; \ + } \ + this->buf = (char *) _src; \ + } #ifdef CCTK_MPI /* internal structure describing a handle for a single CCTK data type */ typedef struct { - int vtypesize; /* variable type's size in bytes */ - MPI_Datatype mpi_type; /* corresponding MPI datatype */ - int num_arrays; /* number of in/out arrays */ - int array; /* iterator for num_arrays */ - void *send_buffer; /* communication send buffer for this type */ - void *recv_buffer; /* communication receive buffer for this type */ - char *buffer; /* work pointer for send_buffer */ + int vtypesize; /* variable type's size in bytes */ + MPI_Datatype mpitype; /* corresponding MPI datatype */ + int N_arrays; /* number of in/out arrays */ + void *sendbuf; /* communication send buffer for this type */ + void *recvbuf; /* communication receive buffer for this type */ + char *buf; /* work pointer for sendbuf */ } type_desc_t; #endif /* prototypes of routines defined in this source file */ -static int PUGHInterp_CheckArguments (cGH *GH, - int num_dims, - int num_points, - int num_in_arrays, - int num_out_arrays, - const int interp_coord_array_types[]); +static int CheckArguments (const cGH *GH, + int N_dims, + int N_points, + int N_input_arrays, + int N_output_arrays, + const int interp_coord_array_types[]); #ifdef CCTK_MPI -static int GetLocalCoords (cGH *GH, - int num_points, - const char *coord_system, +static int GetLocalCoords (const cGH *GH, + int N_points, + const char *coord_system_name, const pGExtras *extras, const CCTK_REAL *coords[], - int *num_local_points, + int *N_local_points, CCTK_REAL **local_coords); #endif @@ -111,55 +107,55 @@ static int GetLocalCoords (cGH *GH, @vtype int @vio in @endvar - @var coord_system + @var coord_system_name @vdesc name of coordinate system to use for interpolation @vtype const char * @vio in @endvar - @var num_points + @var N_points @vdesc number of points to be interpolated on this processor @vtype int @vio in @endvar - @var num_in_array_indices + @var N_input_arrays @vdesc number of input arrays (given by their indices) to interpolate from @vtype int @vio in @endvar - @var num_out_arrays + @var N_output_arrays @vdesc number of output arrays to interpolate to @vtype int @vio in @endvar @var interp_coord_arrays @vdesc coordinates of points to interpolate at - @vtype void *const [num_dims] + @vtype void *const [N_dims] @vio in @endvar @var interp_coord_array_types @vdesc CCTK data type of coordinate arrays - @vtype int [num_dims] + @vtype int [N_dims] @vio in @endvar - @var in_arrays + @var input_arrays @vdesc list of input arrays to interpolate on - @vtype void *[num_in_arrays] + @vtype void *[N_input_arrays] @vio in @endvar - @var in_array_types + @var input_array_types @vdesc CCTK data types of input arrays - @vtype int [num_in_arrays] + @vtype int [N_input_arrays] @vio in @endvar - @var out_arrays + @var output_arrays @vdesc list of output arrays to interpolate to - @vtype void *const [num_out_arrays] + @vtype void *const [N_output_arrays] @vio out @endvar - @var out_array_types + @var output_array_types @vdesc CCTK data types of output arrays - @vtype int [num_out_arrays] + @vtype int [N_output_arrays] @vio in @endvar @@ -171,55 +167,53 @@ static int GetLocalCoords (cGH *GH, @@*/ int PUGHInterp_InterpGV (cGH *GH, int order, - const char *coord_system, - int num_points, - int num_in_array_indices, - int num_out_arrays, + const char *coord_system_name, + int N_points, + int N_input_arrays, + int N_output_arrays, const void *const interp_coord_arrays[], const int interp_coord_array_types[], - const int in_array_indices[], - void *const out_arrays[], - const int out_array_types[]) + const int input_array_indices[], + void *const output_arrays[], + const int output_array_types[]) { const CCTK_REAL **data; CCTK_REAL *interp_local_coords; CCTK_REAL *origin, *delta; int nprocs; - int timelevel; - int dim, num_dims, *dims; + int dim, N_dims, *dims; int array; int point; int retval; - const void **in_arrays; + const void **input_arrays; pGH *pughGH; pGExtras *extras; cGroupDynamicData group_data; #ifdef CCTK_MPI - int num_local_points; + int N_local_points; int offset; int proc; int type; - void **local_out_arrays; + void **local_output_arrays; pughInterpGH *myGH; - int max_type; - type_desc_t *type_desc; + int maxtype; + type_desc_t *this, *type_desc; #endif /* get dimensionality of the coordinate system */ - num_dims = CCTK_CoordSystemDim (coord_system); - if (num_dims <= 0) + N_dims = CCTK_CoordSystemDim (coord_system_name); + if (N_dims <= 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Cannot get dimensions of coordinate system '%s'", - coord_system); + coord_system_name); return (-1); } /* check other arguments */ - retval = PUGHInterp_CheckArguments (GH, num_dims, num_points, - num_in_array_indices, num_out_arrays, - interp_coord_array_types); + retval = CheckArguments (GH, N_dims, N_points, N_input_arrays, + N_output_arrays, interp_coord_array_types); if (retval <= 0) { return (retval); @@ -229,46 +223,46 @@ int PUGHInterp_InterpGV (cGH *GH, /* NOTE: getting the dimensions should be a flesh routine as well for now we get the dimensions of every coordinate and take the i'th element - this is inconsistent !! */ - dims = (int *) malloc (2 * num_dims * sizeof (int)); - origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL)); - delta = origin + num_dims; - in_arrays = malloc (num_in_array_indices * sizeof (void *)); - for (dim = 0; dim < num_dims; dim++) + dims = malloc (2 * N_dims * sizeof (int)); + origin = malloc (2 * N_dims * sizeof (CCTK_REAL)); + delta = origin + N_dims; + input_arrays = malloc (N_input_arrays * sizeof (void *)); + for (dim = 0; dim < N_dims; dim++) { - CCTK_GrouplshVI (GH, num_dims, dims + num_dims, - CCTK_CoordIndex (dim + 1, NULL, coord_system)); - dims[dim] = dims[dim + num_dims]; + CCTK_GrouplshVI (GH, N_dims, dims + N_dims, + CCTK_CoordIndex (dim + 1, NULL, coord_system_name)); + dims[dim] = dims[dim + N_dims]; CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL, - coord_system); + coord_system_name); delta[dim] = (delta[dim] - origin[dim]) / dims[dim]; } /* get extension handle for PUGH */ - pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH"); + pughGH = CCTK_GHExtension (GH, "PUGH"); extras = NULL; /* check that the input arrays dimensions match the coordinate system (but their dimensionality can be less) */ - for (array = 0; array < num_in_array_indices; array++) + for (array = 0; array < N_input_arrays; array++) { if (CCTK_GroupDynamicData (GH, - CCTK_GroupIndexFromVarI(in_array_indices[array]), + CCTK_GroupIndexFromVarI(input_array_indices[array]), &group_data) < 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Invalid input array index %d", - in_array_indices[array]); + input_array_indices[array]); retval = -1; continue; } - if (group_data.dim > num_dims) + if (group_data.dim > N_dims) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Input array variable with index %d has more dimensions " "than coordinate system '%s'", - in_array_indices[array], coord_system); + input_array_indices[array], coord_system_name); retval = -1; continue; } @@ -278,7 +272,7 @@ int PUGHInterp_InterpGV (cGH *GH, CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Dimensions of input array variable with index %d " "doesn't match with coordinate system '%s'", - in_array_indices[array], coord_system); + input_array_indices[array], coord_system_name); retval = -1; } @@ -286,50 +280,48 @@ int PUGHInterp_InterpGV (cGH *GH, highest dimension (need this information later on) */ if (! extras || extras->dim < group_data.dim) { - extras = ((pGA *) pughGH->variables[in_array_indices[array]][0])->extras; + extras = ((pGA *) pughGH->variables[input_array_indices[array]][0])->extras; } /* get the data pointer to the input array (use current timelevel) */ - timelevel = 0; - in_arrays[array] = CCTK_VarDataPtrI(GH, timelevel, in_array_indices[array]); + input_arrays[array] = CCTK_VarDataPtrI (GH, 0, input_array_indices[array]); } if (retval < 0) { - free (in_arrays); + free (input_arrays); free (origin); free (dims); return (retval); } /* single-processor case is easy: no communication or buffering, just direct - interpolation of interp_coord_arrays from in_arrays into out_arrays */ + interpolation of interp_coord_arrays from input_arrays into output_arrays */ nprocs = CCTK_nProcs (GH); if (nprocs == 1) { /* sort the individual interpolation coordinate arrays into a single one */ - interp_local_coords = (CCTK_REAL *) malloc (num_dims * num_points * - sizeof (CCTK_REAL)); + interp_local_coords = malloc (N_dims * N_points * sizeof (CCTK_REAL)); data = (const CCTK_REAL **) interp_coord_arrays; - for (point = 0; point < num_points; point++) + for (point = 0; point < N_points; point++) { - for (dim = 0; dim < num_dims; dim++) + for (dim = 0; dim < N_dims; dim++) { *interp_local_coords++ = data[dim][point]; } } - interp_local_coords -= num_dims * num_points; + interp_local_coords -= N_dims * N_points; /* call the interpolator function */ retval = PUGHInterp_Interpolate (order, - num_points, num_dims, num_out_arrays, + N_points, N_dims, N_output_arrays, dims, interp_local_coords, origin, delta, - out_array_types, in_arrays, - out_array_types, out_arrays); + output_array_types, input_arrays, + output_array_types, output_arrays); /* free allocated resources */ free (interp_local_coords); - free (in_arrays); + free (input_arrays); free (origin); free (dims); @@ -358,52 +350,50 @@ int PUGHInterp_InterpGV (cGH *GH, /* get the maximum value of the output array CCTK data types NOTE: we assume that CCTK data types are defined as consecutive positive constants starting from zero */ - max_type = 0; - for (array = 0; array < num_out_arrays; array++) + for (array = maxtype = 0; array < N_output_arrays; array++) { - if (max_type < out_array_types[array]) + if (maxtype < output_array_types[array]) { - max_type = out_array_types[array]; + maxtype = output_array_types[array]; } } /* now allocate an array of structures for all potential types */ - type_desc = (type_desc_t *) calloc (max_type + 1, sizeof (type_desc_t)); + type_desc = calloc (maxtype + 1, sizeof (type_desc_t)); /* count the number of arrays of same type - (the num_arrays element was already initialized to zero by calloc() */ - for (array = 0; array < num_out_arrays; array++) + (the N_arrays element was already initialized to zero by calloc() */ + for (array = 0; array < N_output_arrays; array++) { - type_desc[out_array_types[array]].num_arrays++; + type_desc[output_array_types[array]].N_arrays++; } /* fill in the type description information */ - retval = 0; - for (type = 0; type <= max_type; type++) + for (type = retval = 0, this = type_desc; type <= maxtype; type++, this++) { - if (type_desc[type].num_arrays > 0) + if (this->N_arrays > 0) { /* get the variable type size in bytes */ - type_desc[type].vtypesize = CCTK_VarTypeSize (type); - if (type_desc[type].vtypesize <= 0) + this->vtypesize = CCTK_VarTypeSize (type); + if (this->vtypesize <= 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Invalid variable type %d passed, " "arrays of such type will be skipped during interpolation", type); - type_desc[type].num_arrays = 0; + this->N_arrays = 0; continue; } /* get the MPI data type to use for communicating such a CCTK data type */ - type_desc[type].mpi_type = PUGH_MPIDataType (pughGH, type); - if (type_desc[type].mpi_type == MPI_DATATYPE_NULL) + this->mpitype = PUGH_MPIDataType (pughGH, type); + if (this->mpitype == MPI_DATATYPE_NULL) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "No MPI data type defined for variable type %d, " "arrays of such type will be skipped during interpolation", type); - type_desc[type].num_arrays = 0; + this->N_arrays = 0; continue; } @@ -414,7 +404,7 @@ int PUGHInterp_InterpGV (cGH *GH, /* check that there's at least one array with a valid CCTK data type */ if (retval <= 0) { - free (in_arrays); + free (input_arrays); free (dims); free (origin); free (type_desc); @@ -424,14 +414,14 @@ int PUGHInterp_InterpGV (cGH *GH, /* map the requested points to interpolate at onto the processors they belong to and gather the coordinates of all the points that this processor owns - the number of processor-local points is returned in num_local_points, + the number of processor-local points is returned in N_local_points, their coordinates in interp_local_coords */ - retval = GetLocalCoords (GH, num_points, coord_system, extras, + retval = GetLocalCoords (GH, N_points, coord_system_name, extras, (const CCTK_REAL **) interp_coord_arrays, - &num_local_points, &interp_local_coords); + &N_local_points, &interp_local_coords); if (retval) { - free (in_arrays); + free (input_arrays); free (origin); free (dims); free (type_desc); @@ -444,45 +434,41 @@ int PUGHInterp_InterpGV (cGH *GH, If there are no points to process on this processor set the buffer pointer to an invalid but non-NULL value otherwise we might get trouble with NULL pointers in MPI_Alltoallv () */ - for (type = 0; type <= max_type; type++) + for (type = 0, this = type_desc; type <= maxtype; type++, this++) { - if (type_desc[type].num_arrays > 0 && num_local_points > 0) + if (this->N_arrays > 0 && N_local_points > 0) { - type_desc[type].send_buffer = malloc (num_local_points * - type_desc[type].num_arrays * - type_desc[type].vtypesize); - type_desc[type].buffer = (char *) type_desc[type].send_buffer; - type_desc[type].array = 0; + this->sendbuf = malloc (N_local_points * this->N_arrays *this->vtypesize); + this->buf = this->sendbuf; } else { /* dereferencing such an address should code crash on most systems */ - type_desc[type].send_buffer = (void *) type_desc[type].vtypesize; + this->sendbuf = (void *) this->vtypesize; } } /* get extension handle for interp */ - myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp"); + myGH = CCTK_GHExtension (GH, "PUGHInterp"); - /* allocate new out_arrays array for local interpolation results + /* allocate new output_arrays array for local interpolation results from this processor */ - local_out_arrays = calloc (num_out_arrays, sizeof (void *)); + local_output_arrays = calloc (N_output_arrays, sizeof (void *)); /* now, in a loop over all processors, do the interpolation and put the results in the communication buffer at the proper offset */ for (proc = 0; proc < nprocs; proc++) { - for (type = 0; type <= max_type; type++) + for (type = 0, this = type_desc; type <= maxtype; type++, this++) { - if (type_desc[type].num_arrays > 0) + if (this->N_arrays > 0) { - for (array = 0; array < num_out_arrays; array++) + for (array = 0; array < N_output_arrays; array++) { - if (out_array_types[array] == type) + if (output_array_types[array] == type) { - local_out_arrays[array] = type_desc[type].buffer; - type_desc[type].buffer += myGH->num_points_to[proc] * - type_desc[type].vtypesize; + local_output_arrays[array] = this->buf; + this->buf += myGH->N_points_to[proc] * this->vtypesize; } } } @@ -491,24 +477,24 @@ int PUGHInterp_InterpGV (cGH *GH, /* call the interpolation operator to process all points of all output arrays for this processor */ PUGHInterp_Interpolate (order, - myGH->num_points_to[proc], num_dims, num_out_arrays, + myGH->N_points_to[proc], N_dims, N_output_arrays, dims, interp_local_coords, origin, delta, - out_array_types, in_arrays, - out_array_types, local_out_arrays); + output_array_types, input_arrays, + output_array_types, local_output_arrays); /* have to add offset for this processor to coordinates array */ - interp_local_coords += myGH->num_points_to[proc] * num_dims; + interp_local_coords += myGH->N_points_to[proc] * N_dims; } /* end of loop over all processors */ /* don't need these anymore */ - if (num_local_points > 0) + if (N_local_points > 0) { - interp_local_coords -= num_local_points * num_dims; + interp_local_coords -= N_local_points * N_dims; free (interp_local_coords); } - free (local_out_arrays); - free (in_arrays); + free (local_output_arrays); + free (input_arrays); free (origin); free (dims); @@ -517,143 +503,114 @@ int PUGHInterp_InterpGV (cGH *GH, by other processors Since all the locally computed results are in a single contiguous buffer we need to call MPI_Alltoall() only once for each CCTK data type. */ - for (type = 0; type <= max_type; type++) + for (type = 0, this = type_desc; type <= maxtype; type++, this++) { /* skip unused types */ - if (type_desc[type].num_arrays <= 0) + if (this->N_arrays <= 0) { continue; } /* set up the communication (this is type-independent) */ - myGH->send_count[0] = type_desc[type].num_arrays * myGH->num_points_to[0]; - myGH->recv_count[0] = type_desc[type].num_arrays * myGH->num_points_from[0]; - myGH->send_displ[0] = myGH->recv_displ[0] = 0; + myGH->sendcnt[0] = this->N_arrays * myGH->N_points_to[0]; + myGH->recvcnt[0] = this->N_arrays * myGH->N_points_from[0]; + myGH->senddispl[0] = myGH->recvdispl[0] = 0; for (proc = 1; proc < pughGH->nprocs; proc++) { - myGH->send_count[proc] = type_desc[type].num_arrays * - myGH->num_points_to[proc]; - myGH->recv_count[proc] = type_desc[type].num_arrays * - myGH->num_points_from[proc]; - myGH->send_displ[proc] = myGH->send_displ[proc-1] + - myGH->send_count[proc-1]; - myGH->recv_displ[proc] = myGH->recv_displ[proc-1] + - myGH->recv_count[proc-1]; + myGH->sendcnt[proc] = this->N_arrays * myGH->N_points_to[proc]; + myGH->recvcnt[proc] = this->N_arrays * myGH->N_points_from[proc]; + myGH->senddispl[proc] = myGH->senddispl[proc-1] + myGH->sendcnt[proc-1]; + myGH->recvdispl[proc] = myGH->recvdispl[proc-1] + myGH->recvcnt[proc-1]; } /* allocate buffer for receiving my own requested points */ /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */ - if (num_points > 0) + if (N_points > 0) { - type_desc[type].recv_buffer = malloc (num_points * - type_desc[type].num_arrays * - type_desc[type].vtypesize); + this->recvbuf = malloc (N_points * this->N_arrays * this->vtypesize); } else { /* access to such a fake address should crash the code on most systems */ - type_desc[type].recv_buffer = (void *) type_desc[type].vtypesize; + this->recvbuf = (void *) this->vtypesize; } /* now exchange the data for this CCTK data type */ - CACTUS_MPI_ERROR (MPI_Alltoallv (type_desc[type].send_buffer, - myGH->send_count, myGH->send_displ, - type_desc[type].mpi_type, - type_desc[type].recv_buffer, - myGH->recv_count, myGH->recv_displ, - type_desc[type].mpi_type, + CACTUS_MPI_ERROR (MPI_Alltoallv (this->sendbuf, myGH->sendcnt, + myGH->senddispl, this->mpitype, + this->recvbuf, myGH->recvcnt, + myGH->recvdispl, this->mpitype, pughGH->PUGH_COMM_WORLD)); /* now that the data is sent we don't need the buffer anymore */ - if (num_local_points > 0) + if (N_local_points > 0) { - free (type_desc[type].send_buffer); + free (this->sendbuf); } /* no sort neccessary if there are no points */ - if (num_points <= 0) + if (N_points <= 0) { continue; } /* Go back from processor-sorted data to input-ordered data. - The creation of the indices array above makes this not so bad. */ - offset = 0; - type_desc[type].buffer = (char *) type_desc[type].recv_buffer; - for (proc = 0; proc < nprocs; proc++) + The creation of the indices array above makes this not so bad. */ + this->buf = this->recvbuf; + for (proc = offset = 0; proc < nprocs; proc++) { - for (array = 0; array < num_out_arrays; array++) + for (array = 0; array < N_output_arrays; array++) { - if (out_array_types[array] != type) + if (output_array_types[array] != type) { continue; } /* now do the sorting according to the CCTK data type */ - if (out_array_types[array] == CCTK_VARIABLE_CHAR) + if (output_array_types[array] == CCTK_VARIABLE_CHAR) { - SORT_TYPED_ARRAY (CCTK_BYTE, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_BYTE); } - else if (out_array_types[array] == CCTK_VARIABLE_INT) + else if (output_array_types[array] == CCTK_VARIABLE_INT) { - SORT_TYPED_ARRAY (CCTK_INT, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_INT); } - else if (out_array_types[array] == CCTK_VARIABLE_REAL) + else if (output_array_types[array] == CCTK_VARIABLE_REAL) { - SORT_TYPED_ARRAY (CCTK_REAL, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_REAL); } - else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX) + else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX) { - SORT_TYPED_ARRAY (CCTK_COMPLEX, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_COMPLEX); } #ifdef CCTK_REAL4 - else if (out_array_types[array] == CCTK_VARIABLE_REAL4) + else if (output_array_types[array] == CCTK_VARIABLE_REAL4) { - SORT_TYPED_ARRAY (CCTK_REAL4, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_REAL4); } - else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX8) + else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX8) { - SORT_TYPED_ARRAY (CCTK_COMPLEX8, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_COMPLEX8); } #endif #ifdef CCTK_REAL8 - else if (out_array_types[array] == CCTK_VARIABLE_REAL8) + else if (output_array_types[array] == CCTK_VARIABLE_REAL8) { - SORT_TYPED_ARRAY (CCTK_REAL8, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_REAL8); } - else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX16) + else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX16) { - SORT_TYPED_ARRAY (CCTK_COMPLEX16, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_COMPLEX16); } #endif #ifdef CCTK_REAL16 - else if (out_array_types[array] == CCTK_VARIABLE_REAL16) + else if (output_array_types[array] == CCTK_VARIABLE_REAL16) { - SORT_TYPED_ARRAY (CCTK_REAL16, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_REAL16); } - else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX32) + else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX32) { - SORT_TYPED_ARRAY (CCTK_COMPLEX32, type_desc[type].buffer, - out_arrays[array], myGH->num_points_from[proc], - myGH->indices, offset); + SORT_TYPED_ARRAY (CCTK_COMPLEX32); } #endif else @@ -664,12 +621,12 @@ int PUGHInterp_InterpGV (cGH *GH, } /* end of loop over all output arrays */ /* advance the offset into the communication receive buffer */ - offset += myGH->num_points_from[proc]; + offset += myGH->N_points_from[proc]; } /* end of loop over all processors */ /* this communication receive buffer isn't needed anymore */ - free (type_desc[type].recv_buffer); + free (this->recvbuf); } /* end of loop over all types */ @@ -696,33 +653,33 @@ int PUGHInterp_InterpGV (cGH *GH, @routine GetLocalCoords @date Sun Jul 04 1999 @author Thomas Radke - @desc + @desc Collect the coordinates of all points to be processed locally - into an array coords[num_dims][num_local_points]. + into an array coords[N_dims][N_local_points]. This means for the single-processor case to sort - inCoords1[num_points], inCoords2[npoints], ..., inCoords[npoints] - into coords[num_dims][num_points] + inCoords1[N_points], inCoords2[npoints], ..., inCoords[npoints] + into coords[N_dims][N_points] - where num_points == num_local_points. + where N_points == N_local_points. In the multiprocessor case all processors map their points' coordinates to the processor that owns this point and exchange this information via MPI_Alltoall (). - num_local_points is then the number of all processors' points to be + N_local_points is then the number of all processors' points to be interpolated locally on this processor. This routine returns the number of points to be processed locally and a pointer to the allocated array of their coordinates. - @enddesc + @enddesc @var GH @vdesc Pointer to CCTK grid hierarchy - @vtype cGH * + @vtype const cGH * @vio in @endvar - @var num_points + @var N_points @vdesc number of points to be interpolated on this processor @vtype int @vio in @@ -733,37 +690,36 @@ int PUGHInterp_InterpGV (cGH *GH, @vtype int @vio in @endvar - @var num_dims + @var N_dims @vdesc number of coordinate dimensions for each point @vtype int @vio in @endvar @var coords @vdesc coordinates of each point to be interpolated on this processor - @vtype CCTK_REAL array of size num_dims + @vtype CCTK_REAL array of size N_dims @vio in @endvar - @var num_local_points + @var N_local_points @vdesc number of points to be processed by this processor @vtype int * @vio out @endvar @var local_coords @vdesc coordinates of each point to be processed by this processor - @vtype pointer to CCTK_REAL array of dims[num_dims][num_local_points] + @vtype pointer to CCTK_REAL array of dims[N_dims][N_local_points] @vio out @endvar @@*/ -static int GetLocalCoords (cGH *GH, - int num_points, - const char *coord_system, +static int GetLocalCoords (const cGH *GH, + int N_points, + const char *coord_system_name, const pGExtras *extras, const CCTK_REAL *coords[], - int *num_local_points, + int *N_local_points, CCTK_REAL **local_coords) { - DECLARE_CCTK_PARAMETERS int dim, point; int tmp, proc, nprocs; pGH *pughGH; @@ -773,14 +729,15 @@ static int GetLocalCoords (cGH *GH, CCTK_REAL *proc_coords; #define FUDGE 0.0 + /* get GH extension handles for PUGHInterp and PUGH */ - myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp"); - pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH"); + myGH = CCTK_GHExtension (GH, "PUGHInterp"); + pughGH = CCTK_GHExtension (GH, "PUGH"); - /* This holds the proccessor for *each* of num_points points */ - if (num_points > 0) + /* This holds the proccessor for *each* of N_points points */ + if (N_points > 0) { - myGH->whichproc = (int *) malloc (2 * num_points * sizeof (int)); + myGH->whichproc = malloc (2 * N_points * sizeof (int)); } else { @@ -788,20 +745,20 @@ static int GetLocalCoords (cGH *GH, } /* indices[] is used to make the sorting easier when receiving the output data */ - myGH->indices = myGH->whichproc + num_points; + myGH->indices = myGH->whichproc + N_points; /* initialize whichproc with invalid processor number -1 */ - for (point = 0; point < num_points; point++) + for (point = 0; point < N_points; point++) { myGH->whichproc[point] = -1; } - - /* initialize num_points_from to 0 for counting it up in the following loop */ + + /* initialize N_points_from to 0 for counting it up in the following loop */ nprocs = CCTK_nProcs (GH); - memset (myGH->num_points_from, 0, nprocs * sizeof (CCTK_INT)); + memset (myGH->N_points_from, 0, nprocs * sizeof (CCTK_INT)); /* allocate the ranges for my local coordinates */ - range_min = (CCTK_REAL *) malloc (4 * extras->dim * sizeof (CCTK_REAL)); + range_min = malloc (4 * extras->dim * sizeof (CCTK_REAL)); range_max = range_min + 1*extras->dim; origin = range_min + 2*extras->dim; delta = range_min + 3*extras->dim; @@ -809,7 +766,7 @@ static int GetLocalCoords (cGH *GH, /* get the global origin and delta of the coordinate system */ for (dim = 0; dim < extras->dim; dim++) { - CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system); + CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system_name); delta[dim] = (delta[dim] - origin[dim]) / (extras->nsize[dim]-1); } @@ -820,7 +777,7 @@ static int GetLocalCoords (cGH *GH, { /* compute the coordinate ranges */ /* TODO: use bbox instead -- but the bboxes of other processors - are now known */ + are now known */ int const has_lower = extras->lb[proc][dim] == 0; int const has_upper = extras->ub[proc][dim] == GH->cctk_gsh[dim]-1; range_min[dim] = origin[dim] + (extras->lb[proc][dim] + (!has_lower) * (extras->nghostzones[dim]-0.5) - FUDGE)*delta[dim]; @@ -828,7 +785,7 @@ static int GetLocalCoords (cGH *GH, } /* and now which point will be processed by what processor */ - for (point = 0; point < num_points; point++) + for (point = 0; point < N_points; point++) { /* skip points which have already been located */ if (myGH->whichproc[point] >= 0) @@ -850,7 +807,7 @@ static int GetLocalCoords (cGH *GH, if (tmp == extras->dim) { myGH->whichproc[point] = proc; - myGH->num_points_from[proc]++; + myGH->N_points_from[proc]++; } } } @@ -858,13 +815,12 @@ static int GetLocalCoords (cGH *GH, free (range_min); /* make sure that all points could be mapped onto a processor */ - tmp = 0; - for (point = 0; point < num_points; point++) + for (point = tmp = 0; point < N_points; point++) { if (myGH->whichproc[point] < 0) { int i; - char *msg = (char *) malloc (80 + extras->dim*20); + char *msg = malloc (80 + extras->dim*20); sprintf (msg, "Unable to locate point %d [%f", @@ -889,7 +845,7 @@ static int GetLocalCoords (cGH *GH, return (-1); } - /* Now we want to resolve the num_points_from[]. Currently this is + /* Now we want to resolve the N_points_from[]. Currently this is the form of ( in 2 proc mode ) P1: Num from P1 NFP2 P2: NFP1 NFP2 @@ -901,16 +857,16 @@ static int GetLocalCoords (cGH *GH, Since NTP1 = NFP2 (and this works in more proc mode too) this is an all-to-all communication. */ - CACTUS_MPI_ERROR (MPI_Alltoall (myGH->num_points_from, 1, PUGH_MPI_INT, - myGH->num_points_to, 1, PUGH_MPI_INT, + CACTUS_MPI_ERROR (MPI_Alltoall (myGH->N_points_from, 1, PUGH_MPI_INT, + myGH->N_points_to, 1, PUGH_MPI_INT, pughGH->PUGH_COMM_WORLD)); #ifdef PUGHINTERP_DEBUG for (proc = 0; proc < nprocs; proc++) { printf ("processor %d <-> %d From: %d To: %d\n", - CCTK_MyProc (GH), proc, myGH->num_points_from[proc], - myGH->num_points_to[proc]); + CCTK_MyProc (GH), proc, myGH->N_points_from[proc], + myGH->N_points_to[proc]); } #endif @@ -918,37 +874,34 @@ static int GetLocalCoords (cGH *GH, and how many to send to each processor. So first we have to send the locations to the processors which hold our data. This means I send coords[dim][point] to whichproc[point]. - I have num_points_from[proc] to send to each processor. + I have N_points_from[proc] to send to each processor. */ - + /* This is backwards in the broadcast location; the number of points we are getting is how many everyone else is sending to us, - eg, num_points_to, not how many we get back from everyone else, - eg, num_points_from. The number we are sending, of course, is - all of our locations, eg, num_points */ - *num_local_points = 0; + eg, N_points_to, not how many we get back from everyone else, + eg, N_points_from. The number we are sending, of course, is + all of our locations, eg, N_points */ + *N_local_points = 0; for (proc = 0; proc < nprocs; proc++) { - *num_local_points += myGH->num_points_to[proc]; + *N_local_points += myGH->N_points_to[proc]; } #ifdef PUGHINTERP_DEBUG printf ("processor %d gets %d points in total\n", - CCTK_MyProc (GH), *num_local_points); + CCTK_MyProc (GH), *N_local_points); #endif /* allocate the local coordinates array (sorted in processor order) and the resulting coordinates array that I have to process */ - proc_coords = (CCTK_REAL *) malloc (extras->dim * num_points * - sizeof (CCTK_REAL)); - *local_coords = (CCTK_REAL *) malloc (extras->dim * *num_local_points * - sizeof (CCTK_REAL)); + proc_coords = malloc (extras->dim * N_points * sizeof (CCTK_REAL)); + *local_coords = malloc (extras->dim * *N_local_points * sizeof (CCTK_REAL)); /* now sort my own coordinates as tupels of [extras->dim] */ - tmp = 0; - for (proc = 0; proc < nprocs; proc++) + for (proc = tmp = 0; proc < nprocs; proc++) { - for (point = 0; point < num_points; point++) + for (point = 0; point < N_points; point++) { if (myGH->whichproc[point] == proc) { @@ -964,25 +917,23 @@ static int GetLocalCoords (cGH *GH, /* So load up the send and recv stuff */ /* Send extras->dim elements per data point */ - myGH->send_count[0] = extras->dim * myGH->num_points_from[0]; - myGH->recv_count[0] = extras->dim * myGH->num_points_to[0]; - myGH->send_displ[0] = myGH->recv_displ[0] = 0; + myGH->sendcnt[0] = extras->dim * myGH->N_points_from[0]; + myGH->recvcnt[0] = extras->dim * myGH->N_points_to[0]; + myGH->senddispl[0] = myGH->recvdispl[0] = 0; for (proc = 1; proc < nprocs; proc++) { - myGH->send_count[proc] = extras->dim * myGH->num_points_from[proc]; - myGH->recv_count[proc] = extras->dim * myGH->num_points_to[proc]; - myGH->send_displ[proc] = myGH->send_displ[proc-1] + - myGH->send_count[proc-1]; - myGH->recv_displ[proc] = myGH->recv_displ[proc-1] + - myGH->recv_count[proc-1]; + myGH->sendcnt[proc] = extras->dim * myGH->N_points_from[proc]; + myGH->recvcnt[proc] = extras->dim * myGH->N_points_to[proc]; + myGH->senddispl[proc] = myGH->senddispl[proc-1] + myGH->sendcnt[proc-1]; + myGH->recvdispl[proc] = myGH->recvdispl[proc-1] + myGH->recvcnt[proc-1]; } - + /* Great, and now exchange the coordinates and collect the ones that I have to process in *local_coords[] */ - CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->send_count, - myGH->send_displ, PUGH_MPI_REAL, - *local_coords, myGH->recv_count, - myGH->recv_displ, PUGH_MPI_REAL, + CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->sendcnt, + myGH->senddispl, PUGH_MPI_REAL, + *local_coords, myGH->recvcnt, + myGH->recvdispl, PUGH_MPI_REAL, pughGH->PUGH_COMM_WORLD)); /* don't need this anymore */ @@ -994,7 +945,7 @@ static int GetLocalCoords (cGH *GH, /*@@ - @routine PUGHInterp_CheckArguments + @routine CheckArguments @date Thu 25 Jan 2001 @author Thomas Radke @desc @@ -1007,32 +958,32 @@ static int GetLocalCoords (cGH *GH, @var GH @vdesc Pointer to CCTK grid hierarchy - @vtype cGH * + @vtype const cGH * @vio in @endvar - @var num_dims + @var N_dims @vdesc dimensionality of the underlying grid @vtype int @vio in @endvar - @var num_points + @var N_points @vdesc number of points to interpolate at @vtype int @vio in @endvar - @var num_in_arrays + @var N_input_arrays @vdesc number of passed input arrays @vtype int @vio in @endvar - @var num_out_arrays + @var N_output_arrays @vdesc number of passed input arrays @vtype int @vio in @endvar @var interp_coord_array_types @vdesc types of passed coordinates to interpolate at - @vtype int [num_dims] + @vtype int [N_dims] @vio in @endvar @@ -1043,32 +994,32 @@ static int GetLocalCoords (cGH *GH, -1 for failure (wrong parameters passed or limitations not met) @endreturndesc @@*/ -static int PUGHInterp_CheckArguments (cGH *GH, - int num_dims, - int num_points, - int num_in_arrays, - int num_out_arrays, - const int interp_coord_array_types[]) +static int CheckArguments (const cGH *GH, + int N_dims, + int N_points, + int N_input_arrays, + int N_output_arrays, + const int interp_coord_array_types[]) { int i; /* check for invalid arguments */ - if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0) + if (N_dims < 0 || N_points < 0 || N_input_arrays < 0 || N_output_arrays < 0) { return (-1); } /* check if there's anything to do at all */ - /* NOTE: num_points can be 0 in a collective call */ - if (num_dims == 0 || (CCTK_nProcs (GH) == 1 && num_points == 0) || - num_in_arrays == 0 || num_out_arrays == 0) + /* NOTE: N_points can be 0 in a collective call */ + if (N_dims == 0 || (CCTK_nProcs (GH) == 1 && N_points == 0) || + N_input_arrays == 0 || N_output_arrays == 0) { return (0); } /* for now we can only deal with coordinates of type CCTK_REAL */ - for (i = 0; i < num_dims; i++) + for (i = 0; i < N_dims; i++) { if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL) { @@ -1079,7 +1030,7 @@ static int PUGHInterp_CheckArguments (cGH *GH, /* PUGHInterp's interpolation operators compute one output array per input array */ - if (num_in_arrays != num_out_arrays) + if (N_input_arrays != N_output_arrays) { CCTK_WARN (1, "Number of input arrays must match number of output arrays"); return (-1); diff --git a/src/Startup.c b/src/Startup.c index 4ff970d..19cf382 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -125,13 +125,13 @@ static void *SetupGH (tFleshConfig *config, int convergence_level, cGH *GH) /* allocate once for all fields of same type */ nprocs = CCTK_nProcs (GH); - myGH->send_count = malloc (4 * nprocs * sizeof (int)); - myGH->send_displ = myGH->send_count + nprocs; - myGH->recv_count = myGH->send_displ + nprocs; - myGH->recv_displ = myGH->recv_count + nprocs; + myGH->sendcnt = malloc (4 * nprocs * sizeof (int)); + myGH->senddispl = myGH->sendcnt + nprocs; + myGH->recvcnt = myGH->senddispl + nprocs; + myGH->recvdispl = myGH->recvcnt + nprocs; - myGH->num_points_from = malloc (2 * nprocs * sizeof (CCTK_INT)); - myGH->num_points_to = myGH->num_points_from + nprocs; + myGH->N_points_from = malloc (2 * nprocs * sizeof (CCTK_INT)); + myGH->N_points_to = myGH->N_points_from + nprocs; return (myGH); } @@ -143,7 +143,7 @@ static void *SetupGH (tFleshConfig *config, int convergence_level, cGH *GH) @date Wed 14 Feb 2001 @author Thomas Radke @desc - Wrappers for the different interpolation operators + Wrappers for the different interpolation operators registered for first/second/third order interpolation. These wrappers just call the common interpolation routine passing all arguments plus the interpolation order. @@ -170,7 +170,7 @@ static int InterpGV_1stOrder (cGH *GH, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types)); } - + static int InterpGV_2ndOrder (cGH *GH, const char *coord_system, @@ -188,7 +188,7 @@ static int InterpGV_2ndOrder (cGH *GH, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types)); } - + static int InterpGV_3rdOrder (cGH *GH, const char *coord_system, diff --git a/src/pughInterpGH.h b/src/pughInterpGH.h index 23dc6ef..543d1d5 100644 --- a/src/pughInterpGH.h +++ b/src/pughInterpGH.h @@ -3,28 +3,33 @@ @date Sun Jul 04 1999 @author Thomas Radke @desc - Definition of GH extensions of thorn PUGHInterp and - declaration of function prototypes + Definition of GH extensions structure for thorn PUGHInterp + and declaration of function prototypes. @enddesc - @history - @endhistory + @version $Header$ @@*/ +#ifndef _PUGHINTERP_PUGHINTERP_H_ +#define _PUGHINTERP_PUGHINTERP_H_ 1 + +#ifdef __cplusplus +extern "C" +{ +#endif typedef struct { - int *send_count; /* number of elements to be sent */ - int *send_displ; /* offset of each element to be sent */ - int *recv_count; /* number of elements to be received */ - int *recv_displ; /* offset of each element to be received */ + int *sendcnt; /* number of elements to be sent */ + int *senddispl; /* offset of each element to be sent */ + int *recvcnt; /* number of elements to be received */ + int *recvdispl; /* offset of each element to be received */ - CCTK_INT *num_points_to; /* number of coordinate points to be sent */ - CCTK_INT *num_points_from; /* number of coordinate points to be received */ + CCTK_INT *N_points_to; /* number of coordinate points to be sent */ + CCTK_INT *N_points_from; /* number of coordinate points to be received */ - int *whichproc; /* which processor does point i belong to */ - int *indices; /* indices to sort from processor-sorted into - point-sorted order */ -} pughInterpGH; + int *whichproc; /* which processor does point i belong to */ + int *indices; /* indices to sort from processor-sorted into */ +} pughInterpGH; /* point-sorted order */ /* prototypes of interpolation operators to be registered */ @@ -53,3 +58,9 @@ int PUGHInterp_Interpolate (int order, const void *const in_arrays[ /* num_arrays */ ], const int out_types[ /* num_arrays */ ], void *const out_arrays[ /* num_arrays */ ]); + +#ifdef __cplusplus +} // extern "C" +#endif + +#endif /* _PUGHINTERP_PUGHINTERP_H_ */ -- cgit v1.2.3