PlayaLinearCombinationDecl.hpp

00001 /* @HEADER@ */
00002 //   
00003  /* @HEADER@ */
00004 
00005 #ifndef PLAYA_LINEARCOMBINATIONDECL_HPP
00006 #define PLAYA_LINEARCOMBINATIONDECL_HPP
00007 
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaVectorDecl.hpp"
00010 #include "Teuchos_ScalarTraits.hpp"
00011 
00012 
00013 
00014 namespace Playa
00015 {
00016 
00017 template <class Scalar> class LinearOperator;
00018 
00019 template <class Scalar, int N>
00020 class LCNBase
00021 {
00022 public:
00024   LCNBase(){}
00025 
00027   virtual Vector<Scalar> eval() const = 0 ;
00028   
00030   const Vector<Scalar>& vec(int i) const {return x_[i];}
00031 
00033   const Scalar& coeff(int i) const {return a_[i];}
00034 
00036   void multiply(const Scalar& beta) 
00037     {for (int i=0; i<N; i++) a_[i] *= beta;}
00038 
00040   int size() const {return N;}
00041 
00043   Vector<Scalar> dotStar(const Vector<Scalar>& other) const 
00044     {return eval().dotStar(other);}
00045 
00047   Vector<Scalar> dotSlash(const Vector<Scalar>& other) const 
00048     {return eval().dotSlash(other);}
00049 
00051   Vector<Scalar> reciprocal() const {return eval().reciprocal();}
00052 
00054   Vector<Scalar> abs() const {return eval().abs();} 
00055 
00057   Scalar norm2() const {return eval().norm2();}
00058 
00060   Scalar norm1() const {return eval().norm1();}
00061 
00063   Scalar normInf() const {return eval().normInf();}
00064 
00066   Scalar max() const {return eval().max();}
00067 
00069   Scalar min() const {return eval().min();}
00070 
00071 protected:
00072 
00073   Vector<Scalar> x_[N];
00074   Scalar a_[N];
00075 };
00076 
00077 template <class Scalar, int N>
00078 class LCN : public LCNBase<Scalar, N>
00079 {
00080 public:
00082   LCN(){}
00083   
00085   void set(int i, const Scalar& a, const Vector<Scalar>& x)
00086     {
00087       this->a_[i] = a;
00088       this->x_[i] = x;
00089     }
00090 
00092   LCN<Scalar, N> operator-() const 
00093     {
00094       LCN<Scalar, N> rtn=*this; 
00095       rtn.multiply(-1.0);
00096       return rtn;
00097     }
00098 
00100   Vector<Scalar> eval() const 
00101     {
00102       Vector<Scalar> rtn = this->x_[0].copy();
00103       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00104       for (int i=1; i<N; i++)
00105       {
00106         rtn += this->a_[i]*this->x_[i];
00107       }
00108       return rtn;
00109     }
00110 };
00111 
00112 
00113 //template<>
00114 template <class Scalar>
00115 class LCN<Scalar, 1> : public LCNBase<Scalar, 1>
00116 {
00117 public:
00119   LCN(const Scalar& a, const Vector<Scalar>& x) : LCNBase<Scalar, 1>()
00120     {
00121       this->x_[0] = x;
00122       this->a_[0] = a;
00123     }
00124     
00126   LCN<Scalar, 1> operator-() const 
00127     {
00128       LCN<Scalar, 1> rtn=*this; 
00129       rtn.multiply(-1.0);
00130       return rtn;
00131     }
00132 
00134   Vector<Scalar> eval() const 
00135     {
00136       Vector<Scalar> rtn = this->x_[0].copy();
00137       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00138       return rtn;
00139     }
00140 };
00141 
00142 
00143 //template<>
00144 template <class Scalar>
00145 class LCN<Scalar, 2> : public LCNBase<Scalar, 2>
00146 {
00147 public:
00149   LCN(const Scalar& a, const Vector<Scalar>& x,
00150     const Scalar& b, const Vector<Scalar>& y) : LCNBase<Scalar, 2>()
00151     {
00152       this->x_[0] = x;
00153       this->a_[0] = a;
00154       this->x_[1] = y;
00155       this->a_[1] = b;
00156     }
00157 
00159   LCN(const LCN<Scalar, 1>& ax,
00160     const Scalar& b, const Vector<Scalar>& y) : LCNBase<Scalar, 2>()
00161     {
00162       this->x_[0] = ax.vec(0);
00163       this->a_[0] = ax.coeff(0);
00164       this->x_[1] = y;
00165       this->a_[1] = b;
00166     }
00167 
00169   LCN(const LCN<Scalar, 1>& ax,
00170     const LCN<Scalar, 1>& by) : LCNBase<Scalar, 2>()
00171     {
00172       this->x_[0] = ax.vec(0);
00173       this->a_[0] = ax.coeff(0);
00174       this->x_[1] = by.vec(0);
00175       this->a_[1] = by.coeff(0);
00176     }
00177     
00178 
00180   LCN(const Scalar& a, const Vector<Scalar>& x,
00181     const LCN<Scalar, 1>& by) : LCNBase<Scalar, 2>()
00182     {
00183       this->x_[0] = x;
00184       this->a_[0] = a;
00185       this->x_[1] = by.vec(0);
00186       this->a_[1] = by.coeff(0);
00187     }
00188     
00190   LCN<Scalar, 2> operator-() const 
00191     {
00192       LCN<Scalar, 2> rtn=*this; 
00193       rtn.multiply(-1.0);
00194       return rtn;
00195     }
00196 
00198   Vector<Scalar> eval() const 
00199     {
00200       Vector<Scalar> rtn = this->x_[0].copy();
00201       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00202       rtn += this->a_[1]*this->x_[1];
00203       return rtn;
00204     }
00205 };
00206 
00207 
00208 //template<>
00209 template <class Scalar>
00210 class LCN<Scalar, 3> : public LCNBase<Scalar, 3>
00211 {
00212 public:
00214   LCN(const LCN<Scalar, 1>& ax,
00215     const LCN<Scalar, 2>& bycz)
00216     {
00217       this->x_[0] = ax.vec(0);
00218       this->a_[0] = ax.coeff(0);
00219       this->x_[1] = bycz.vec(0);
00220       this->a_[1] = bycz.coeff(0);
00221       this->x_[2] = bycz.vec(1);
00222       this->a_[2] = bycz.coeff(1);
00223     }
00224     
00225 
00227   LCN(const LCN<Scalar, 2>& axby,
00228     const LCN<Scalar, 1>& cz)
00229     {
00230       this->x_[0] = axby.vec(0);
00231       this->a_[0] = axby.coeff(0);
00232       this->x_[1] = axby.vec(1);
00233       this->a_[1] = axby.coeff(1);
00234       this->x_[2] = cz.vec(0);
00235       this->a_[2] = cz.coeff(0);
00236     }
00237 
00239   LCN<Scalar, 3> operator-() const 
00240     {
00241       LCN<Scalar, 3> rtn=*this; 
00242       rtn.multiply(-1.0);
00243       return rtn;
00244     }
00245 
00247   Vector<Scalar> eval() const 
00248     {
00249       Vector<Scalar> rtn = this->x_[0].copy();
00250       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00251       rtn += this->a_[1]*this->x_[1] + this->a_[2]*this->x_[2];
00252       return rtn;
00253     }
00254 };
00255 
00256 
00257 
00258 
00259 
00260 /*======================================================================
00261  *
00262  *    Write a LC description
00263  *
00264  *======================================================================*/
00265 
00266 template <class Scalar, int N> inline
00267 std::ostream& operator<<(std::ostream& os, const LCN<Scalar, N>& lc)
00268 {
00269   os << "{(size=" << lc.size() << ") ";
00270   for (int i=0; i<lc.size(); i++)
00271   {
00272     if (i != 0) os << ", ";
00273     os << "(" << lc.coeff(i) << ", " << lc.vec(i).description() << ")";
00274   }
00275   os << "}";
00276   return os;
00277 }
00278 
00279 /*======================================================================
00280  *
00281  *    operator times vector or LC
00282  *
00283  *======================================================================*/
00284 
00286 template <class Scalar> 
00287 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00288   const Vector<Scalar>& x);
00289 
00290 
00292 template <class Scalar, int N> 
00293 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00294   const LCN<Scalar, N>& x);
00295   
00297 template <class Scalar> 
00298 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00299   const LCN<Scalar, 1>& x);
00300 
00301 
00302 /*======================================================================
00303  *
00304  *    scalar times vector
00305  *
00306  *======================================================================*/
00307 
00309 template <class Scalar> 
00310 LCN<Scalar, 1> operator*(const Scalar& alpha, 
00311   const Vector<Scalar>& x);
00312 
00314 template <class Scalar> 
00315 LCN<Scalar, 1> operator*(const Vector<Scalar>& x, 
00316   const Scalar& alpha);
00317 
00318 
00320 template <class Scalar> 
00321 LCN<Scalar, 1> operator/(const Vector<Scalar>& x, 
00322   const Scalar& alpha);
00323 
00324 
00325 /*======================================================================
00326  *
00327  *    scalar times LC
00328  *
00329  *======================================================================*/
00330 
00331 
00333 template <class Scalar, int N> 
00334 LCN<Scalar, N> operator*(const LCN<Scalar, N>& lc, const Scalar& beta);
00335 
00337 template <class Scalar, int N> 
00338 LCN<Scalar, N> operator*(const Scalar& beta, const LCN<Scalar, N>& lc);
00339 
00341 template <class Scalar, int N> 
00342 LCN<Scalar, N> operator/(const LCN<Scalar, N>& lc, const Scalar& beta);
00343 
00344 
00345   
00346 /*======================================================================
00347  *
00348  *  vector plus/minus vector
00349  *
00350  *======================================================================*/
00351 
00352 
00354 template <class Scalar> 
00355 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, 
00356   const Vector<Scalar>& y);
00357 
00359 template <class Scalar> 
00360 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, 
00361   const Vector<Scalar>& y);
00362 
00363 /*======================================================================
00364  *
00365  *  LC plus LC
00366  *
00367  *======================================================================*/
00368 
00370 template <class Scalar, int N, int M> 
00371 LCN<Scalar, N+M> operator+(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g);
00372 
00373 
00375 template <class Scalar, int N> 
00376 LCN<Scalar, N+1> operator+(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00377 
00378 
00380 template <class Scalar, int N> 
00381 LCN<Scalar, N+1> operator+(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00382 
00383 
00385 template <class Scalar> 
00386 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00387 
00389 template <class Scalar> 
00390 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00391 
00392 
00394 template <class Scalar> 
00395 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00396 
00397 
00399 template <class Scalar> 
00400 LCN<Scalar, 3> operator+(const LCN<Scalar, 1>& ax, 
00401   const LCN<Scalar, 2>& bycz);
00402 
00404 template <class Scalar> 
00405 LCN<Scalar, 3> operator+(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00406 
00407 
00409 template <class Scalar> 
00410 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, 
00411   const LCN<Scalar, 1>& cz);
00412 
00413 
00415 template <class Scalar> 
00416 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, 
00417   const Vector<Scalar>& z);
00418 
00419 
00420 /*======================================================================
00421  *
00422  *  LC minus LC
00423  *
00424  *======================================================================*/
00425 
00426 
00428 template <class Scalar, int N, int M> 
00429 LCN<Scalar, N+M> operator-(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g);
00430 
00431 
00433 template <class Scalar, int N> 
00434 LCN<Scalar, N+1> operator-(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00435 
00437 template <class Scalar, int N> 
00438 LCN<Scalar, N+1> operator-(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00439 
00441 template <class Scalar> 
00442 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00443 
00445 template <class Scalar> 
00446 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00447 
00449 template <class Scalar> 
00450 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00451 
00452 
00454 template <class Scalar> 
00455 LCN<Scalar, 3> operator-(const LCN<Scalar, 1>& ax, 
00456   const LCN<Scalar, 2>& bycz);
00457 
00459 template <class Scalar> 
00460 LCN<Scalar, 3> operator-(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00461 
00463 template <class Scalar> 
00464 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const LCN<Scalar, 1>& cz);
00465 
00467 template <class Scalar> 
00468 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const Vector<Scalar>& z);
00469 
00470 
00471 /*======================================================================
00472  *
00473  *  Operations on LC
00474  *
00475  *======================================================================*/
00476 
00478 template <class Scalar, int N> 
00479 Scalar norm1(const LCN<Scalar, N>& lc) ;
00480 
00482 template <class Scalar, int N> 
00483 Scalar norm2(const LCN<Scalar, N>& lc) ;
00484 
00486 template <class Scalar, int N> 
00487 Scalar normInf(const LCN<Scalar, N>& lc) ;
00488 
00489 
00491 template <class Scalar, int N> 
00492 Vector<Scalar> abs(const LCN<Scalar, N>& lc) ;
00493 
00495 template <class Scalar, int N> 
00496 Scalar min(const LCN<Scalar, N>& lc) ;
00497 
00498 
00500 template <class Scalar, int N> 
00501 Scalar max(const LCN<Scalar, N>& lc) ;
00502 
00504 template <class Scalar, int N> 
00505 Vector<Scalar> reciprocal(const LCN<Scalar, N>& lc) ;
00506 
00507 
00508 
00509 
00510 }
00511 
00512 
00513 
00514 #endif