diff options
Diffstat (limited to 'src/pugh/pGF_Reduction.c')
| -rw-r--r-- | src/pugh/pGF_Reduction.c | 198 |
1 files changed, 198 insertions, 0 deletions
diff --git a/src/pugh/pGF_Reduction.c b/src/pugh/pGF_Reduction.c new file mode 100644 index 0000000..bc8e211 --- /dev/null +++ b/src/pugh/pGF_Reduction.c @@ -0,0 +1,198 @@ + /*@@ + @file pGF_Reduction.c + @date Thu Apr 3 11:54:53 1997 + @author Paul Walker + @desc + Various MPI reduction operators for maxima, minima, norms + and the like on grid functions. These are all pretty + straighforward thanks to the miracle of MPI_Allreduce. + @enddesc + @version $Id$ + @@*/ + +#include "pugh.h" + +static char *rcsid = "$Id$"; + + /*@@ + @routine pGF_MaxVal + @date Thu Apr 3 11:55:25 1997 + @author Paul Walker + @desc + Returns the maximum value of a distributed grid function. + @enddesc +@@*/ + + +Double pGF_MaxVal(pGH *GH, pGF *GF) { + Double res, lres; + int i; + + if (!GF->storage) return GF->maxval; + + res = GF->data[0]; + for (i=1;i<GH->npoints;i++) + res = (GF->data[i] > res? GF->data[i] : res); + + + /* So now res is the local max */ +#ifdef MPI + lres = res; + CACTUS_MPI_ERROR(MPI_Allreduce(&lres,&res,1,PUGH_MPI_TYPE,MPI_MAX,GH->PUGH_COMM_WORLD)); +#endif + + return res; +} + + + /*@@ + @routine pGF_MinVal + @date Thu Apr 3 11:55:48 1997 + @author Paul Walker + @desc + Returns the minimum value of a distributed grid function. + @enddesc +@@*/ + +Double pGF_MinVal(pGH *GH, pGF *GF) { + Double res, lres; + int i; + + if (!GF->storage) return GF->minval; + + res = GF->data[0]; + for (i=1;i<GH->npoints;i++) + res = (GF->data[i] < res ? GF->data[i] : res); + + /* So now res is the local max */ +#ifdef MPI + lres = res; + CACTUS_MPI_ERROR(MPI_Allreduce(&lres,&res,1,PUGH_MPI_TYPE,MPI_MIN,GH->PUGH_COMM_WORLD)); +#endif + + return res; +} + + + + /*@@ + @routine pGF_Norm1 + @date Thu Apr 3 11:56:05 1997 + @author Paul Walker + @desc + Returns the average of a distributed grid function. + @enddesc +@@*/ + +Double pGF_Norm1(pGH *GH, pGF *GF) { + Double res, lres; + int i, ii, jj, kk; + int is,ie,js,je,ks,ke; + int np,tnp; + + if (!GF->storage) return GF->norm1; + + +#define ABS(x) ((x)<0?-(x):(x)) + +#ifndef MPI + res = 0.0; + for (i=0;i<GH->npoints;i++) res += ABS(GF->data[i]); + tnp = GH->npoints; + +#else + res = 0.0; + np = 0; + for (kk=GH->ownership[GF->stagger][2][0]; + kk < GH->ownership[GF->stagger][2][1]; + kk ++) + for (jj = GH->ownership[GF->stagger][1][0]; + jj < GH->ownership[GF->stagger][1][1]; + jj ++) + for (ii = GH->ownership[GF->stagger][0][0]; + ii < GH->ownership[GF->stagger][0][1]; + ii ++) { + res += ABS(GF->data[DATINDEX(GH,ii,jj,kk)]); + np++; + } + lres = res; + CACTUS_MPI_ERROR(MPI_Allreduce(&lres,&res,1,PUGH_MPI_TYPE,MPI_SUM,GH->PUGH_COMM_WORLD)); + CACTUS_MPI_ERROR(MPI_Allreduce(&np,&tnp,1,MPI_INT,MPI_SUM,GH->PUGH_COMM_WORLD)); +#endif + + res = res/tnp; + return res; +} + + + /*@@ + @routine pGF_Norm2 + @date Thu Apr 3 11:56:36 1997 + @author Paul Walker + @desc + Returns the "norm2" of a distributed grid function. That + is, returns $\sqrt{\Sigma (a_i * a_i) / np}$. + @enddesc +@@*/ + + +Double pGF_Norm2(pGH *GH, pGF *GF) { + Double res, lres; + int i, ii, jj, kk; + int is,ie,js,je,ks,ke; + int np, tnp; + + if (!GF->storage) return GF->norm2; + +#ifndef MPI + res = 0.0; + for (i=0;i<GH->npoints;i++) res += GF->data[i]*GF->data[i]; + tnp = GH->npoints; + +#else + /* MPI version depends heavily on the ownership array */ + res = 0.0; + np = 0; + + for (kk=GH->ownership[GF->stagger][2][0]; + kk < GH->ownership[GF->stagger][2][1]; + kk ++) + for (jj = GH->ownership[GF->stagger][1][0]; + jj < GH->ownership[GF->stagger][1][1]; + jj ++) + for (ii = GH->ownership[GF->stagger][0][0]; + ii < GH->ownership[GF->stagger][0][1]; + ii ++) { + res += GF->data[DATINDEX(GH,ii,jj,kk)] * + GF->data[DATINDEX(GH,ii,jj,kk)]; + np ++; + } + lres = res; + CACTUS_MPI_ERROR(MPI_Allreduce(&lres,&res,1,PUGH_MPI_TYPE,MPI_SUM,GH->PUGH_COMM_WORLD)); + CACTUS_MPI_ERROR(MPI_Allreduce(&np,&tnp,1,MPI_INT,MPI_SUM,GH->PUGH_COMM_WORLD)); +#endif + + res = sqrt(res/tnp); + return res; +} + + /*@@ + @routine pGF_AllNorms + @date Mon Jun 30 10:44:57 1997 + @author Paul Walker + @desc + Returns all 4 norms into a Double pointer in the order + max min norm1 norm2 + <p> + This is purely a convenience function. + @enddesc + @calls pGF_MaxVal, pGF_MinVal, pGF_Norm1, pGF_Norm2 +@@*/ + +void pGF_AllNorms(pGH *GH, pGF *GF, Double *res) { + res[0] = pGF_MaxVal(GH,GF); + res[1] = pGF_MinVal(GH,GF); + res[2] = pGF_Norm1(GH,GF); + res[3] = pGF_Norm2(GH,GF); +} + |