fmm3d_eval_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 
00021 using std::cerr;
00022 using std::endl;
00023 
00024 // ---------------------------------------------------------------------- 
00025 #undef __FUNCT__
00026 #define __FUNCT__ "FMM3d_MPI::evaluate"
00027 int FMM3d_MPI::evaluate(Vec srcDen, Vec trgVal)
00028 {
00029   //begin  //ebiLogInfo( "multiply.............");
00030   //-----------------------------------
00031   //cerr<<"fmm src and trg numbers "<<pglbnum(_srcPos)<<" "<<pglbnum(_trgPos)<<endl;
00032   int tmp;
00033   pC( VecGetSize(srcDen,&tmp) );  pA(tmp==srcDOF()*procGlbNum(_srcPos));
00034   pC( VecGetSize(trgVal,&tmp) );  pA(tmp==trgDOF()*procGlbNum(_trgPos));
00035   
00036   int srcDOF = this->srcDOF();
00037   int trgDOF = this->trgDOF();
00038   
00039   //1. zero out vecs.  This includes all global, contributor, user, evaluator vectors.
00040   PetscScalar zero=0.0;
00041   pC( VecSet(trgVal, zero) );
00042   pC( VecSet(_glbSrcExaDen, zero) );
00043   pC( VecSet(_glbSrcUpwEquDen, zero) );
00044   pC( VecSet(_ctbSrcExaDen, zero) );
00045   pC( VecSet(_ctbSrcUpwEquDen, zero) );
00046   pC( VecSet(_ctbSrcUpwChkVal, zero) );
00047   pC( VecSet(_usrSrcExaDen, zero) );
00048   pC( VecSet(_usrSrcUpwEquDen, zero) );
00049   pC( VecSet(_evaTrgExaVal, zero) );  
00050   pC( VecSet(_evaTrgDwnEquDen, zero) );
00051   pC( VecSet(_evaTrgDwnChkVal, zero) );
00052   
00053   clock_t ck0, ck1;
00054   vector<int> ordVec;
00055   pC( _let->upwOrderCollect(ordVec) ); //BOTTOM UP collection of nodes
00056 
00057   //2. for contributors, load exact densities
00058   int procLclStart, procLclEnd; _let->procLclRan(_srcPos, procLclStart, procLclEnd);
00059   double* darr; pC( VecGetArray(srcDen, &darr) );
00060   for(int i=0; i<ordVec.size(); i++) {
00061          int gNodeIdx = ordVec[i];
00062          if(_let->node(gNodeIdx).tag() & LET_CBTRNODE) {
00063                 if(_let->terminal(gNodeIdx)==true) {
00064                   DblNumVec ctbSrcExaDen(this->ctbSrcExaDen(gNodeIdx));
00065                   vector<int>& curVecIdxs = _let->node(gNodeIdx).ctbSrcOwnVecIdxs();
00066                   for(int k=0; k<curVecIdxs.size(); k++) {
00067                          int poff = curVecIdxs[k] - procLclStart;
00068                          for(int d=0; d<srcDOF; d++) {
00069                                 ctbSrcExaDen(k*srcDOF+d) = darr[poff*srcDOF+d];
00070                          }
00071                   }
00072                 }
00073          }
00074   }
00075   pC( VecRestoreArray(srcDen, &darr) );
00076     
00077   //SCATTER
00078   pC( VecScatterBegin( _ctbSrcExaDen, _glbSrcExaDen,    ADD_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaDen) );
00079   
00080   //3. up computation
00081     for(int i=0; i<ordVec.size(); i++) {
00082          int gNodeIdx = ordVec[i];
00083          if( _let->node(gNodeIdx).tag() & LET_CBTRNODE) {
00084                 if(_let->depth(gNodeIdx)>=0) {
00085                   DblNumVec ctbSrcUpwChkValgNodeIdx(ctbSrcUpwChkVal(gNodeIdx));
00086                   DblNumVec ctbSrcUpwEquDengNodeIdx(ctbSrcUpwEquDen(gNodeIdx));
00087                   if(_let->terminal(gNodeIdx)==true) {
00088                          //S2M
00089                          pC( SrcEqu2UpwChk_dgemv(ctbSrcExaPos(gNodeIdx), ctbSrcExaNor(gNodeIdx), _let->center(gNodeIdx), _let->radius(gNodeIdx), ctbSrcExaDen(gNodeIdx), ctbSrcUpwChkValgNodeIdx) );
00090                   } else {
00091                          //M2M
00092                          for(int a=0; a<2; a++) for(int b=0; b<2; b++) for(int c=0; c<2; c++) {
00093                                 Index3 idx(a,b,c);
00094                                 int chi = _let->child(gNodeIdx, idx);
00095                                 if(_let->node(chi).tag() & LET_CBTRNODE) {
00096                                   pC( _matmgnt->UpwEqu2UpwChk_dgemv(_let->depth(chi)+_rootLevel, idx, ctbSrcUpwEquDen(chi), ctbSrcUpwChkValgNodeIdx) );
00097                                 }
00098                          }
00099                   }
00100                   //M2M
00101                   pC( _matmgnt->UpwChk2UpwEqu_dgemv(_let->depth(gNodeIdx)+_rootLevel, ctbSrcUpwChkValgNodeIdx, ctbSrcUpwEquDengNodeIdx) );
00102                 }
00103          }
00104   }
00105   
00106   //4. vectbscatters
00107   //SCATTER
00108   pC( VecScatterBegin( _ctbSrcUpwEquDen, _glbSrcUpwEquDen,    ADD_VALUES, SCATTER_FORWARD, _ctb2GlbSrcUpwEquDen) );
00109   pC( VecScatterEnd(   _ctbSrcExaDen, _glbSrcExaDen,    ADD_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaDen) );
00110   pC( VecScatterEnd(   _ctbSrcUpwEquDen, _glbSrcUpwEquDen,    ADD_VALUES, SCATTER_FORWARD, _ctb2GlbSrcUpwEquDen) );
00111   
00112   pC( VecScatterBegin( _glbSrcExaDen, _usrSrcExaDen, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaDen) );
00113   pC( VecScatterBegin( _glbSrcUpwEquDen, _usrSrcUpwEquDen, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcUpwEquDen) );
00114   pC( VecScatterEnd(   _glbSrcExaDen, _usrSrcExaDen, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaDen) );
00115   
00116   ordVec.clear();  pC( _let->dwnOrderCollect(ordVec) );
00117   //U
00118   for(int i=0; i<ordVec.size(); i++) {
00119          int gNodeIdx = ordVec[i];
00120          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00121                 if( _let->terminal(gNodeIdx)==true ) { //terminal
00122                   DblNumVec evaTrgExaValgNodeIdx(evaTrgExaVal(gNodeIdx));
00123                   DblNumMat evaTrgExaPosgNodeIdx(evaTrgExaPos(gNodeIdx));
00124                   for(vector<int>::iterator vi=_let->node(gNodeIdx).Unodes().begin(); vi!=_let->node(gNodeIdx).Unodes().end(); vi++) {
00125                          //S2T
00126                          pC( SrcEqu2TrgChk_dgemv(usrSrcExaPos(*vi), usrSrcExaNor(*vi), evaTrgExaPosgNodeIdx, usrSrcExaDen(*vi), evaTrgExaValgNodeIdx) );
00127                   }
00128                 }
00129          }
00130   }
00131   pC( VecScatterEnd(   _glbSrcUpwEquDen, _usrSrcUpwEquDen, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcUpwEquDen) );
00132   
00133   //V
00134   for(int i=0; i<ordVec.size(); i++) {
00135          int gNodeIdx = ordVec[i];
00136          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE ) { //evaluator          
00137                 Point3 gNodeIdxctr(_let->center(gNodeIdx));
00138                 double D = 2.0 * _let->radius(gNodeIdx);
00139                 
00140                 DblNumVec evaTrgDwnChkVal(this->evaTrgDwnChkVal(gNodeIdx));
00141                 for(vector<int>::iterator vi=_let->node(gNodeIdx).Vnodes().begin(); vi!=_let->node(gNodeIdx).Vnodes().end(); vi++) {
00142                   Point3 victr(_let->center(*vi));
00143                   Index3 idx;             for(int d=0; d<dim(); d++)                     idx(d) = int(round( (victr[d]-gNodeIdxctr[d])/D ));
00144                   
00145                   Node& srcnode = node(*vi);
00146                   Node& trgnode = node(gNodeIdx);
00147                   if(srcnode.vLstOthCnt()==0) {
00148                          srcnode.effDen().resize( _matmgnt->effDatSze(UE) );                     setvalue(srcnode.effDen(), 0.0);//1. resize effDen
00149                          pC( _matmgnt->plnDen2EffDen(_let->depth(gNodeIdx)+_rootLevel, usrSrcUpwEquDen(*vi),  srcnode.effDen()) );                       //2. transform from UpwEquDen to effDen
00150                   }
00151                   if(trgnode.vLstInCnt()==0) {
00152                          trgnode.effVal().resize( _matmgnt->effDatSze(DC) );                     setvalue(trgnode.effVal(), 0.0);                        //1. resize effVal
00153                   }
00154                   //M2L           
00155                   pC( _matmgnt->UpwEqu2DwnChk_dgemv(_let->depth(gNodeIdx)+_rootLevel, idx, srcnode.effDen(), trgnode.effVal()) );
00156                   
00157                   srcnode.vLstOthCnt()++;
00158                   trgnode.vLstInCnt()++;
00159                   if(srcnode.vLstOthCnt()==srcnode.vLstOthNum()) {
00160                          srcnode.effDen().resize(0);                     //1. resize effDen to 0
00161                          srcnode.vLstOthCnt()=0;
00162                   }
00163                   if(trgnode.vLstInCnt()==trgnode.vLstInNum()) {
00164                          pC( _matmgnt->effVal2PlnVal(_let->depth(gNodeIdx)+_rootLevel, trgnode.effVal(), evaTrgDwnChkVal) );                     //1. transform from effval to DwnChkVal
00165                          trgnode.effVal().resize(0); //2. resize effVal to 0
00166                          trgnode.vLstInCnt()=0;
00167                   }
00168                 }
00169          }
00170   }
00171 
00172   //W
00173   for(int i=0; i<ordVec.size(); i++) {
00174          int gNodeIdx = ordVec[i];
00175          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00176                 if( _let->terminal(gNodeIdx)==true ) {
00177                   DblNumVec evaTrgExaVal_gNodeIdx(this->evaTrgExaVal(gNodeIdx));
00178                   for(vector<int>::iterator vi=_let->node(gNodeIdx).Wnodes().begin(); vi!=_let->node(gNodeIdx).Wnodes().end(); vi++) {
00179                          if(_let->terminal(*vi) && _let->node(*vi).usrSrcExaNum()*srcDOF<_matmgnt->plnDatSze(UE)) { //use Exa instead
00180                                 //S2T
00181                                 pC( SrcEqu2TrgChk_dgemv(usrSrcExaPos(*vi), usrSrcExaNor(*vi), evaTrgExaPos(gNodeIdx), usrSrcExaDen(*vi), evaTrgExaVal_gNodeIdx) );
00182                          } else {
00183                                 //M2T
00184                                 int vni = *vi;          
00185                                 pC( UpwEqu2TrgChk_dgemv(_let->center(vni), _let->radius(vni), evaTrgExaPos(gNodeIdx), usrSrcUpwEquDen(*vi), evaTrgExaVal_gNodeIdx) );
00186                          }
00187                   }
00188                 }
00189          }
00190   }
00191 
00192   //X
00193   for(int i=0; i<ordVec.size(); i++) {
00194          int gNodeIdx = ordVec[i];
00195          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00196                 DblNumVec evaTrgExaVal_gNodeIdx(evaTrgExaVal(gNodeIdx));
00197                 DblNumVec evaTrgDwnChkVal_gNodeIdx(evaTrgDwnChkVal(gNodeIdx));
00198                 for(vector<int>::iterator vi=_let->node(gNodeIdx).Xnodes().begin(); vi!=_let->node(gNodeIdx).Xnodes().end(); vi++) {
00199                   if(_let->terminal(gNodeIdx) && _let->node(gNodeIdx).evaTrgExaNum()*trgDOF<_matmgnt->plnDatSze(DC)) { //use Exa instead
00200                          pC( SrcEqu2TrgChk_dgemv(usrSrcExaPos(*vi), usrSrcExaNor(*vi), evaTrgExaPos(gNodeIdx), usrSrcExaDen(*vi), evaTrgExaVal_gNodeIdx) );
00201                   } else {
00202                          //S2L
00203                          pC( SrcEqu2DwnChk_dgemv(usrSrcExaPos(*vi), usrSrcExaNor(*vi), _let->center(gNodeIdx), _let->radius(gNodeIdx), usrSrcExaDen(*vi), evaTrgDwnChkVal_gNodeIdx) );
00204                   }
00205                 }
00206          }
00207   }
00208 
00209   //7. combine
00210   for(int i=0; i<ordVec.size(); i++) {
00211          int gNodeIdx = ordVec[i];
00212          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE ) { //evaluator  
00213                 if(_let->depth(gNodeIdx)>=3) {
00214                   int pargNodeIdx = _let->parent(gNodeIdx);     
00215                   Index3 chdidx( _let->path2Node(gNodeIdx)-2 * _let->path2Node(pargNodeIdx) );
00216                   //L2L
00217                   DblNumVec evaTrgDwnChkVal_gNodeIdx(evaTrgDwnChkVal(gNodeIdx));
00218                   pC( _matmgnt->DwnEqu2DwnChk_dgemv(_let->depth(pargNodeIdx)+_rootLevel, chdidx, evaTrgDwnEquDen(pargNodeIdx), evaTrgDwnChkVal_gNodeIdx) );
00219                 }
00220                 if(_let->depth(gNodeIdx)>=2) {
00221                   //L2L
00222                   DblNumVec evaTrgDwnEquDen_gNodeIdx(evaTrgDwnEquDen(gNodeIdx));
00223                   pC( _matmgnt->DwnChk2DwnEqu_dgemv(_let->depth(gNodeIdx)+_rootLevel, evaTrgDwnChkVal(gNodeIdx), evaTrgDwnEquDen_gNodeIdx) );
00224                 }
00225                 if(_let->terminal(gNodeIdx)) {
00226                   //L2T
00227                   DblNumVec evaTrgExaVal_gNodeIdx(evaTrgExaVal(gNodeIdx));
00228                   pC( DwnEqu2TrgChk_dgemv(_let->center(gNodeIdx), _let->radius(gNodeIdx), evaTrgExaPos(gNodeIdx), evaTrgDwnEquDen(gNodeIdx), evaTrgExaVal_gNodeIdx) );
00229                 }
00230          }
00231   }
00232   
00233   //8. save tdtExaVal
00234   _let->procLclRan(_trgPos, procLclStart, procLclEnd);
00235   double* varr; pC( VecGetArray(trgVal, &varr) );
00236   for(int i=0; i<ordVec.size(); i++) {
00237          int gNodeIdx = ordVec[i];
00238          if( _let->node(gNodeIdx).tag() & LET_EVTRNODE ) {
00239                 if( _let->terminal(gNodeIdx)==true ) {
00240                   DblNumVec evaTrgExaVal(this->evaTrgExaVal(gNodeIdx));
00241                   vector<int>& curVecIdxs = _let->node(gNodeIdx).evaTrgOwnVecIdxs();
00242                   for(int k=0; k<curVecIdxs.size(); k++) {
00243                          int poff = curVecIdxs[k] - procLclStart;
00244                          for(int d=0; d<trgDOF; d++) {
00245                                 varr[poff*trgDOF+d] = evaTrgExaVal(k*trgDOF+d);
00246                          }
00247                   }
00248                 }
00249          }
00250   }
00251   pC( VecRestoreArray(trgVal, &varr) );
00252   
00253   //----------------
00254   pC( MPI_Barrier(mpiComm()) );  //check vLstInCnt, vLstOthCnt
00255   return(0);
00256 }
00257 
00258 
00259 

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