# PlayaLinearCombinationDecl.hpp

```00001 /* @HEADER@ */
00002 //
00004
00005 #ifndef PLAYA_LINEARCOMBINATIONDECL_HPP
00006 #define PLAYA_LINEARCOMBINATIONDECL_HPP
00007
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
```