/*@@ @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_Reduction_c) /* some useful macros */ #ifdef MIN #undef MIN #endif #ifdef MAX #undef MAX #endif #ifdef ABS #undef ABS #endif #ifdef SQR #undef SQR #endif #define MIN(x, y) ((x) < (y) ? (x) : (y)) #define MAX(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 @@*/ void PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [], int nArrays, void *inArrays [], int inType, int nOutVals, void *outVals, int outType) { 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; } nPoints = 1; } if ((elemSize = getSizeofElement (outType)) < 0) { CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type"); return; } 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); } } free (buffer); } /*@@ @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 @@*/ void PUGH_ReduceNorm1 (cGH *GH, int proc, int nOutVals, int outType, void *outVals, int numinvars, int *varlist) { int i, elemSize; pGH *pughGH; CCTK_REAL intermediate_result; outVals = NULL; if (nOutVals != 1) { CCTK_WARN (1, "PUGH_ReduceNorm1 accepts only one output value per " "input array"); return; } if ((elemSize = getSizeofElement (outType)) < 0) { CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type"); return; } pughGH = PUGH_pGH (GH); for (i = 0; i < numinvars; i++) { int result; result = 1; if (CCTK_GroupDimFromVarI(varlist[i])==1) { 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); } } } /*********************************************************************/ /* 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); }