fmm3d_check.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.hpp"
00019 #include "common/vecmatop.hpp"
00020 
00021 using std::cerr;
00022 using std::endl;
00023 
00024 /* ********************************************************************** */
00025 int FMM3d::check(const DblNumVec& srcden, DblNumVec& trgVal, int numChk, double& rerr)
00026 {
00027   iA(_trgPos->n()>0); //have things to check
00028   int dim = this->dim();
00029   int srcDOF = this->srcDOF();
00030   int trgDOF = this->trgDOF();
00031   int srcNum = _srcPos->n();
00032   int trgNum = _trgPos->n();   //int lclnum = this->lclnum();
00033   
00034   //1. select point to check
00035   vector<int> chkVec(numChk);
00036   for(int k=0; k<numChk; k++) {
00037          chkVec[k] = int( floor(drand48()*trgNum) );     iA(chkVec[k]>=0 && chkVec[k]<trgNum);
00038   }
00039   DblNumMat ChkPos(dim, numChk);
00040   DblNumVec chkVal(numChk*trgDOF);
00041   DblNumVec chkDen(numChk*trgDOF);
00042   for(int k=0; k<numChk; k++) {
00043          for(int i=0; i<dim; i++)               ChkPos(i, k) = (*_trgPos)(i, chkVec[k]);
00044          for(int i=0; i<trgDOF;i++)             chkVal(k*trgDOF+i) = trgVal(chkVec[k]*trgDOF+i);
00045   }
00046   
00047   //2. compute and gather
00048   DblNumMat inter(trgDOF, srcNum*srcDOF);
00049   for(int i=0; i<numChk; i++) {
00050          DblNumMat oneChkPos(dim, 1, false, ChkPos.clmdata(i));
00051          DblNumVec onechkDen(trgDOF, false, chkDen.data()+i*trgDOF);
00052          iC( _knl.kernel(*_srcPos, *_srcNor, oneChkPos, inter) );
00053          iC( dgemv(1.0, inter, srcden, 0.0, onechkDen) );
00054   }
00055   
00056   //3. distribute to individual
00057   for(int k=0; k<numChk; k++)
00058          for(int i=0; i<trgDOF; i++)
00059                 chkDen(k*trgDOF+i) -= chkVal(k*trgDOF+i);
00060   
00061   double vn = 0;  double en = 0;
00062   for(int k=0; k<numChk; k++)
00063          for(int i=0; i<trgDOF; i++) {
00064                 vn += chkVal(k*trgDOF+i) * chkVal(k*trgDOF+i);
00065                 en += chkDen(k*trgDOF+i) * chkDen(k*trgDOF+i);
00066          }
00067   vn = sqrt(vn);
00068   en = sqrt(en);
00069   
00070   //cerr<<"relative error: "<<en/vn<<endl;  cerr<<"error norm    : "<<en<<endl;  cerr<<"solution norm : "<<vn<<endl;
00071   rerr = en/vn; //relative error
00072   
00073   return 0;
00074 }
00075 
00076 
00077 

Generated on Sun Dec 4 19:24:39 2005 for fmm3d by  doxygen 1.4.5