fmm3d_check_mpi.cpp

Go to the documentation of this file.
00001 /* Kernel Independent Fast Multipole Method
00002    Copyright (C) 2004 Lexing Ying, New York University
00003 
00004 This program is free software; you can redistribute it and/or modify
00005 it under the terms of the GNU General Public License as published by
00006 the Free Software Foundation; either version 2, or (at your option)
00007 any later version.
00008 
00009 This program is distributed in the hope that it will be useful, but WITHOUT
00010 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00011 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00012 for more details.
00013 
00014 You should have received a copy of the GNU General Public License
00015 along with this program; see the file COPYING.  If not, write to the Free
00016 Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
00017 02111-1307, USA.  */
00018 #include "fmm3d_mpi.hpp"
00019 #include "common/vecmatop.hpp"
00020 
00027 /* ********************************************************************** */
00028 /* Check takes the source densities, target values as solved for and a value
00029  * numChk, which is the point being checked.  The value rerr (relative error)
00030  * is returned.
00031  */
00032 #undef __FUNCT__
00033 #define __FUNCT__ "FMM3d_MPI::check"
00034 int FMM3d_MPI::check(Vec srcDen, Vec trgVal, int numChk, double& rerr)
00035 {
00036   //begin
00037   pA(procLclNum(_trgPos)>0); //have things to check
00038   //pC( PetscPrintf(mpiComm(), "checking.........\n") );
00039   int dim = this->dim();
00040   int srcDOF = this->srcDOF();
00041   int trgDOF = this->trgDOF();
00042   int slclnum = procLclNum(_srcPos);
00043   int tlclnum = procLclNum(_trgPos);   //int lclnum = this->lclnum();
00044   
00045   //1. select point to check
00046   int ttlChk = numChk*mpiSize();
00047   vector<int> chkVec(numChk);
00048   for(int k=0; k<numChk; k++) {
00049          chkVec[k] = int( floor(drand48()*tlclnum) );
00050          pA(chkVec[k]>=0 && chkVec[k]<tlclnum);
00051   }
00052   pC( MPI_Barrier(mpiComm()) );
00053   
00054   Vec chkPos; pC( VecCreateMPI(mpiComm(), numChk*dim,  ttlChk*dim,  &chkPos) ); //check position
00055   Vec chkVal; pC( VecCreateMPI(mpiComm(), numChk*trgDOF, ttlChk*trgDOF, &chkVal) ); //check value
00056   Vec chkDen; pC( VecCreateMPI(mpiComm(), numChk*trgDOF, ttlChk*trgDOF, &chkDen) ); //check density
00057   
00058   double* posArr; pC( VecGetArray(_trgPos, &posArr) );
00059   double* valArr; pC( VecGetArray( trgVal, &valArr) );
00060   double* chkPosArr; pC( VecGetArray(chkPos, &chkPosArr) );
00061   double* chkValArr; pC( VecGetArray(chkVal, &chkValArr) );
00062   for(int k=0; k<numChk; k++) {
00063          for(int i=0; i<dim; i++)
00064                 chkPosArr[ k*dim+i  ] = posArr[ chkVec[k]*dim+i  ];
00065          for(int i=0; i<trgDOF;i++)
00066                 chkValArr[ k*trgDOF+i ] = valArr[ chkVec[k]*trgDOF+i ];
00067   } 
00068   pC( VecRestoreArray(_trgPos, &posArr) );
00069   pC( VecRestoreArray( trgVal, &valArr) );
00070   pC( VecRestoreArray(chkPos, &chkPosArr) );
00071   pC( VecRestoreArray(chkVal, &chkValArr) );
00072   
00073   //2. compute and gather
00074   Vec allChkPos;
00075   {
00076          VecScatter ctx;  pC( VecScatterCreateToAll(chkPos, &ctx, &allChkPos) );
00077          pC( VecScatterBegin(chkPos, allChkPos, INSERT_VALUES, SCATTER_FORWARD, ctx) );
00078          pC( VecScatterEnd(  chkPos, allChkPos, INSERT_VALUES, SCATTER_FORWARD, ctx) );
00079          pC( VecScatterDestroy(ctx) );
00080   }
00081   Vec allChkDen; pC( VecCreateSeq(PETSC_COMM_SELF, ttlChk*trgDOF, &allChkDen) );
00082   Vec glbChkDen; pC( VecCreateSeq(PETSC_COMM_SELF, ttlChk*trgDOF, &glbChkDen) );
00083   
00084   double* lclSrcPosArr; pC( VecGetArray(_srcPos, &lclSrcPosArr) );
00085   double* lclSrcNorArr; pC( VecGetArray(_srcNor, &lclSrcNorArr) );
00086   double* lclSrcDenArr; pC( VecGetArray( srcDen, &lclSrcDenArr) );
00087   double* allChkPosArr; pC( VecGetArray(allChkPos, &allChkPosArr) );
00088   double* allChkDenArr; pC( VecGetArray(allChkDen, &allChkDenArr) );
00089   double* glbChkDenArr; pC( VecGetArray(glbChkDen, &glbChkDenArr) );
00090   
00091   DblNumMat lclSrcPosMat(dim, slclnum, false, lclSrcPosArr);
00092   DblNumMat lclSrcNorMat(dim, slclnum, false, lclSrcNorArr);
00093   DblNumVec lclSrcDenVec(srcDOF*slclnum, false, lclSrcDenArr);
00094   DblNumMat allChkPosMat(dim,  ttlChk, false, allChkPosArr);
00095   DblNumVec allChkDenVec(trgDOF* ttlChk, false, allChkDenArr);
00096   
00097   DblNumMat inter(trgDOF, slclnum*srcDOF);
00098   for(int i=0; i<ttlChk; i++) {
00099          DblNumMat oneChkPosMat(dim, 1, false, allChkPosMat.clmdata(i));
00100          DblNumVec oneChkDenVec(trgDOF, false, allChkDenVec.data()+i*trgDOF);
00101          pC( _knl.buildKnlIntCtx(lclSrcPosMat, lclSrcNorMat, oneChkPosMat, inter) );
00102          pC( dgemv(1.0, inter, lclSrcDenVec, 0.0, oneChkDenVec) );
00103   }
00104   
00105   pC( MPI_Barrier(mpiComm()) );
00106   pC( MPI_Allreduce(allChkDenArr, glbChkDenArr, ttlChk*trgDOF, MPI_DOUBLE, MPI_SUM, mpiComm()) );
00107   
00108   //3. distribute to individual
00109   double* lclChkDenArr; pC( VecGetArray(chkDen, &lclChkDenArr) );
00110   void* dst = (void*)lclChkDenArr;
00111   void* src = (void*)(glbChkDenArr + mpiRank()*numChk*trgDOF);  
00112   memcpy(dst, src, sizeof(double)*numChk*trgDOF);
00113   pC( VecRestoreArray(chkDen, &lclChkDenArr) );
00114   
00115   pC( VecRestoreArray(_srcPos, &lclSrcPosArr) );
00116   pC( VecRestoreArray(_srcNor, &lclSrcNorArr) );
00117   pC( VecRestoreArray( srcDen, &lclSrcDenArr) );
00118   pC( VecRestoreArray(allChkPos, &allChkPosArr) );
00119   pC( VecRestoreArray(allChkDen, &allChkDenArr) );
00120   pC( VecRestoreArray(glbChkDen, &glbChkDenArr) );
00121   
00122   pC( VecDestroy(allChkPos) );
00123   pC( VecDestroy(allChkDen) );
00124   pC( VecDestroy(glbChkDen) );
00125   
00126   PetscScalar mone = -1.0;
00127   pC( VecAXPY(chkVal, mone, chkDen) ); //chkVal now is the error - HARPER:  In correct order??
00128   double nerr, npot;
00129   pC( VecNorm(chkVal, NORM_2, &nerr) );
00130   pC( VecNorm(chkDen, NORM_2, &npot) );
00131   
00132   rerr = nerr/npot;//pC( PetscPrintf(mpiComm(), "Error %e, %e %e\n", nerr/npot, nerr, npot) );
00133   
00134   pC( VecDestroy(chkPos) );
00135   pC( VecDestroy(chkVal) );
00136   pC( VecDestroy(chkDen) );
00137   
00138   return(0);
00139 }
00140 
00141 
00142 

Generated on Sun Dec 4 21:12:40 2005 for fmm3d_mpi by  doxygen 1.4.5