00001
00002
00003
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:
00023
00024 LCNBase(){}
00025
00026
00027 virtual Vector<Scalar> eval() const = 0 ;
00028
00029
00030 const Vector<Scalar>& vec(int i) const {return x_[i];}
00031
00032
00033 const Scalar& coeff(int i) const {return a_[i];}
00034
00035
00036 void multiply(const Scalar& beta)
00037 {for (int i=0; i<N; i++) a_[i] *= beta;}
00038
00039
00040 int size() const {return N;}
00041
00042
00043 Vector<Scalar> dotStar(const Vector<Scalar>& other) const
00044 {return eval().dotStar(other);}
00045
00046
00047 Vector<Scalar> dotSlash(const Vector<Scalar>& other) const
00048 {return eval().dotSlash(other);}
00049
00050
00051 Vector<Scalar> reciprocal() const {return eval().reciprocal();}
00052
00053
00054 Vector<Scalar> abs() const {return eval().abs();}
00055
00056
00057 Scalar norm2() const {return eval().norm2();}
00058
00059
00060 Scalar norm1() const {return eval().norm1();}
00061
00062
00063 Scalar normInf() const {return eval().normInf();}
00064
00065
00066 Scalar max() const {return eval().max();}
00067
00068
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:
00081
00082 LCN(){}
00083
00084
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
00091
00092 LCN<Scalar, N> operator-() const
00093 {
00094 LCN<Scalar, N> rtn=*this;
00095 rtn.multiply(-1.0);
00096 return rtn;
00097 }
00098
00099
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
00114 template <class Scalar>
00115 class LCN<Scalar, 1> : public LCNBase<Scalar, 1>
00116 {
00117 public:
00118
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
00125
00126 LCN<Scalar, 1> operator-() const
00127 {
00128 LCN<Scalar, 1> rtn=*this;
00129 rtn.multiply(-1.0);
00130 return rtn;
00131 }
00132
00133
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
00144 template <class Scalar>
00145 class LCN<Scalar, 2> : public LCNBase<Scalar, 2>
00146 {
00147 public:
00148
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
00158
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
00168
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
00179
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
00189
00190 LCN<Scalar, 2> operator-() const
00191 {
00192 LCN<Scalar, 2> rtn=*this;
00193 rtn.multiply(-1.0);
00194 return rtn;
00195 }
00196
00197
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
00209 template <class Scalar>
00210 class LCN<Scalar, 3> : public LCNBase<Scalar, 3>
00211 {
00212 public:
00213
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
00226
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
00238
00239 LCN<Scalar, 3> operator-() const
00240 {
00241 LCN<Scalar, 3> rtn=*this;
00242 rtn.multiply(-1.0);
00243 return rtn;
00244 }
00245
00246
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
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
00282
00283
00284
00285
00286 template <class Scalar>
00287 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00288 const Vector<Scalar>& x);
00289
00290
00291
00292 template <class Scalar, int N>
00293 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00294 const LCN<Scalar, N>& x);
00295
00296
00297 template <class Scalar>
00298 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00299 const LCN<Scalar, 1>& x);
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 template <class Scalar>
00310 LCN<Scalar, 1> operator*(const Scalar& alpha,
00311 const Vector<Scalar>& x);
00312
00313
00314 template <class Scalar>
00315 LCN<Scalar, 1> operator*(const Vector<Scalar>& x,
00316 const Scalar& alpha);
00317
00318
00319
00320 template <class Scalar>
00321 LCN<Scalar, 1> operator/(const Vector<Scalar>& x,
00322 const Scalar& alpha);
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 template <class Scalar, int N>
00334 LCN<Scalar, N> operator*(const LCN<Scalar, N>& lc, const Scalar& beta);
00335
00336
00337 template <class Scalar, int N>
00338 LCN<Scalar, N> operator*(const Scalar& beta, const LCN<Scalar, N>& lc);
00339
00340
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
00349
00350
00351
00352
00353
00354 template <class Scalar>
00355 LCN<Scalar, 2> operator+(const Vector<Scalar>& x,
00356 const Vector<Scalar>& y);
00357
00358
00359 template <class Scalar>
00360 LCN<Scalar, 2> operator-(const Vector<Scalar>& x,
00361 const Vector<Scalar>& y);
00362
00363
00364
00365
00366
00367
00368
00369
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
00374
00375 template <class Scalar, int N>
00376 LCN<Scalar, N+1> operator+(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00377
00378
00379
00380 template <class Scalar, int N>
00381 LCN<Scalar, N+1> operator+(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00382
00383
00384
00385 template <class Scalar>
00386 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00387
00388
00389 template <class Scalar>
00390 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00391
00392
00393
00394 template <class Scalar>
00395 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00396
00397
00398
00399 template <class Scalar>
00400 LCN<Scalar, 3> operator+(const LCN<Scalar, 1>& ax,
00401 const LCN<Scalar, 2>& bycz);
00402
00403
00404 template <class Scalar>
00405 LCN<Scalar, 3> operator+(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00406
00407
00408
00409 template <class Scalar>
00410 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby,
00411 const LCN<Scalar, 1>& cz);
00412
00413
00414
00415 template <class Scalar>
00416 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby,
00417 const Vector<Scalar>& z);
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
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
00432
00433 template <class Scalar, int N>
00434 LCN<Scalar, N+1> operator-(const LCN<Scalar, N>& f, const Vector<Scalar>& g);
00435
00436
00437 template <class Scalar, int N>
00438 LCN<Scalar, N+1> operator-(const Vector<Scalar>& f, const LCN<Scalar, N>& g);
00439
00440
00441 template <class Scalar>
00442 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x);
00443
00444
00445 template <class Scalar>
00446 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc);
00447
00448
00449 template <class Scalar>
00450 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by);
00451
00452
00453
00454 template <class Scalar>
00455 LCN<Scalar, 3> operator-(const LCN<Scalar, 1>& ax,
00456 const LCN<Scalar, 2>& bycz);
00457
00458
00459 template <class Scalar>
00460 LCN<Scalar, 3> operator-(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz);
00461
00462
00463 template <class Scalar>
00464 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const LCN<Scalar, 1>& cz);
00465
00466
00467 template <class Scalar>
00468 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const Vector<Scalar>& z);
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 template <class Scalar, int N>
00479 Scalar norm1(const LCN<Scalar, N>& lc) ;
00480
00481
00482 template <class Scalar, int N>
00483 Scalar norm2(const LCN<Scalar, N>& lc) ;
00484
00485
00486 template <class Scalar, int N>
00487 Scalar normInf(const LCN<Scalar, N>& lc) ;
00488
00489
00490
00491 template <class Scalar, int N>
00492 Vector<Scalar> abs(const LCN<Scalar, N>& lc) ;
00493
00494
00495 template <class Scalar, int N>
00496 Scalar min(const LCN<Scalar, N>& lc) ;
00497
00498
00499
00500 template <class Scalar, int N>
00501 Scalar max(const LCN<Scalar, N>& lc) ;
00502
00503
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