matmgnt3d.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 "matmgnt3d.hpp"
00019 #include "common/vecmatop.hpp"
00020 
00021 using std::cerr;
00022 using std::endl;
00023 
00024 // ---------------------------------------------------------------------- 
00025 double MatMgnt3d::_wsbuf[16384];
00026 
00027 // ---------------------------------------------------------------------- 
00028 MatMgnt3d::MatMgnt3d(): _np(6)
00029 {
00030 }
00031 // ---------------------------------------------------------------------- 
00032 MatMgnt3d::~MatMgnt3d()
00033 {
00034 }
00035 
00036 double MatMgnt3d::alt()
00037 {
00038   return pow(0.1, _np+1);
00039 }
00040 // ---------------------------------------------------------------------- 
00041 // ---------------------------------------------------------------------- 
00042 int MatMgnt3d::setup()
00043 {
00044   //--------------------------------------------------------
00045   _hom = _knl.homogeneous();
00046   if(_hom==true) {
00047          _knl.homogeneousDeg(_degVec); iA(_degVec.size()==srcDOF());
00048   }
00049   iC( samPosCal(_np,   1.0, _samPos[UE]) );
00050   iC( samPosCal(_np+2, 3.0, _samPos[UC]) );
00051   iC( samPosCal(_np,   3.0, _samPos[DE]) );
00052   iC( samPosCal(_np,   1.0, _samPos[DC]) );
00053   
00054   iC( regPosCal(_np,   1.0, _regPos    ) ); //only one regPos
00055   //--------------------------------------------------------
00056   return (0);
00057 }
00058 // ---------------------------------------------------------------------- 
00059 int MatMgnt3d::report()
00060 {
00061   cerr << "matrix map size"<<endl;
00062   cerr << _upwEqu2UpwChk.size() << " ";
00063   cerr << _upwChk2UpwEqu.size() << " ";
00064   cerr << _dwnChk2DwnEqu.size() << " ";
00065   cerr << _dwnEqu2DwnChk.size() << " ";
00066   cerr << _upwEqu2DwnChk.size() << endl;
00067   return (0);
00068 }
00069 
00070 // ---------------------------------------------------------------------- 
00071 int MatMgnt3d::plnDatSze(int tp)
00072 {
00073   if(tp==UE || tp==DE)
00074          return _samPos[tp].n()*srcDOF();
00075   else
00076          return _samPos[tp].n()*trgDOF();
00077 }
00078 // ---------------------------------------------------------------------- 
00079 int MatMgnt3d::effDatSze(int tp)
00080 {
00081   //int et = eq().et();
00082   int effNum = (2*_np+2)*(2*_np)*(2*_np);
00083   if(tp==UE || tp==DE)
00084          return effNum*srcDOF();
00085   else
00086          return effNum*trgDOF();
00087 }
00088 
00089 // ---------------------------------------------------------------------- 
00090 int MatMgnt3d::UpwChk2UpwEqu_dgemv(int l, const DblNumVec& chk, DblNumVec& den)
00091 {
00092   DblNumMat& _UC2UE = (_hom==true) ? _upwChk2UpwEqu[0] : _upwChk2UpwEqu[l];
00093   double R = (_hom==true) ? 1 : 1.0/pow(2.0,l);
00094   //---------compute matrix
00095   if(_UC2UE.m()==0) {    //cerr<<"UpwChk2UpwEqu compute"<<endl;
00096          //set matrix
00097          DblNumMat ud2c(plnDatSze(UC), plnDatSze(UE));
00098          DblNumMat chkPos(dim(),samPos(UC).n());         clear(chkPos);  iC( daxpy(R, samPos(UC), chkPos) ); //scale
00099          DblNumMat denPos(dim(),samPos(UE).n());         clear(denPos);  iC( daxpy(R, samPos(UE), denPos) ); //scale
00100          
00101          iC( _knl.kernel(denPos, denPos, chkPos, ud2c) );
00102          _UC2UE.resize(plnDatSze(UE), plnDatSze(UC));
00103          iC( pinv(ud2c, alt(), _UC2UE) );
00104   }
00105   //---------matvec
00106   if(_hom==true) {
00107          //matvec
00108          int srcDOF = this->srcDOF();
00109          DblNumVec tmpDen(srcDOF*samPos(UE).n(), false, _wsbuf);         clear(tmpDen);
00110          iC( dgemv(1.0, _UC2UE, chk, 1.0, tmpDen) );
00111          //scale
00112          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, - l*_degVec[s]);
00113          int cnt = 0;
00114          for(int i=0; i < samPos(UE).n(); i++)
00115                 for(int s=0; s < srcDOF; s++) {
00116                   den(cnt) = den(cnt) + tmpDen(cnt) * sclvec[s];
00117                   cnt++;
00118                 }
00119   } else {
00120          iC( dgemv(1.0, _UC2UE, chk, 1.0, den) );
00121   }
00122   return (0);
00123 }
00124 // ---------------------------------------------------------------------- 
00125 int MatMgnt3d::UpwEqu2UpwChk_dgemv(int l, Index3 idx, const DblNumVec& den, DblNumVec& chk)
00126 {
00127   NumTns<DblNumMat>& _UE2UC = (_hom==true) ? _upwEqu2UpwChk[0] : _upwEqu2UpwChk[l];
00128   double R = (_hom==true) ? 1 : 1.0/pow(2.0, l);  
00129   if(_UE2UC.m()==0)      _UE2UC.resize(2,2,2);
00130   DblNumMat& _UE2UCii = _UE2UC(idx(0), idx(1), idx(2));
00131   //---------compute matrix
00132   if(_UE2UCii.m()==0) {  //cerr<<"UpwEqu2UpwChk compute"<<endl;
00133          _UE2UCii.resize(plnDatSze(UC), plnDatSze(UE)); //_memused[1] += plnnum(UC)*dof()*plnnum(UE)*dof()*sizeof(double);
00134          DblNumMat chkPos(dim(),samPos(UC).n());         clear(chkPos);  iC( daxpy(2.0*R, samPos(UC), chkPos) ); //scale
00135          DblNumMat denPos(dim(),samPos(UE).n());         clear(denPos);  iC( daxpy(R, samPos(UE), denPos) ); //scale
00136          for(int i=0; i<dim(); i++) for(int j=0; j<samPos(UE).n(); j++) denPos(i,j) = denPos(i,j) + (2*idx(i)-1)*R;//shift
00137          
00138          iC( _knl.kernel(denPos, denPos, chkPos, _UE2UCii) );
00139   }
00140   //---------matvec
00141   if(_hom==true) {
00142          int srcDOF = this->srcDOF();
00143          DblNumVec tmpDen(srcDOF*samPos(UE).n(), false, _wsbuf);         clear(tmpDen);
00144          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, l*_degVec[s]);
00145          int cnt = 0;
00146          for(int i=0; i<samPos(UE).n(); i++)
00147                 for(int s=0; s<srcDOF; s++) {
00148                   tmpDen(cnt) = den(cnt) * sclvec[s];
00149                   cnt++;
00150                 }
00151          iC( dgemv(1.0, _UE2UCii, tmpDen, 1.0, chk) );
00152   } else {
00153          iC( dgemv(1.0, _UE2UCii, den, 1.0, chk) );
00154   }
00155   return (0);
00156 }
00157 
00158 // ---------------------------------------------------------------------- 
00159 int MatMgnt3d::DwnChk2DwnEqu_dgemv(int l, const DblNumVec& chk, DblNumVec& den)
00160 {
00161   DblNumMat& _DC2DE = (_hom==true) ? _dwnChk2DwnEqu[0]: _dwnChk2DwnEqu[l];
00162   double R = (_hom==true) ? 1 : 1.0/pow(2.0,l);
00163   //---------compute matrix
00164   if(_DC2DE.m()==0) {    //cerr<<"DwnChk2DwnEqu compute"<<endl;
00165          DblNumMat dd2c(plnDatSze(DC), plnDatSze(DE));
00166          DblNumMat chkPos(dim(),samPos(DC).n());                clear(chkPos);   iC( daxpy(R, samPos(DC), chkPos) ); //scale
00167          DblNumMat denPos(dim(),samPos(DE).n());                clear(denPos);   iC( daxpy(R, samPos(DE), denPos) ); //scale
00168          
00169          iC( _knl.kernel(denPos, denPos, chkPos, dd2c) );//matrix
00170          _DC2DE.resize(plnDatSze(DE), plnDatSze(DC)); //_memused[2] += plndnenum()*dof()*plndncnum()*dof()*sizeof(double);
00171          iC( pinv(dd2c, alt(), _DC2DE) );
00172   }
00173   //---------matvec
00174   if(_hom==true) {
00175          int srcDOF = this->srcDOF();
00176          DblNumVec tmpDen(srcDOF*samPos(DE).n(), false, _wsbuf);         clear(tmpDen);
00177          iC( dgemv(1.0, _DC2DE, chk, 1.0, tmpDen) );
00178          //scale
00179          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, - l*_degVec[s]);
00180          int cnt = 0;
00181          for(int i=0; i<samPos(DE).n(); i++)
00182                 for(int s=0; s<srcDOF; s++) {
00183                   den(cnt) = den(cnt) + tmpDen(cnt) * sclvec[s];
00184                   cnt++;
00185                 }
00186   } else {
00187          iC( dgemv(1.0, _DC2DE, chk, 1.0, den) );
00188   }
00189   return 0;
00190 }
00191 // ---------------------------------------------------------------------- 
00192 int MatMgnt3d::DwnEqu2DwnChk_dgemv(int l, Index3 idx, const DblNumVec& den, DblNumVec& chk)
00193 {
00194   NumTns<DblNumMat>& _DE2DC = (_hom==true) ? _dwnEqu2DwnChk[0] : _dwnEqu2DwnChk[l];
00195   double R = (_hom==true) ? 1 : 1.0/pow(2.0, l);  //OffTns<DblNumMat>& _DwnEqu2DwnChk = _de2DwnChk[l];  double R       = 1.0/pow(2.0, l);
00196   if(_DE2DC.m()==0)      _DE2DC.resize(2,2,2);
00197   DblNumMat& _DE2DCii = _DE2DC(idx[0], idx[1], idx[2]);
00198   
00199   //---------compute matrix
00200   if(_DE2DCii.m()==0) {  //cerr<<"DwnEqu2DwnChk compute"<<endl;
00201          _DE2DCii.resize(plnDatSze(DC), plnDatSze(DE)); //_memused[3] += plndncnum()*dof()*plndnenum()*dof()*sizeof(double);
00202          DblNumMat denPos(dim(),samPos(DE).n());                  clear(denPos);         iC( daxpy(R, samPos(DE), denPos) ); //scale
00203          DblNumMat chkPos(dim(),samPos(DC).n());                  clear(chkPos);         iC( daxpy(0.5*R, samPos(DC), chkPos) ); //scale
00204          for(int i=0; i<dim(); i++) for(int j=0; j<samPos(DC).n(); j++) chkPos(i,j) = chkPos(i,j) + (double(idx(i))-0.5)*R;
00205          
00206          iC( _knl.kernel(denPos, denPos, chkPos, _DE2DCii) );
00207   }
00208   //---------matvec
00209   if(_hom==true) {
00210          int srcDOF = this->srcDOF();
00211          DblNumVec tmpDen(srcDOF*samPos(DE).n(), false, _wsbuf);         clear(tmpDen);
00212          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, l*_degVec[s]);
00213          int cnt = 0;
00214          for(int i=0; i<samPos(DE).n(); i++)
00215                 for(int s=0; s<srcDOF; s++) {
00216                   tmpDen(cnt) = den(cnt) * sclvec[s];
00217                   cnt++;
00218                 }
00219          iC( dgemv(1.0, _DE2DCii, tmpDen, 1.0, chk) );
00220   } else {
00221          iC( dgemv(1.0, _DE2DCii, den, 1.0, chk) );
00222   }
00223   return (0);
00224 }
00225 // ---------------------------------------------------------------------- 
00226 int MatMgnt3d::plnDen2EffDen(int l, const DblNumVec& plnDen, DblNumVec& effDen)
00227 {
00228   DblNumVec regDen(regPos().n()*srcDOF()); clear(regDen);
00229   //iC( samDen2RegularDen(plnDen, regDen) );
00230   //rfftwnd_real_to_complex(_forplan, srcDOF(), regDen.data(), srcDOF(), 1, (fftw_complex*)(effDen.data()), srcDOF(), 1);
00231   if(_hom==true) {
00232          int srcDOF = this->srcDOF();
00233          DblNumVec tmpDen(srcDOF*samPos(UE).n(), false, _wsbuf);         clear(tmpDen);
00234          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, l*_degVec[s]);
00235          int cnt = 0;
00236          for(int i=0; i<samPos(UE).n(); i++)
00237                 for(int s=0; s<srcDOF; s++) {
00238                   tmpDen(cnt) = plnDen(cnt) * sclvec[s];
00239                   cnt++;
00240                 }
00241          iC( samDen2RegDen(tmpDen, regDen) );
00242   } else {
00243          iC( samDen2RegDen(plnDen, regDen) );
00244   }
00245   
00246   int nnn[3];  nnn[0] = 2*_np;  nnn[1] = 2*_np;  nnn[2] = 2*_np;
00247   fftw_plan forplan = fftw_plan_many_dft_r2c(3,nnn,srcDOF(), regDen.data(),NULL, srcDOF(),1, (fftw_complex*)(effDen.data()),NULL, srcDOF(),1, FFTW_ESTIMATE);
00248   fftw_execute(forplan);
00249   fftw_destroy_plan(forplan);
00250   //rfftwnd_real_to_complex(_forplan, srcDOF(), regDen.data(), srcDOF(), 1, (fftw_complex*)(effDen.data()), srcDOF(), 1);
00251   
00252   return (0);
00253 }
00254 // ---------------------------------------------------------------------- 
00255 int MatMgnt3d::samDen2RegDen(const DblNumVec& samDen, DblNumVec& regDen)
00256 {
00257   int np = _np;
00258   int rgnum = 2*np;
00259   int srcDOF = this->srcDOF();
00260   int cnt=0;
00261   //the order of iterating is the same as SampleGrid
00262   for(int i=0; i<np; i++)
00263          for(int j=0; j<np; j++)
00264                 for(int k=0; k<np; k++) {
00265                   if(i==0 || i==np-1 || j==0 || j==np-1 || k==0 || k==np-1) {
00266                          //the position is fortran style
00267                          int rgoff = (k+np/2)*rgnum*rgnum + (j+np/2)*rgnum + (i+np/2);
00268                          for(int f=0; f<srcDOF; f++) {
00269                                 regDen(srcDOF*rgoff + f) += samDen(srcDOF*cnt + f);
00270                          }
00271                          cnt++;
00272                   }
00273                 }  //iC( PetscLogFlops(np*np*np*dof) );
00274   return 0;
00275 }
00276 // ---------------------------------------------------------------------- 
00277 int MatMgnt3d::effVal2PlnVal(int level, const DblNumVec& effVal, DblNumVec& plnVal)
00278 {
00279   DblNumVec regVal(regPos().n()*trgDOF());
00280   int nnn[3];  nnn[0] = 2*_np;  nnn[1] = 2*_np;  nnn[2] = 2*_np;
00281   fftw_plan invplan = fftw_plan_many_dft_c2r(3,nnn,trgDOF(), (fftw_complex*)(effVal.data()),NULL, trgDOF(),1, regVal.data(),NULL, trgDOF(),1, FFTW_ESTIMATE);
00282   fftw_execute(invplan);
00283   fftw_destroy_plan(invplan);  //rfftwnd_complex_to_real(_invplan, trgDOF(), (fftw_complex*)(effVal.data()), trgDOF(), 1, regVal.data(), trgDOF(), 1);
00284   
00285   iC( regVal2SamVal(regVal, plnVal) );
00286   return (0);
00287 }
00288 // ---------------------------------------------------------------------- 
00289 int MatMgnt3d::regVal2SamVal(const DblNumVec& regVal, DblNumVec& samVal)
00290 {
00291   int np = _np;
00292   int rgnum = 2*np;
00293   int trgDOF = this->trgDOF();
00294   int cnt=0;
00295   //the order of iterating is the same as SampleGrid
00296   for(int i=0; i<np; i++)
00297          for(int j=0; j<np; j++)
00298                 for(int k=0; k<np; k++) {
00299                   if(i==0 || i==np-1 || j==0 || j==np-1 || k==0 || k==np-1) {
00300                          //the position is fortran style
00301                          int rgoff = (k+np/2)*rgnum*rgnum + (j+np/2)*rgnum + (i+np/2);
00302                          for(int f=0; f<trgDOF; f++) {
00303                                 samVal(trgDOF*cnt + f) += regVal(trgDOF*rgoff + f);
00304                          }
00305                          cnt++;
00306                   }
00307                 }
00308   return 0;
00309 }
00310 // ---------------------------------------------------------------------- 
00311 int MatMgnt3d::UpwEqu2DwnChk_dgemv(int l, Index3 idx, const DblNumVec& effDen, DblNumVec& effVal)
00312 {
00313   OffTns<DblNumMat>& _UpwEqu2DwnChk = (_hom==true) ? _upwEqu2DwnChk[0] : _upwEqu2DwnChk[l];
00314   double R = (_hom==true) ? 1.0 : 1.0/pow(2.0, l); //OffTns< DblNumMat >& _UpwEqu2DwnChk = _upwEqu2DwnChk[l];  double R       = 1.0/pow(2.0, l);
00315   if(_UpwEqu2DwnChk.m()==0)      _UpwEqu2DwnChk.resize(7,7,7,-3,-3,-3);
00316   //DblNumMat* _UpwEqu2DwnChkii = new DblNumMat;
00317   DblNumMat& _UpwEqu2DwnChkii = _UpwEqu2DwnChk(idx[0], idx[1], idx[2]);
00318   
00319   int effNum = (2*_np+2)*(2*_np)*(2*_np);
00320   int srcDOF = this->srcDOF();
00321   int trgDOF = this->trgDOF();
00322   
00323   //---------compute matrix
00324   if(_UpwEqu2DwnChkii.m()==0) { //compute it if necessary
00325          //-----------------------       //cerr<<"UpwEqu2DwnChk(FFT) compute"<<endl;     //COMPUTE FFT   //Index3 idx = ii.idx();        
00326          iA( idx.linfty()>1 );
00327          DblNumMat denPos(dim(),1);      for(int i=0; i<dim(); i++)             denPos(i,0) = double(idx(i))*2.0*R; //shift
00328          DblNumMat chkPos(dim(),regPos().n());   clear(chkPos);  iC( daxpy(R, regPos(), chkPos) );
00329          DblNumMat tt(regPos().n()*trgDOF, srcDOF);
00330          iC( _knl.kernel(denPos, denPos, chkPos, tt) );
00331          // move data to tmp
00332          DblNumMat tmp(trgDOF,regPos().n()*srcDOF);
00333          for(int k=0; k<regPos().n();k++) {
00334                 for(int i=0; i<trgDOF; i++)
00335                   for(int j=0; j<srcDOF; j++) {
00336                          tmp(i,j+k*srcDOF) = tt(i+k*trgDOF,j);
00337                   }
00338          }
00339          _UpwEqu2DwnChkii.resize(trgDOF*srcDOF, effNum); //_memused[4] += dof*dof*effNum(UE)*sizeof(double);
00340          //forward FFT from tmp to _UpwEqu2DwnChkii;
00341          
00342          int nnn[3];  nnn[0] = 2*_np;  nnn[1] = 2*_np;  nnn[2] = 2*_np;
00343          fftw_plan forplan = fftw_plan_many_dft_r2c(3,nnn,srcDOF*trgDOF, tmp.data(),NULL, srcDOF*trgDOF,1, (fftw_complex*)(_UpwEqu2DwnChkii.data()),NULL, srcDOF*trgDOF,1, FFTW_ESTIMATE);
00344          fftw_execute(forplan);
00345          fftw_destroy_plan(forplan);
00346 
00347          //rfftwnd_real_to_complex(_forplan, srcDOF*trgDOF, tmp.data(), srcDOF*trgDOF, 1, (fftw_complex*)(_UpwEqu2DwnChkii.data()), srcDOF*trgDOF, 1);
00348   }
00349   //---------matvec
00350   /*
00351   //scale
00352   DblNumVec tmpDen(effDen.m(), false, _wsbuf);  clear(tmpDen);
00353   fftw_complex* denptr = NULL;// = (fftw_complex*)(tmpDen.data());
00354   if(_hom==true) {
00355          vector<double> sclvec(srcDOF);  for(int s=0; s<srcDOF; s++)            sclvec[s] = pow(2.0, l*_degVec[s]);
00356          int cnt = 0;
00357          for(int i=0; i<effDen.m()/srcDOF; i++)
00358                 for(int s=0; s<srcDOF; s++) {
00359                   tmpDen(cnt) = effDen(cnt) * sclvec[s];                cnt++;
00360                 }
00361          denptr = (fftw_complex*)(tmpDen.data());
00362   } else {
00363          denptr = (fftw_complex*)(effDen.data());
00364   }
00365   */
00366   //fft mult
00367   double nrmfc = 1.0/double(regPos().n());
00368   fftw_complex* matPtr  = (fftw_complex*)(_UpwEqu2DwnChkii.data());
00369   fftw_complex* denPtr = (fftw_complex*)(effDen.data());
00370   fftw_complex* chkPtr   = (fftw_complex*)(effVal.data());
00371   int matStp  = srcDOF*trgDOF;
00372   int denStp = srcDOF;
00373   int chkStp   = trgDOF;
00374   
00375   double newalpha = nrmfc;
00376   for(int i=0; i<trgDOF; i++)
00377          for(int j=0; j<srcDOF; j++) {
00378                 int matOff = j*trgDOF + i;
00379                 int denOff = j;
00380                 int chkOff = i;
00381                 iC( cptwvv(effNum/2, newalpha, matPtr+matOff, matStp, denPtr+denOff, denStp, chkPtr+chkOff, chkStp) );
00382          }
00383   return (0);
00384 }
00385 
00386 // ---------------------------------------------------------------------- 
00387 int MatMgnt3d::localPos(int tp, Point3 center, double radius, DblNumMat& positions)
00388 {
00389   const DblNumMat& bas = samPos(tp);
00390   positions.resize(dim(), bas.n());
00391   for(int i=0; i<dim(); i++)
00392          for(int j=0; j<positions.n(); j++)
00393                 positions(i,j) = center(i) + radius * bas(i,j);
00394   return (0);
00395 }
00396 
00397 // ---------------------------------------------------------------------- 
00398 int MatMgnt3d::samPosCal(int np, double R, DblNumMat& positions)
00399 {
00400   int n = np*np*np - (np-2)*(np-2)*(np-2);
00401   positions.resize(dim(),n);
00402   double step = 2.0/(np-1);
00403   double init = -1.0;
00404   int cnt = 0;
00405   for(int i=0; i<np; i++)
00406          for(int j=0; j<np; j++)
00407                 for(int k=0; k<np; k++) {
00408                   if(i==0 || i==np-1 || j==0 || j==np-1 || k==0 || k==np-1) {
00409                          double x = init + i*step;
00410                          double y = init + j*step;
00411                          double z = init + k*step;
00412                          positions(0,cnt) = R*x;
00413                          positions(1,cnt) = R*y;
00414                          positions(2,cnt) = R*z;
00415                          cnt++;
00416                   }
00417                 }
00418   iA(cnt==n);
00419   return 0;
00420 }
00421 
00422 // ---------------------------------------------------------------------- 
00423 int MatMgnt3d::regPosCal(int np, double R, DblNumMat& positions)
00424 {
00425   int n = 2*np*2*np*2*np;
00426   positions.resize(dim(), n);
00427   double step = 2.0/(np-1);
00428   int cnt = 0;
00429   for(int k=0; k<2*np; k++)
00430          for(int j=0; j<2*np; j++)
00431                 for(int i=0; i<2*np; i++) {
00432                   int gi = (i<np) ? i : i-2*np;
00433                   int gj = (j<np) ? j : j-2*np;
00434                   int gk = (k<np) ? k : k-2*np;
00435                   positions(0, cnt) = R * gi*step;
00436                   positions(1, cnt) = R * gj*step;
00437                   positions(2, cnt) = R * gk*step;
00438                   cnt ++;
00439                 }
00440   iA(cnt==n);
00441   return 0;
00442 }
00443 
00444 // ---------------------------------------------------------------------- 
00445 inline int MatMgnt3d::cptwvv(int n, double alpha, fftw_complex* x, int incrementX, fftw_complex* y, int incrementY, fftw_complex* z, int incrementZ)
00446 {
00447   for(int i=0; i<n; i++) {
00448          (*z)[0] += alpha * ( (*x)[0] * (*y)[0] - (*x)[1] * (*y)[1]);
00449          (*z)[1] += alpha * ( (*x)[0] * (*y)[1] + (*x)[1] * (*y)[0]);
00450          x = x + incrementX;
00451          y = y + incrementY;
00452          z = z + incrementZ;
00453   }  //iC( PetscLogFlops( 10*n ) );
00454   return 0;
00455 }
00456 
00457 //-----------------------------------------------------------------------------
00458 //-----------------------------------------------------------------------------
00459 //-----------------------------------------------------------------------------
00460 //-----------------------------------------------------------------------------
00461 vector<MatMgnt3d> MatMgnt3d::_mmvec;
00462 
00463 MatMgnt3d* MatMgnt3d::getmmptr(Kernel3d knl, int np)
00464 {
00465   for(int i=0; i<_mmvec.size(); i++)
00466          if(_mmvec[i].knl()==knl && _mmvec[i].np()==np)
00467                 return &(_mmvec[i]);
00468   
00469   _mmvec.push_back( MatMgnt3d() );
00470   int last = _mmvec.size()-1;
00471   MatMgnt3d* tmp = &(_mmvec[last]); //get the last one
00472   tmp->knl() = knl;  tmp->np() = np;
00473   tmp->setup();
00474   return tmp;
00475 }

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