From 648dec4fea36e44e38c9acd00d17def2ff74f208 Mon Sep 17 00:00:00 2001 From: allen Date: Wed, 19 Jul 2000 12:43:53 +0000 Subject: Registered reduction routines return an error code, and split reduction.c into several files git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGH/trunk@243 b61c5cb5-eaca-4651-9a7a-d64986f99364 --- src/ReductionNorm1.c | 650 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 650 insertions(+) create mode 100644 src/ReductionNorm1.c (limited to 'src/ReductionNorm1.c') diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c new file mode 100644 index 0000000..94a9d3f --- /dev/null +++ b/src/ReductionNorm1.c @@ -0,0 +1,650 @@ + /*@@ + @file Reduction_nm1.c + @date Thu Apr 3 11:54:53 1997 + @author Thomas Radke, Paul Walker + @desc + Various MPI reduction operators. + @enddesc + @version $Header$ + @@*/ + +#include +#include +#include + +#include "cctk.h" +#include "pugh.h" + +static char *rcsid = "$Id$"; + +CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm1_c) + +/* some useful macros */ +#ifdef MIN +#undef MIN +#endif +#ifdef ABS +#undef ABS +#endif +#ifdef SQR +#undef SQR +#endif +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#define ABS(x) ((x) < 0 ?-(x) : (x)) +#define SQR(x) ((x) * (x)) + + +/* local function prototypes */ +int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal); +int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints, + void *inArray, int inType, CCTK_REAL outVals []); +void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n); +int getSizeofElement (int varType); + + +/*@@ + @routine PUGH_ReduceArrayNorm1 + @author Thomas Radke + @date 19 Aug 1999 + @desc + Registered PUGH reduction routine for computing the "norm1" + of a set of arrays. + The arrays are described by their dimensions and variable type. + If the number of requested output values equals the number of + elements in an input array, global reduction is done element-wise. + Otherwise nOutVals must be 1, and global reduction is done on the + results of the local reductions. + Type casting of the result is provided by specifying the + requested output datatype. The intermediate reduction value + is always computed as a CCTK_REAL value internally. + @enddesc + @history + @endhistory + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH + @vio in + @endvar + @var proc + @vdesc processor that should receive the result of operation + (negative value means all processors receive the result) + @vtype int + @vio in + @endvar + @var nDims + @vdesc number of dimensions of input arrays + @vtype int + @vio in + @endvar + @var dims + @vdesc dimensions of input arrays + @vtype int * + @vio in + @endvar + @var nArrays + @vdesc number of input arrays + @vtype int + @vio in + @endvar + @var inArrays + @vdesc field of input arrays + @vtype void ** + @vio in + @endvar + @var inType + @vdesc (common) variable type of input arrays + @vtype int + @vio in + @endvar + @var nOutVals + @vdesc number of values per output array + @vtype int + @vio in + @endvar + @var outVals + @vdesc pointer to buffer holding the output values + @vtype void * + @vio in + @endvar + @var outType + @vdesc (common) variable type of output arrays + @vtype int + @vio in + @endvar +@@*/ +int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [], + int nArrays, void *inArrays [], int inType, + int nOutVals, void *outVals, int outType) +{ + int retval = 1; + int i, elemSize, nPoints; + CCTK_REAL *buffer; + + + nPoints = dims [0]; + for (i = 1; i < nDims; i++) + { + nPoints *= dims [i]; + } + + if (nOutVals != 1) + { + if (nOutVals != nPoints) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "PUGH_ReduceArrayNorm1: Don't know how to reduce " + "a %d-dimensional array with " + "%d elements to an output array of %d elements", + nDims, nPoints, nOutVals); + return -1; + } + nPoints = 1; + } + + if ((elemSize = getSizeofElement (outType)) < 0) + { + CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type"); + return -1; + } + + buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL)); + + for (i = 0; i < nArrays; i++) + { + + /* do the reduction on input array [i] */ + if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType, + buffer) == 0) + { + + /* type-cast the result to the requested datatype */ + if (proc < 0 || proc == CCTK_MyProc (GH)) + { + copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize, + outType, nOutVals); + retval = 0; + } + } + else + { + retval = -1; + } + } + + free (buffer); + return retval; +} + + + + +/*@@ + @routine PUGH_ReduceNorm1 + @author Thomas Radke + @date 23 Apr 1999 + @desc + PUGH reduction routine to be registered for calculating the norm1. + It just goes through the list of variables and calls the appropriate + grouptype reduction routine. + Global reduction is always done on the results of the local + reductions, so nOutVals must be 1. + @enddesc + @history + @endhistory + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH + @vio in + @endvar + @var proc + @vdesc processor that should receive the result of operation + (negative value means all processors receive the result) + @vtype int + @vio in + @endvar + @var nOutVars + @vdesc number of requested output elements + @vtype int + @vio in + @endvar + @var outType + @vdesc type of return value + @vtype int + @vio in + @endvar + @var outVals + @vdesc array for the result values to be stored + @vtype void * + @vio in + @endvar + @var numinvars + @vdesc number of variables in the list + @vtype int + @vio in + @endvar + @var varlist + @vdesc list of variables (given as their global indices) + @vtype int * + @vio in + @endvar +@@*/ +int PUGH_ReduceNorm1 (cGH *GH, int proc, + int nOutVals, + int outType, + void *outVals, + int numinvars, + int *varlist) +{ + int retval=1; + int i, elemSize; + pGH *pughGH; + CCTK_REAL intermediate_result; + + if (nOutVals != 1) + { + CCTK_WARN (1, "PUGH_ReduceNorm1: Only one output value per " + "input array allowed"); + return -2; + } + + if ((elemSize = getSizeofElement (outType)) < 0) + { + CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type"); + return -3; + } + + pughGH = PUGH_pGH (GH); + + for (i = 0; i < numinvars; i++) + { + int result; + + result = 1; + + if (CCTK_GroupDimFromVarI(varlist[i])<3) + { + CCTK_VWarn(1,__LINE__,__FILE__,"PUGH", + "PUGH_ReduceNorm1: Reducing dimension %d not yet supported", + CCTK_GroupDimFromVarI(varlist[i])); + } + else + { + + switch (CCTK_GroupTypeFromVarI (varlist [i])) + { + case GROUP_GF: + result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result); + break; + case GROUP_ARRAY: + CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported"); + break; + case GROUP_SCALAR: + CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense"); + break; + default: + CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type"); + break; + } + } + + if (result == 0 && (proc < 0 || proc == pughGH->myproc)) + { + copy_real_to_outtype (&intermediate_result, (char *) outVals + + i*nOutVals*elemSize, outType, nOutVals); + retval = 0; + } + else if (result == 0) + { + retval = -1; + } + } + + return retval; +} + + + +/*********************************************************************/ +/* local functions */ +/*********************************************************************/ + + +/*@@ + @routine PUGH_Norm1_Array + @date Aug 19 1999 + @author Thomas Radke + @desc Returns the "norm1" of a distributed array with + 'nPoints' elements. Global reduction is done element-wise + (nOutVals == 1) or on the results of the local reductions. + + "norm1" is defined as $\Sigma |a_i| / np$. + @enddesc +@@*/ +int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints, + void *inArray, int inType, CCTK_REAL outVals []) +{ + int i, j; + pGH *pughGH = PUGH_pGH (GH); + + + for (i = 0; i < nOutVals; i++) + { + outVals [i] = 0.0; + for (j = 0; j < nPoints; j++) + { + if (inType == CCTK_VARIABLE_CHAR) /* CCTK_CHAR is unsigned */ + { + outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j]; + } + else if (inType == CCTK_VARIABLE_INT) + { + outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]); + } +#ifdef CCTK_INT2 + else if (inType == CCTK_VARIABLE_INT2) + { + outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]); + } +#endif +#ifdef CCTK_INT4 + else if (inType == CCTK_VARIABLE_INT4) + { + outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]); + } +#endif +#ifdef CCTK_INT8 + else if (inType == CCTK_VARIABLE_INT8) + { + outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]); + } +#endif + else if (inType == CCTK_VARIABLE_REAL) + { + outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]); + } +#ifdef CCTK_REAL4 + else if (inType == CCTK_VARIABLE_REAL4) + { + outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]); + } +#endif +#ifdef CCTK_REAL8 + else if (inType == CCTK_VARIABLE_REAL8) + { + outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]); + } +#endif +#ifdef CCTK_REAL16 + else if (inType == CCTK_VARIABLE_REAL16) + { + outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]); + } +#endif + else if (inType == CCTK_VARIABLE_COMPLEX) + { + outVals [i] += + sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+ + SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) ); + } + else + { + CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type"); + return (1); + } + } + } + +#ifdef CCTK_MPI + + if (pughGH->nprocs > 1) + { + CCTK_REAL *localOutVals; + + localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL)); + + /* outVals[] contains now the local sum */ + memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL)); + if (proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals, + PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals, + PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD)); + } + free (localOutVals); + + nPoints *= pughGH->nprocs; + } + +#endif + + if (proc < 0 || proc == pughGH->myproc) + { + for (i = 0; i < nOutVals; i++) + { + outVals [i] /= nPoints; + } + } + + return (0); +} + + + +/*@@ + @routine PUGH_Norm1_GF + @date Thu Apr 3 11:56:05 1997 + @author Paul Walker + @desc + Returns the average of a distributed grid function. + @enddesc +@@*/ +int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal) +{ + int timelevel; + int ii; + CCTK_INT tnp; + pGH *pughGH; + pGA *GA; +#ifdef CCTK_MPI + int jj, kk; + CCTK_INT np; + CCTK_REAL local_result; +#endif + + timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1; + if (timelevel > 0) + { + timelevel--; + } + + pughGH = PUGH_pGH (GH); + GA = ((pGA ***) pughGH->variables) [index][timelevel]; + + *outVal = 0.0; + +#ifndef CCTK_MPI + for (ii = 0; ii < GA->extras->npoints; ii++) + { + if (GA->vtype == CCTK_VARIABLE_CHAR) + { + *outVal += ((CCTK_CHAR *) GA->data) [ii]; /* CCTK_CHAR is unsigned */ + } + else if (GA->vtype == CCTK_VARIABLE_INT) + { + *outVal += ABS (((CCTK_INT *) GA->data) [ii]); + } +#ifdef CCTK_INT2 + else if (GA->vtype == CCTK_VARIABLE_INT2) + { + *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]); + } +#endif +#ifdef CCTK_INT4 + else if (GA->vtype == CCTK_VARIABLE_INT4) + { + *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]); + } +#endif +#ifdef CCTK_INT8 + else if (GA->vtype == CCTK_VARIABLE_INT8) + { + *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]); + } +#endif + else if (GA->vtype == CCTK_VARIABLE_REAL) + { + *outVal += ABS (((CCTK_REAL *) GA->data) [ii]); + } +#ifdef CCTK_REAL4 + else if (GA->vtype == CCTK_VARIABLE_REAL4) + { + *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]); + } +#endif +#ifdef CCTK_REAL8 + else if (GA->vtype == CCTK_VARIABLE_REAL8) + { + *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]); + } +#endif +#ifdef CCTK_REAL16 + else if (GA->vtype == CCTK_VARIABLE_REAL16) + { + *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]); + } +#endif + else if (GA->vtype == CCTK_VARIABLE_COMPLEX) + { + *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+ + SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) ); + } + else + { + CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type"); + return (1); + } + } + + tnp = GA->extras->npoints; + +#else + + /* MPI version depends heavily on the ownership array */ + np = 0; + local_result = 0.0; + + for (kk=GA->extras->ownership[GA->stagger][0][2]; + kk < GA->extras->ownership[GA->stagger][1][2]; + kk ++) + { + for (jj = GA->extras->ownership[GA->stagger][0][1]; + jj < GA->extras->ownership[GA->stagger][1][1]; + jj ++) + { + for (ii = GA->extras->ownership[GA->stagger][0][0]; + ii < GA->extras->ownership[GA->stagger][1][0]; + ii++, np++) + { + if (GA->vtype == CCTK_VARIABLE_CHAR) + { + local_result += /* CCTK_CHAR is unsigned */ + ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]; + } + else if (GA->vtype == CCTK_VARIABLE_INT) + { + local_result += + ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#ifdef CCTK_INT2 + else if (GA->vtype == CCTK_VARIABLE_INT2) + { + local_result += + ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif +#ifdef CCTK_INT4 + else if (GA->vtype == CCTK_VARIABLE_INT4) + { + local_result += + ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif +#ifdef CCTK_INT8 + else if (GA->vtype == CCTK_VARIABLE_INT8) + { + local_result += + ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif + else if (GA->vtype == CCTK_VARIABLE_REAL) + { + local_result += + ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#ifdef CCTK_REAL4 + else if (GA->vtype == CCTK_VARIABLE_REAL4) + { + local_result += + ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif +#ifdef CCTK_REAL8 + else if (GA->vtype == CCTK_VARIABLE_REAL8) + { + local_result += + ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif +#ifdef CCTK_REAL16 + else if (GA->vtype == CCTK_VARIABLE_REAL16) + { + local_result += + ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]); + } +#endif + else if (GA->vtype == CCTK_VARIABLE_COMPLEX) + { + local_result += + sqrt( SQR((((CCTK_COMPLEX *)GA->data) + [DATINDEX (GA->extras,ii,jj,kk)]).Re)+ + SQR((((CCTK_COMPLEX *)GA->data) + [DATINDEX (GA->extras,ii,jj,kk)]).Im) ); + } + else + { + CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type"); + return (1); + } + } + } + } + + if (proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL, + MPI_SUM, pughGH->PUGH_COMM_WORLD)); + CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM, + pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL, + MPI_SUM, proc, pughGH->PUGH_COMM_WORLD)); + CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM, + proc, pughGH->PUGH_COMM_WORLD)); + } +#endif + + if (proc < 0 || proc == pughGH->myproc) + { + *outVal /= tnp; + } + + return (0); +} + + -- cgit v1.2.3