vecmatop.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 "blas.h"
00019 #include "lapack.h"
00020 #include "svdrep.hpp"
00021 #include "numvec.hpp"
00022 #include "nummat.hpp"
00023 #include "numtns.hpp"
00024 
00025 #include "vecmatop.hpp"
00026 
00027 using std::cerr;
00028 
00029 // ---------------------------------------------------------------------- 
00030 int dscal(double alpha, DblNumVec& X)
00031 {
00032   dscal( X.m(), alpha, X.data() );
00033   return 0;
00034 }
00035 // ---------------------------------------------------------------------- 
00036 int dscal(int n, double alpha, double* X)
00037 {
00038   int incx = 1;
00039   DSCAL(&n, &alpha, X, &incx);
00040   return 0;
00041 }
00042 // ---------------------------------------------------------------------- 
00043 int daxpy(double a, const DblNumVec& X, DblNumVec& Y)
00044 {
00045   assert( X.m() == Y.m() );
00046   daxpy(X.m(), a, X.data(), Y.data());
00047   return 0;
00048 }
00049 // ---------------------------------------------------------------------- 
00050 int daxpy(double a, const DblNumMat& X, DblNumMat& Y)
00051 {
00052   assert( X.m() == Y.m() );  assert( X.n() == Y.n() );
00053   iC( daxpy(X.m()*X.n(), a, X.data(), Y.data()) );
00054   return 0;
00055 }
00056 // ---------------------------------------------------------------------- 
00057 int daxpy(int n, double a, double* X, double* Y)
00058 {
00059   int incx = 1;  int incy = 1;
00060   DAXPY(&n, &a, X, &incx, Y, &incy);
00061   return 0;
00062 }
00063 //Y <- a M X + b Y
00064 // ---------------------------------------------------------------------- 
00065 int dgemm(double alpha, const DblNumMat& A, const DblNumMat& B, double beta, DblNumMat& C)
00066 {
00067   assert( A.m() == C.m() );  assert( A.n() == B.m() );  assert( B.n() == C.n() );
00068   iC( dgemm(C.m(), C.n(), A.n(), alpha, A.data(), B.data(), beta, C.data()) );
00069   return 0;
00070 }
00071 // ---------------------------------------------------------------------- 
00072 int dgemm(int m, int n, int k, double alpha, double* A, double* B, double beta, double* C)
00073 {
00074   char transa = 'N';
00075   char transb = 'N';
00076   assert(m!=0 && n!=0 && k!=0);
00077   DGEMM(&transa, &transb, &m, &n, &k,
00078                   &alpha, A, &m, B, &k, &beta, C, &m);
00079   return 0;
00080 }
00081 // ---------------------------------------------------------------------- 
00082 int dger(double alpha, const DblNumVec& X, const DblNumVec& Y, DblNumMat& A)
00083 {
00084   assert(X.m() == A.m());
00085   assert(Y.m() == A.n());
00086   iC( dger(A.m(), A.n(), alpha, X.data(), Y.data(), A.data()) );
00087   return 0;
00088 }
00089 // ---------------------------------------------------------------------- 
00090 int dger(int m, int n, double alpha, double* X, double* Y, double* A)
00091 {
00092   assert(m!=0 && n!=0);
00093   int incx = 1;  int incy = 1;
00094   DGER(&m, &n, &alpha, X, &incx, Y, &incy, A, &m);
00095   return 0;
00096 }
00097 //Y <- a M X + b Y
00098 // ---------------------------------------------------------------------- 
00099 int dgemv(double alpha, const DblNumMat& A, const DblNumVec& X, double beta, DblNumVec& Y)
00100 {
00101   assert(Y.m() == A.m());
00102   assert(A.n() == X.m());
00103   iC( dgemv(A.m(), A.n(), alpha, A.data(), X.data(), beta, Y.data()) );
00104   return 0;
00105 }
00106 // ---------------------------------------------------------------------- 
00107 int dgemv(int m, int n, double alpha, double* A, double* X, double beta, double* Y)
00108 {
00109   char trans = 'N';
00110   assert(m!=0 && n!=0);
00111   int incx = 1;
00112   int incy = 1;
00113   DGEMV(&trans, &m, &n, &alpha, A, &m, X, &incx, &beta, Y, &incy);
00114   return 0;
00115 }
00116 // ---------------------------------------------------------------------- 
00117 int tran(const DblNumMat& M, DblNumMat& R)
00118 {
00119   assert(R.m()==M.n() && R.n()==M.m());  //R.resize(M.n(), M.m());
00120   for(int i=0; i<M.m(); i++)
00121          for(int j=0; j<M.n(); j++)
00122                 R(j,i) = M(i,j);
00123   return 0;
00124 }
00125 // ----------------------------------------------------------------------
00126 int pinv(const DblNumMat& M, double epsilon, DblNumMat& R)
00127 {
00128   assert(M.m() == R.n());  assert(M.n() == R.m());
00129   SVDRep svd;
00130   iC( svd.construct(epsilon, M) );
00131   //invert Svd    //cerr<<svd.S()(0)<<" "<<svd.S()(svd.S().m()-1)<<" "<<svd.S()(0)/svd.S()(svd.S().m()-1)<<endl;
00132   double cutoff = epsilon * svd.S()(0);  //double cutoff = epsilon;
00133   for(int i=0; i<svd.S().m(); i++) {
00134     if( svd.S()(i) >= cutoff) {
00135       svd.S()(i) = 1.0/(svd.S()(i));
00136          } else {
00137                 assert(0);      //svd.S()(i) = 0.0;
00138          }
00139   }
00140   DblNumMat UT(svd.U().n(),  svd.U().m());
00141   DblNumMat V( svd.VT().n(), svd.VT().m());
00142   iC( tran(svd.U(), UT) );
00143   iC( tran(svd.VT(), V) );
00144   for(int i=0; i<V.m(); i++)
00145     for(int j=0; j<V.n(); j++) {
00146       V(i,j) = V(i,j) * svd.S()(j);
00147          }
00148   
00149   char transa = 'N';
00150   char transb = 'N';
00151   double alpha = 1.0;
00152   double beta = 0.0;
00153   int m = V.m();
00154   int n = UT.n();
00155   int k = V.n();
00156   DGEMM(&transa, &transb, &m, &n, &k, &alpha,
00157                   V.data(), &m, UT.data(), &k, 
00158                   &beta, R.data(), &m);  
00159   
00160   return 0;
00161 }
00162 
00163 // ---------------------------------------------------------------------- 
00164 int inv(const DblNumMat& M, DblNumMat& R) //Gaussian Elimination
00165 {
00166   //OR pinv(M, 0.0, R);
00167   assert(M.m()==M.n() && R.m()==R.n() && M.m()==R.m());
00168   memcpy(R.data(), M.data(), M.m()*M.n()*sizeof(double));
00169   int info;
00170   int m = M.m();
00171   int* ipiv = new int[m];
00172   DGETRF(&m, &m, R.data(), &m, ipiv, &info); assert(info==0);
00173   int lwork = m;
00174   double* work = new double[lwork];
00175   DGETRI(&m, R.data(), &m, ipiv, work, &lwork, &info);  assert(info==0);
00176   delete [] ipiv;
00177   delete [] work;
00178   return 0;
00179 }
00180 // -------------------------------------------------------------------------------------------------------------------------------------
00181 // ---------------------------------------------------------------------- 
00182 int spev1d(int evflag, int dmflag, int dof, double* data, int m, double e, int i, double u, double* res)
00183 {
00184   int is[4];
00185   if(dmflag==DMFLAG_PERIOD) {
00186          for(int k=0; k<4; k++) is[k]=(i+k-1 + m) % m;
00187   } else {
00188          assert(i>=1 && i<=m-3);
00189          for(int k=0; k<4; k++) is[k]=(i+k-1);
00190   }
00191   DblNumMat M(dof,m,false,data); //assert(M.n()==n && M.m()==res.m()); //double dof = M.m();  //int cnt = 0;
00192   //---------------------------
00193   if(evflag & EVFLAG_VL) {
00194          double scl = 1.0;
00195          double us[4]; iC( spcoef(EVFLAG_VL, u, us) );
00196          for(int d=0; d<dof; d++)               res[d] = 0;
00197          for(int d=0; d<dof; d++) 
00198                 for(int a=0; a<4; a++) {
00199                   res[d] += us[a] * M(d,is[a]);
00200                 }
00201          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00202          res+=dof;
00203   }
00204   //---------------------------
00205   if(evflag & EVFLAG_FD) {
00206          double scl = double(m) / e;
00207          double us[4]; iC( spcoef(EVFLAG_FD, u, us) );
00208          for(int d=0; d<dof; d++)               res[d] = 0;
00209          for(int d=0; d<dof; d++)
00210                 for(int a=0; a<4; a++) {
00211                   res[d] += us[a] * M(d,is[a]);
00212                 }
00213          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00214          res+=dof;
00215   }
00216   //---------------------------
00217   if(evflag & EVFLAG_SD) {
00218          double scl = double(m*m)/(e*e);
00219          double us[4]; iC( spcoef(EVFLAG_SD, u, us) );
00220          for(int d=0; d<dof; d++)               res[d] = 0;
00221          for(int d=0; d<dof; d++)
00222                 for(int a=0; a<4; a++) {
00223                   res[d] += us[a] * M(d,is[a]);
00224                 }
00225          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00226          res+=dof;
00227   }
00228   return 0;
00229 }
00230 // ---------------------------------------------------------------------- 
00231 int spev2d(int evflag, int dmflag, int dof, double* data, int* mn, double* ef, int* ij, double* uv, double* res)
00232 {
00233   int m = mn[0];  int n = mn[1];
00234   double e = ef[0];  double f = ef[1];
00235   int i = ij[0];  int j = ij[1];
00236   double u = uv[0];  double v = uv[1];
00237   
00238   int is[4]; int js[4];
00239   if(dmflag==DMFLAG_PERIOD) {
00240          for(int k=0; k<4; k++) is[k]=(i+k-1 + m) % m;
00241          for(int k=0; k<4; k++) js[k]=(j+k-1 + n) % n;
00242   } else {
00243          assert(i>=1 && i<=m-3);
00244          for(int k=0; k<4; k++) is[k]=(i+k-1);
00245          assert(j>=1 && j<=n-3);
00246          for(int k=0; k<4; k++) js[k]=(j+k-1);
00247   }
00248   DblNumMat M(dof,m*n,false,data);
00249   double scl;
00250   double us[4], vs[4];
00251   //---------------------------
00252   if(evflag & EVFLAG_VL) {
00253          scl = 1.0;
00254          iC( spcoef(EVFLAG_VL, u, us) );
00255          iC( spcoef(EVFLAG_VL, v, vs) );
00256          for(int d=0; d<dof; d++)               res[d] = 0;
00257          for(int a=0; a<4; a++)
00258                 for(int b=0; b<4; b++) {
00259                   double coef = us[a]*vs[b]; 
00260                   for(int d=0; d<dof; d++)
00261                          res[d] += coef * M(d, is[a]+js[b]*m);
00262                 }
00263          for(int d=0; d<dof; d++)               res[d] *= scl;
00264          res+=dof;
00265   }
00266   //---------------------------
00267   if(evflag & EVFLAG_FD) {
00268          scl = double(m)/e;
00269          iC( spcoef(EVFLAG_FD, u, us) );
00270          iC( spcoef(EVFLAG_VL, v, vs) );
00271          for(int d=0; d<dof; d++)               res[d] = 0;
00272          for(int a=0; a<4; a++)
00273                 for(int b=0; b<4; b++) {
00274                   double coef = us[a]*vs[b];
00275                   for(int d=0; d<dof; d++)
00276                          res[d] += coef * M(d, is[a]+js[b]*m);
00277                 }
00278          for(int d=0; d<dof; d++)               res[d] *= scl;
00279          res+=dof;
00280          //...
00281          scl = double(n)/f;
00282          iC( spcoef(EVFLAG_VL, u, us) );
00283          iC( spcoef(EVFLAG_FD, v, vs) );
00284          for(int d=0; d<dof; d++)               res[d] = 0;
00285          for(int a=0; a<4; a++)
00286                 for(int b=0; b<4; b++) {
00287                   double coef = us[a]*vs[b]; 
00288                   for(int d=0; d<dof; d++)
00289                          res[d] += coef * M(d, is[a]+js[b]*m);
00290                 }
00291          for(int d=0; d<dof; d++)               res[d] *= scl;
00292          res+=dof;
00293   }
00294   //---------------------------
00295   if(evflag & EVFLAG_SD) {
00296          scl = double(m*m)/(e*e);
00297          iC( spcoef(EVFLAG_SD, u, us) );
00298          iC( spcoef(EVFLAG_VL, v, vs) );
00299          for(int d=0; d<dof; d++)               res[d] = 0;
00300          for(int a=0; a<4; a++)
00301                 for(int b=0; b<4; b++) {
00302                   double coef = us[a]*vs[b]; 
00303                   for(int d=0; d<dof; d++)
00304                          res[d] += coef * M(d, is[a]+js[b]*m);
00305                 }
00306          for(int d=0; d<dof; d++)               res[d] *= scl;
00307          res+=dof;
00308          //...
00309          scl = double(m*n)/(e*f);
00310          iC( spcoef(EVFLAG_FD, u, us) );
00311          iC( spcoef(EVFLAG_FD, v, vs) );
00312          for(int d=0; d<dof; d++)               res[d] = 0;
00313          for(int a=0; a<4; a++)
00314                 for(int b=0; b<4; b++) {
00315                   double coef = us[a]*vs[b]; 
00316                   for(int d=0; d<dof; d++)
00317                          res[d] += coef * M(d, is[a]+js[b]*m);
00318                 }
00319          for(int d=0; d<dof; d++)               res[d] *= scl;
00320          res+=dof;
00321          //...
00322          scl = double(n*n)/(f*f);
00323          iC( spcoef(EVFLAG_VL, u, us) );
00324          iC( spcoef(EVFLAG_SD, v, vs) );
00325          for(int d=0; d<dof; d++)               res[d] = 0;
00326          for(int a=0; a<4; a++)
00327                 for(int b=0; b<4; b++) {
00328                   double coef = us[a]*vs[b]; 
00329                   for(int d=0; d<dof; d++)
00330                          res[d] += coef * M(d, is[a]+js[b]*m);
00331                 }
00332          for(int d=0; d<dof; d++)               res[d] *= scl;
00333          res+=dof;
00334   }
00335   return 0;
00336 }
00337 // ---------------------------------------------------------------------- 
00338 int spcoef(int evflag, double u, double* us)
00339 {
00340   double u1 = u;
00341   double u2 = u*u;
00342   double u3 = u*u*u;
00343   if(       evflag==EVFLAG_VL) {
00344          us[0] = (  1 - 3*u1 + 3*u2 -   u3)/6.0;
00345          us[1] = (  4        - 6*u2 + 3*u3)/6.0;
00346          us[2] = (  1 + 3*u1 + 3*u2 - 3*u3)/6.0;
00347          us[3] = (                  +   u3)/6.0;
00348   } else if(evflag==EVFLAG_FD) {
00349          us[0] = (- 3 + 6*u1 - 3*u2)/6.0;
00350          us[1] = (    -12*u1 + 9*u2)/6.0;
00351          us[2] = (  3 + 6*u1 - 9*u2)/6.0;
00352          us[3] = (             3*u2)/6.0;
00353   } else if(evflag==EVFLAG_SD) {
00354          us[0] = (  6 - 6*u1 ) / 6.0;
00355          us[1] = (-12 +18*u1 ) / 6.0;
00356          us[2] = (  6 -18*u1 ) / 6.0;
00357          us[3] = (      6*u1 ) / 6.0;    //assert(0); //TODO;
00358   }  
00359   return 0;
00360 }
00361 // -------------------------------------------------------------------------------------------------------------------------------------
00362 // ---------------------------------------------------------------------- 
00363 int lagev1d(int evflag, int dmflag, int dof, double* data, int m, double e, int i, double u, double* res)
00364 {
00365   int is[4];
00366   if(dmflag==DMFLAG_PERIOD) {
00367          for(int k=0; k<4; k++) is[k]=(i+k-1 + m) % m;
00368   } else {
00369          assert(i>=1 && i<=m-3);
00370          for(int k=0; k<4; k++) is[k]=(i+k-1);
00371   }
00372   DblNumMat M(dof,m,false,data); //assert(M.n()==n && M.m()==res.m()); //double dof = M.m();  //int cnt = 0;
00373   //---------------------------
00374   if(evflag & EVFLAG_VL) {
00375          double scl = 1.0;
00376          double us[4]; iC( lagcoef(EVFLAG_VL, u, us) );
00377          for(int d=0; d<dof; d++)               res[d] = 0;
00378          for(int d=0; d<dof; d++) 
00379                 for(int a=0; a<4; a++) {
00380                   res[d] += us[a] * M(d,is[a]);
00381                 }
00382          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00383          res+=dof; //cnt ++;
00384   }
00385   //---------------------------
00386   if(evflag & EVFLAG_FD) {
00387          double scl = double(m) / e;
00388          double us[4]; iC( lagcoef(EVFLAG_FD, u, us) );
00389          for(int d=0; d<dof; d++)               res[d] = 0;
00390          for(int d=0; d<dof; d++)
00391                 for(int a=0; a<4; a++) {
00392                   res[d] += us[a] * M(d,is[a]);
00393                 }
00394          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00395          res+=dof; //cnt ++;
00396   }
00397   //---------------------------
00398   if(evflag & EVFLAG_SD) {
00399          double scl = double(m*m)/(e*e);
00400          double us[4]; iC( lagcoef(EVFLAG_SD, u, us) );
00401          for(int d=0; d<dof; d++)               res[d] = 0;
00402          for(int d=0; d<dof; d++)
00403                 for(int a=0; a<4; a++) {
00404                   res[d] += us[a] * M(d,is[a]);
00405                 }
00406          for(int d=0; d<dof; d++)               res[d] *= scl; //scaling
00407          res+=dof; //cnt ++;
00408   }
00409   return 0;
00410 }
00411 // ---------------------------------------------------------------------- 
00412 int lagev2d(int evflag, int dmflag, int dof, double* data, int* mn, double* ef, int* ij, double* uv, double* res)
00413 {
00414   int m = mn[0];  int n = mn[1];
00415   double e = ef[0];  double f = ef[1];
00416   int i = ij[0];  int j = ij[1];
00417   double u = uv[0];  double v = uv[1];
00418   
00419   int is[4]; int js[4];
00420   if(dmflag==DMFLAG_PERIOD) {
00421          for(int k=0; k<4; k++) is[k]=(i+k-1 + m) % m;
00422          for(int k=0; k<4; k++) js[k]=(j+k-1 + n) % n;
00423   } else {
00424          assert(i>=1 && i<=m-3);
00425          for(int k=0; k<4; k++) is[k]=(i+k-1);
00426          assert(j>=1 && j<=n-3);
00427          for(int k=0; k<4; k++) js[k]=(j+k-1);
00428   }
00429   DblNumMat M(dof,m*n,false,data);
00430   double scl; 
00431   double us[4], vs[4];
00432   //---------------------------
00433   if(evflag & EVFLAG_VL) {
00434          scl = 1.0;
00435          iC( lagcoef(EVFLAG_VL, u, us) );
00436          iC( lagcoef(EVFLAG_VL, v, vs) );
00437          for(int d=0; d<dof; d++)               res[d] = 0;
00438          for(int a=0; a<4; a++)
00439                 for(int b=0; b<4; b++) {
00440                   double coef = us[a]*vs[b]; 
00441                   for(int d=0; d<dof; d++)
00442                          res[d] += coef * M(d, is[a]+js[b]*m);
00443                 }
00444          for(int d=0; d<dof; d++)               res[d] *= scl;
00445          res+=dof;
00446   }
00447   //---------------------------
00448   if(evflag & EVFLAG_FD) {
00449          scl = double(m)/e;
00450          iC( lagcoef(EVFLAG_FD, u, us) );
00451          iC( lagcoef(EVFLAG_VL, v, vs) );
00452          for(int d=0; d<dof; d++)               res[d] = 0;
00453          for(int a=0; a<4; a++)
00454                 for(int b=0; b<4; b++) {
00455                   double coef = us[a]*vs[b];
00456                   for(int d=0; d<dof; d++)
00457                          res[d] += coef * M(d, is[a]+js[b]*m);
00458                 }
00459          for(int d=0; d<dof; d++)               res[d] *= scl;
00460          res+=dof;
00461          //...
00462          scl = double(n)/f;
00463          iC( lagcoef(EVFLAG_VL, u, us) );
00464          iC( lagcoef(EVFLAG_FD, v, vs) );
00465          for(int d=0; d<dof; d++)               res[d] = 0;
00466          for(int a=0; a<4; a++)
00467                 for(int b=0; b<4; b++) {
00468                   double coef = us[a]*vs[b]; 
00469                   for(int d=0; d<dof; d++)
00470                          res[d] += coef * M(d, is[a]+js[b]*m);
00471                 }
00472          for(int d=0; d<dof; d++)               res[d] *= scl;
00473          res+=dof;
00474   }
00475   //---------------------------
00476   if(evflag & EVFLAG_SD) {
00477          scl = double(m*m)/(e*e);
00478          iC( lagcoef(EVFLAG_SD, u, us) );
00479          iC( lagcoef(EVFLAG_VL, v, vs) );
00480          for(int d=0; d<dof; d++)               res[d] = 0;
00481          for(int a=0; a<4; a++)
00482                 for(int b=0; b<4; b++) {
00483                   double coef = us[a]*vs[b]; 
00484                   for(int d=0; d<dof; d++)
00485                          res[d] += coef * M(d, is[a]+js[b]*m);
00486                 }
00487          for(int d=0; d<dof; d++)               res[d] *= scl;
00488          res+=dof;
00489          //...
00490          scl = double(m*n)/(e*f);
00491          iC( lagcoef(EVFLAG_FD, u, us) );
00492          iC( lagcoef(EVFLAG_FD, v, vs) );
00493          for(int d=0; d<dof; d++)               res[d] = 0;
00494          for(int a=0; a<4; a++)
00495                 for(int b=0; b<4; b++) {
00496                   double coef = us[a]*vs[b]; 
00497                   for(int d=0; d<dof; d++)
00498                          res[d] += coef * M(d, is[a]+js[b]*m);
00499                 }
00500          for(int d=0; d<dof; d++)               res[d] *= scl;
00501          res+=dof;
00502          //...
00503          scl = double(n*n)/(f*f);
00504          iC( lagcoef(EVFLAG_VL, u, us) );
00505          iC( lagcoef(EVFLAG_SD, v, vs) );
00506          for(int d=0; d<dof; d++)               res[d] = 0;
00507          for(int a=0; a<4; a++)
00508                 for(int b=0; b<4; b++) {
00509                   double coef = us[a]*vs[b]; 
00510                   for(int d=0; d<dof; d++)
00511                          res[d] += coef * M(d, is[a]+js[b]*m);
00512                 }
00513          for(int d=0; d<dof; d++)               res[d] *= scl;
00514          res+=dof;
00515   }
00516   return 0;
00517 }
00518 // ---------------------------------------------------------------------- 
00519 int lagev3d(int evflag, int dmflag, int dof, double* data, int* mno, double* efg, int* ijk, double* uvw, double* res)
00520 {
00521   assert(evflag==EVFLAG_VL);
00522   
00523   int m = mno[0];  int n = mno[1];  int o = mno[2];
00524   double e = efg[0];  double f = efg[1];  double g = efg[2];
00525   int i = ijk[0];  int j = ijk[1];  int k = ijk[2];
00526   double u = uvw[0];  double v = uvw[1];  double w = uvw[2];
00527   
00528   int is[4]; int js[4];  int ks[4];
00529   if(dmflag==DMFLAG_PERIOD) {
00530          for(int h=0; h<4; h++) is[h]=(i+h-1 + m) % m;
00531          for(int h=0; h<4; h++) js[h]=(j+h-1 + n) % n;
00532          for(int h=0; h<4; h++) ks[h]=(k+h-1 + o) % o;
00533   } else {
00534          assert(i>=1 && i<=m-3);
00535          for(int h=0; h<4; h++) is[h]=(i+h-1);
00536          assert(j>=1 && j<=n-3);
00537          for(int h=0; h<4; h++) js[h]=(j+h-1);
00538          assert(k>=1 && k<=o-3);
00539          for(int h=0; h<4; h++) ks[h]=(k+h-1);
00540   }
00541   DblNumMat M(dof,m*n*o,false,data);
00542   double scl; 
00543   double us[4], vs[4], ws[4];
00544   //---------------------------
00545   if(evflag & EVFLAG_VL) {
00546          scl = 1.0;
00547          iC( lagcoef(EVFLAG_VL, u, us) );
00548          iC( lagcoef(EVFLAG_VL, v, vs) );
00549          iC( lagcoef(EVFLAG_VL, w, ws) );
00550          for(int d=0; d<dof; d++)               res[d] = 0;
00551          for(int a=0; a<4; a++)
00552                 for(int b=0; b<4; b++)
00553                   for(int c=0; c<4; c++) {
00554                          double coef = us[a]*vs[b]*ws[c]; 
00555                          for(int d=0; d<dof; d++)
00556                                 res[d] += coef * M(d, is[a]+js[b]*m+ks[c]*m*n);
00557                   }
00558          for(int d=0; d<dof; d++)               res[d] *= scl;
00559          res+=dof;
00560   }
00561   return 0;
00562 }
00563 // ---------------------------------------------------------------------- 
00564 int lagcoef(int evflag, double u, double* us)
00565 {
00566   //double u1 = u;
00567   double u2 = u*u;
00568   double u3 = u*u*u;
00569   if(       evflag==EVFLAG_VL) {
00570          us[0] = (u3-3*u2+2*u)/(-6.0);
00571          us[1] = (u3-2*u2-u+2)/(2.0);
00572          us[2] = (u3-u2-2*u)  /(-2.0);
00573          us[3] = (u3-u)       /(6.0);
00574   } else if(evflag==EVFLAG_FD) {
00575          us[0] = (3*u2-6*u+2)/(-6.0);
00576          us[1] = (3*u2-4*u-1)/(2.0);
00577          us[2] = (3*u2-2*u-2)/(-2.0);
00578          us[3] = (3*u2-1)    /(6.0);
00579   } else if(evflag==EVFLAG_SD) {
00580          assert(0); //TODO
00581   }
00582   return 0;
00583 }

Generated on Sun Dec 4 18:13:13 2005 for common by  doxygen 1.4.5