Sundance 2.1

Go to the documentation of this file.
00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                 Playa: Programmable Linear Algebra
00005 //                 Copyright 2012 Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Kevin Long (kevin.long@ttu.edu)
00038 // 
00039 
00040 /* @HEADER@ */
00041 
00042 #ifndef PLAYA_LINEARCOMBINATIONDECL_HPP
00043 #define PLAYA_LINEARCOMBINATIONDECL_HPP
00044 
00045 #include "PlayaDefs.hpp"
00046 #include "PlayaVectorDecl.hpp"
00047 #include "Teuchos_ScalarTraits.hpp"
00048 
00049 
00050 
00051 namespace Playa
00052 {
00053 
00054 template <class Scalar> class LinearOperator;
00055 
00056 template <class Scalar, int N>
00057 class LCNBase
00058 {
00059 public:
00060   /** */
00061   LCNBase(){}
00062 
00063   /** */
00064   virtual Vector<Scalar> eval() const = 0 ;
00065   
00066   /** */
00067   const Vector<Scalar>& vec(int i) const {return x_[i];}
00068 
00069   /** */
00070   const Scalar& coeff(int i) const {return a_[i];}
00071 
00072   /** Multiply through by a scalar constant */
00073   void multiply(const Scalar& beta) 
00074     {for (int i=0; i<N; i++) a_[i] *= beta;}
00075 
00076   /** */
00077   int size() const {return N;}
00078 
00079   /** Element-by-element product (Matlab dot-star operator) */
00080   Vector<Scalar> dotStar(const Vector<Scalar>& other) const 
00081     {return eval().dotStar(other);}
00082 
00083   /** Element-by-element division (Matlab dot-slash operator) */
00084   Vector<Scalar> dotSlash(const Vector<Scalar>& other) const 
00085     {return eval().dotSlash(other);}
00086 
00087   /** Return element-by-element reciprocal as a new vector */
00088   Vector<Scalar> reciprocal() const {return eval().reciprocal();}
00089 
00090   /** Return element-by-element absolute value as a new vector */
00091   Vector<Scalar> abs() const {return eval().abs();} 
00092 
00093   /** Return Euclidean norm */
00094   Scalar norm2() const {return eval().norm2();}
00095 
00096   /** Return 1-norm */
00097   Scalar norm1() const {return eval().norm1();}
00098 
00099   /** Return infinity norm */
00100   Scalar normInf() const {return eval().normInf();}
00101 
00102   /** Return max element */
00103   Scalar max() const {return eval().max();}
00104 
00105   /** Return min element */
00106   Scalar min() const {return eval().min();}
00107 
00108 protected:
00109 
00110   Vector<Scalar> x_[N];
00111   Scalar a_[N];
00112 };
00113 
00114 template <class Scalar, int N>
00115 class LCN : public LCNBase<Scalar, N>
00116 {
00117 public:
00118   /** */
00119   LCN(){}
00120   
00121   /** */
00122   void set(int i, const Scalar& a, const Vector<Scalar>& x)
00123     {
00124       this->a_[i] = a;
00125       this->x_[i] = x;
00126     }
00127 
00128   /** Unary minus */
00129   LCN<Scalar, N> operator-() const 
00130     {
00131       LCN<Scalar, N> rtn=*this; 
00132       rtn.multiply(-1.0);
00133       return rtn;
00134     }
00135 
00136   /** */
00137   Vector<Scalar> eval() const 
00138     {
00139       Vector<Scalar> rtn = this->x_[0].copy();
00140       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00141       for (int i=1; i<N; i++)
00142       {
00143         rtn += this->a_[i]*this->x_[i];
00144       }
00145       return rtn;
00146     }
00147 };
00148 
00149 
00150 //template<>
00151 template <class Scalar>
00152 class LCN<Scalar, 1> : public LCNBase<Scalar, 1>
00153 {
00154 public:
00155   /** */
00156   LCN(const Scalar& a, const Vector<Scalar>& x) : LCNBase<Scalar, 1>()
00157     {
00158       this->x_[0] = x;
00159       this->a_[0] = a;
00160     }
00161     
00162   /** Unary minus */
00163   LCN<Scalar, 1> operator-() const 
00164     {
00165       LCN<Scalar, 1> rtn=*this; 
00166       rtn.multiply(-1.0);
00167       return rtn;
00168     }
00169 
00170   /** */
00171   Vector<Scalar> eval() const 
00172     {
00173       Vector<Scalar> rtn = this->x_[0].copy();
00174       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00175       return rtn;
00176     }
00177 };
00178 
00179 
00180 //template<>
00181 template <class Scalar>
00182 class LCN<Scalar, 2> : public LCNBase<Scalar, 2>
00183 {
00184 public:
00185   /** */
00186   LCN(const Scalar& a, const Vector<Scalar>& x,
00187     const Scalar& b, const Vector<Scalar>& y) : LCNBase<Scalar, 2>()
00188     {
00189       this->x_[0] = x;
00190       this->a_[0] = a;
00191       this->x_[1] = y;
00192       this->a_[1] = b;
00193     }
00194 
00195   /** */
00196   LCN(const LCN<Scalar, 1>& ax,
00197     const Scalar& b, const Vector<Scalar>& y) : LCNBase<Scalar, 2>()
00198     {
00199       this->x_[0] = ax.vec(0);
00200       this->a_[0] = ax.coeff(0);
00201       this->x_[1] = y;
00202       this->a_[1] = b;
00203     }
00204 
00205   /** */
00206   LCN(const LCN<Scalar, 1>& ax,
00207     const LCN<Scalar, 1>& by) : LCNBase<Scalar, 2>()
00208     {
00209       this->x_[0] = ax.vec(0);
00210       this->a_[0] = ax.coeff(0);
00211       this->x_[1] = by.vec(0);
00212       this->a_[1] = by.coeff(0);
00213     }
00214     
00215 
00216   /** */
00217   LCN(const Scalar& a, const Vector<Scalar>& x,
00218     const LCN<Scalar, 1>& by) : LCNBase<Scalar, 2>()
00219     {
00220       this->x_[0] = x;
00221       this->a_[0] = a;
00222       this->x_[1] = by.vec(0);
00223       this->a_[1] = by.coeff(0);
00224     }
00225     
00226   /** Unary minus */
00227   LCN<Scalar, 2> operator-() const 
00228     {
00229       LCN<Scalar, 2> rtn=*this; 
00230       rtn.multiply(-1.0);
00231       return rtn;
00232     }
00233 
00234   /** */
00235   Vector<Scalar> eval() const 
00236     {
00237       Vector<Scalar> rtn = this->x_[0].copy();
00238       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00239       rtn += this->a_[1]*this->x_[1];
00240       return rtn;
00241     }
00242 };
00243 
00244 
00245 //template<>
00246 template <class Scalar>
00247 class LCN<Scalar, 3> : public LCNBase<Scalar, 3>
00248 {
00249 public:
00250   /** */
00251   LCN(const LCN<Scalar, 1>& ax,
00252     const LCN<Scalar, 2>& bycz)
00253     {
00254       this->x_[0] = ax.vec(0);
00255       this->a_[0] = ax.coeff(0);
00256       this->x_[1] = bycz.vec(0);
00257       this->a_[1] = bycz.coeff(0);
00258       this->x_[2] = bycz.vec(1);
00259       this->a_[2] = bycz.coeff(1);
00260     }
00261     
00262 
00263   /** */
00264   LCN(const LCN<Scalar, 2>& axby,
00265     const LCN<Scalar, 1>& cz)
00266     {
00267       this->x_[0] = axby.vec(0);
00268       this->a_[0] = axby.coeff(0);
00269       this->x_[1] = axby.vec(1);
00270       this->a_[1] = axby.coeff(1);
00271       this->x_[2] = cz.vec(0);
00272       this->a_[2] = cz.coeff(0);
00273     }
00274 
00275   /** Unary minus */
00276   LCN<Scalar, 3> operator-() const 
00277     {
00278       LCN<Scalar, 3> rtn=*this; 
00279       rtn.multiply(-1.0);
00280       return rtn;
00281     }
00282 
00283   /** */
00284   Vector<Scalar> eval() const 
00285     {
00286       Vector<Scalar> rtn = this->x_[0].copy();
00287       if (this->a_[0] != 1.0) rtn.scale(this->a_[0]);
00288       rtn += this->a_[1]*this->x_[1] + this->a_[2]*this->x_[2];
00289       return rtn;
00290     }
00291 };
00292 
00293 
00294 
00295 
00296 
00297 /*======================================================================
00298  *
00299  *    Write a LC description
00300  *
00301  *======================================================================*/
00302 
00303 template <class Scalar, int N> inline
00304 std::ostream& operator<<(std::ostream& os, const LCN<Scalar, N>& lc)
00305 {
00306   os << "{(size=" << lc.size() << ") ";
00307   for (int i=0; i<lc.size(); i++)
00308   {
00309     if (i != 0) os << ", ";
00310     os << "(" << lc.coeff(i) << ", " << lc.vec(i).description() << ")";
00311   }
00312   os << "}";
00313   return os;
00314 }
00315 
00316 /*======================================================================
00317  *
00318  *    operator times vector or LC
00319  *
00320  *======================================================================*/
00321 
00322 /** */
00323 template <class Scalar> 
00324 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00325   const Vector<Scalar>& x);
00326 
00327 
00328 /** */
00329 template <class Scalar, int N> 
00330 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00331   const LCN<Scalar, N>& x);
00332   
00333 /** */
00334 template <class Scalar> 
00335 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00336   const LCN<Scalar, 1>& x);
00337 
00338 
00339 /*======================================================================
00340  *
00341  *    scalar times vector
00342  *
00343  *======================================================================*/
00344 
00345 /** \relates LCN scalar * vec */
00346 template <class Scalar> 
00347 LCN<Scalar, 1> operator*(const Scalar& alpha, 
00348   const Vector<Scalar>& x);
00349 
00350 /** \relates LCN vec * scalar */
00351 template <class Scalar> 
00352 LCN<Scalar, 1> operator*(const Vector<Scalar>& x, 
00353   const Scalar& alpha);
00354 
00355 
00356 /** \relates LCN vec / scalar */
00357 template <class Scalar> 
00358 LCN<Scalar, 1> operator/(const Vector<Scalar>& x, 
00359   const Scalar& alpha);
00360 
00361 
00362 /*======================================================================
00363  *
00364  *    scalar times LC
00365  *
00366  *======================================================================*/
00367 
00368 
00369 /** */
00370 template <class Scalar, int N> 
00371 LCN<Scalar, N> operator*(const LCN<Scalar, N>& lc, const Scalar& beta);
00372 
00373 /** */
00374 template <class Scalar, int N> 
00375 LCN<Scalar, N> operator*(const Scalar& beta, const LCN<Scalar, N>& lc);
00376 
00377 /** */
00378 template <class Scalar, int N> 
00379 LCN<Scalar, N> operator/(const LCN<Scalar, N>& lc, const Scalar& beta);
00380 
00381 
00382   
00383 /*======================================================================
00384  *
00385  *  vector plus/minus vector
00386  *
00387  *======================================================================*/
00388 
00389 
00390 /** */
00391 template <class Scalar> 
00392 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, 
00393   const Vector<Scalar>& y);
00394 
00395 /** */
00396 template <class Scalar> 
00397 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, 
00398   const Vector<Scalar>& y);
00399 
00400 /*======================================================================
00401  *
00402  *  LC plus LC
00403  *
00404  *======================================================================*/
00405 
00406 /** */
00407 template <class Scalar, int N, int M> 
00408 LCN<Scalar, N+M> operator+(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g);
00409 
00410 
00411 /** */
00412 template <class Scalar, int N> 
00413 LCN<Scalar, N+1> operator+(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00414 
00415 
00416 /** */
00417 template <class Scalar, int N> 
00418 LCN<Scalar, N+1> operator+(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00419 
00420 
00421 /** */
00422 template <class Scalar> 
00423 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00424 
00425 /** */
00426 template <class Scalar> 
00427 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00428 
00429 
00430 /** */
00431 template <class Scalar> 
00432 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00433 
00434 
00435 /** */
00436 template <class Scalar> 
00437 LCN<Scalar, 3> operator+(const LCN<Scalar, 1>& ax, 
00438   const LCN<Scalar, 2>& bycz);
00439 
00440 /** */
00441 template <class Scalar> 
00442 LCN<Scalar, 3> operator+(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00443 
00444 
00445 /** */
00446 template <class Scalar> 
00447 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, 
00448   const LCN<Scalar, 1>& cz);
00449 
00450 
00451 /** */
00452 template <class Scalar> 
00453 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, 
00454   const Vector<Scalar>& z);
00455 
00456 
00457 /*======================================================================
00458  *
00459  *  LC minus LC
00460  *
00461  *======================================================================*/
00462 
00463 
00464 /** */
00465 template <class Scalar, int N, int M> 
00466 LCN<Scalar, N+M> operator-(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g);
00467 
00468 
00469 /** */
00470 template <class Scalar, int N> 
00471 LCN<Scalar, N+1> operator-(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00472 
00473 /** */
00474 template <class Scalar, int N> 
00475 LCN<Scalar, N+1> operator-(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00476 
00477 /** */
00478 template <class Scalar> 
00479 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00480 
00481 /** */
00482 template <class Scalar> 
00483 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00484 
00485 /** */
00486 template <class Scalar> 
00487 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00488 
00489 
00490 /** */
00491 template <class Scalar> 
00492 LCN<Scalar, 3> operator-(const LCN<Scalar, 1>& ax, 
00493   const LCN<Scalar, 2>& bycz);
00494 
00495 /** */
00496 template <class Scalar> 
00497 LCN<Scalar, 3> operator-(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00498 
00499 /** */
00500 template <class Scalar> 
00501 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const LCN<Scalar, 1>& cz);
00502 
00503 /** */
00504 template <class Scalar> 
00505 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const Vector<Scalar>& z);
00506 
00507 
00508 /*======================================================================
00509  *
00510  *  Operations on LC
00511  *
00512  *======================================================================*/
00513 
00514 /** */
00515 template <class Scalar, int N> 
00516 Scalar norm1(const LCN<Scalar, N>& lc) ;
00517 
00518 /** */
00519 template <class Scalar, int N> 
00520 Scalar norm2(const LCN<Scalar, N>& lc) ;
00521 
00522 /** */
00523 template <class Scalar, int N> 
00524 Scalar normInf(const LCN<Scalar, N>& lc) ;
00525 
00526 
00527 /** */
00528 template <class Scalar, int N> 
00529 Vector<Scalar> abs(const LCN<Scalar, N>& lc) ;
00530 
00531 /** */
00532 template <class Scalar, int N> 
00533 Scalar min(const LCN<Scalar, N>& lc) ;
00534 
00535 
00536 /** */
00537 template <class Scalar, int N> 
00538 Scalar max(const LCN<Scalar, N>& lc) ;
00539 
00540 /** */
00541 template <class Scalar, int N> 
00542 Vector<Scalar> reciprocal(const LCN<Scalar, N>& lc) ;
00543 
00544 
00545 
00546 
00547 }
00548 
00549 
00550 
00551 #endif