int FMM3d::check const DblNumVec &  srcDen,
DblNumVec &  trgVal,
int  numChk,
double &  relativeErr
 

Check relative error - in fmm_check.cpp

Definition at line 25 of file fmm3d_check.cpp.

References KnlMat3d::_srcPos, KnlMat3d::_trgPos, KnlMat3d::dim(), KnlMat3d::srcDOF(), and KnlMat3d::trgDOF().

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 }


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