00001
00002
00003
00004
00005 #ifndef PLAYA_LINEARCOMBINATIONIMPL_HPP
00006 #define PLAYA_LINEARCOMBINATIONIMPL_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaOut.hpp"
00010 #include "PlayaTabs.hpp"
00011 #include "PlayaLinearCombinationDecl.hpp"
00012 #include "PlayaLinearOperatorDecl.hpp"
00013 #include "Teuchos_ScalarTraits.hpp"
00014
00015
00016 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00017 #include "PlayaVectorImpl.hpp"
00018 #include "PlayaLinearOperatorImpl.hpp"
00019 #include "PlayaSimpleTransposedOpImpl.hpp"
00020 #endif
00021
00022
00023 namespace Playa
00024 {
00025
00026 using Playa::Out;
00027 using Playa::Tabs;
00028 using std::endl;
00029
00030
00031
00032
00033
00034
00035 template <class Scalar>
00036 template <int N>
00037 inline
00038 Vector<Scalar>::Vector(const LCN<Scalar, N>& x)
00039 : Playa::Handle<VectorBase<Scalar> >()
00040 {
00041 this->ptr() = x.eval().ptr();
00042 }
00043
00044
00045
00046
00047
00048 template <class Scalar> inline
00049 Vector<Scalar>& Vector<Scalar>::operator=(const LCN<Scalar, 1>& lc)
00050 {
00051 const Vector<Scalar>& other = lc.vec(0);
00052 const Scalar& alpha = lc.coeff(0);
00053 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00054 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00055
00056 if (this->ptr().get() == other.ptr().get())
00057 {
00058
00059
00060 if (alpha == zero) this->zero();
00061 else if (alpha != one) this->scale(alpha);
00062 }
00063 else if (this->ptr().get() != 0 && this->space() == other.space())
00064 {
00065
00066
00067 this->ptr()->update(alpha, other.ptr().get(), zero);
00068 }
00069 else
00070 {
00071
00072
00073 Vector<Scalar> cp = other.copy();
00074 this->ptr() = cp.ptr();
00075 this->scale(alpha);
00076
00077 }
00078 return *this;
00079 }
00080
00081
00082 template <class Scalar> inline
00083 Vector<Scalar>& Vector<Scalar>::operator=(const LCN<Scalar, 2>& lc)
00084 {
00085 const Vector<Scalar>& x = lc.vec(0);
00086 const Scalar& alpha = lc.coeff(0);
00087 const Vector<Scalar>& y = lc.vec(1);
00088 const Scalar& beta = lc.coeff(1);
00089 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00090
00091
00092 TEUCHOS_TEST_FOR_EXCEPTION(!y.space().isCompatible(x.space()),
00093 std::runtime_error,
00094 "Spaces x=" << x.space() << " and y="
00095 << y.space() << " are not compatible in operator=(a*x + b*y)");
00096
00097 if (this->ptr().get() != 0 && this->space() == x.space())
00098 {
00099
00100
00101
00102
00103 this->update(alpha, x, beta, y, zero);
00104 }
00105 else
00106 {
00107
00108
00109 Vector<Scalar> e = lc.eval();
00110 this->ptr() = e.ptr();
00111 }
00112
00113 return *this;
00114 }
00115
00116 template <class Scalar> inline
00117 Vector<Scalar>& Vector<Scalar>::operator=(const LCN<Scalar, 3>& lc)
00118 {
00119 const Vector<Scalar>& x = lc.vec(0);
00120 const Scalar& alpha = lc.coeff(0);
00121 const Vector<Scalar>& y = lc.vec(1);
00122 const Scalar& beta = lc.coeff(1);
00123 const Vector<Scalar>& z = lc.vec(2);
00124 const Scalar& gamma = lc.coeff(2);
00125 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00126
00127 TEUCHOS_TEST_FOR_EXCEPTION(!y.space().isCompatible(x.space()),
00128 std::runtime_error,
00129 "Spaces x=" << x.space() << " and y="
00130 << y.space() << " are not compatible in operator=(a*x + b*y + c*z)");
00131
00132 TEUCHOS_TEST_FOR_EXCEPTION(!z.space().isCompatible(x.space()),
00133 std::runtime_error,
00134 "Spaces x=" << x.space() << " and z="
00135 << z.space() << " are not compatible in operator=(a*x + b*y + c*z)");
00136
00137 if (this->ptr().get() != 0 && this->space() == x.space())
00138 {
00139
00140
00141
00142
00143 this->update(alpha, x, beta, y, gamma, z, zero);
00144 }
00145 else
00146 {
00147
00148
00149 Vector e = lc.eval();
00150 this->ptr() = e.ptr();
00151 }
00152
00153 return *this;
00154 }
00155
00156 template <class Scalar>
00157 template <int N> inline
00158 Vector<Scalar>& Vector<Scalar>::operator=(const LCN<Scalar, N>& lc)
00159 {
00160 Vector e = lc.eval();
00161 this->ptr() = e.ptr();
00162
00163 return *this;
00164 }
00165
00166
00167
00168
00169
00170
00171
00172
00173 template <class Scalar> inline
00174 Vector<Scalar>& Vector<Scalar>::operator+=(const LCN<Scalar, 1>& lc)
00175 {
00176 const Vector<Scalar>& other = lc.vec(0);
00177 const Scalar& alpha = lc.coeff(0);
00178 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00179
00180 TEUCHOS_TEST_FOR_EXCEPTION(!this->space().isCompatible(other.space()),
00181 std::runtime_error,
00182 "Spaces this=" << this->space() << " and other="
00183 << other.space() << " are not compatible in operator+=()");
00184
00185 this->ptr()->update(alpha, other.ptr().get(), one);
00186
00187 return *this;
00188 }
00189
00190
00191 template <class Scalar> inline
00192 Vector<Scalar>& Vector<Scalar>::operator+=(const LCN<Scalar, 2>& lc)
00193 {
00194 const Vector<Scalar>& x = lc.vec(0);
00195 const Scalar& alpha = lc.coeff(0);
00196 const Vector<Scalar>& y = lc.vec(1);
00197 const Scalar& beta = lc.coeff(1);
00198 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00199
00200 TEUCHOS_TEST_FOR_EXCEPTION(!this->space().isCompatible(x.space()),
00201 std::runtime_error,
00202 "Spaces this=" << this->space() << " and other="
00203 << x.space() << " are not compatible in operator+=()");
00204
00205 this->update(alpha, x, beta, y, one);
00206
00207 return *this;
00208 }
00209
00210
00211 template <class Scalar>
00212 template <int N> inline
00213 Vector<Scalar>& Vector<Scalar>::operator+=(const LCN<Scalar, N>& lc)
00214 {
00215 Vector<Scalar> e = lc.eval();
00216 *this += e;
00217 return *this;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 template <class Scalar> inline
00229 Vector<Scalar>& Vector<Scalar>::operator-=(const LCN<Scalar, 1>& lc)
00230 {
00231 const Vector<Scalar>& other = lc.vec(0);
00232 const Scalar& alpha = -lc.coeff(0);
00233 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00234
00235 TEUCHOS_TEST_FOR_EXCEPTION(!this->space().isCompatible(other.space()),
00236 std::runtime_error,
00237 "Spaces this=" << this->space() << " and other="
00238 << other.space() << " are not compatible in operator+=()");
00239
00240 this->ptr()->update(alpha, other.ptr().get(), one);
00241
00242 return *this;
00243 }
00244
00245
00246 template <class Scalar> inline
00247 Vector<Scalar>& Vector<Scalar>::operator-=(const LCN<Scalar, 2>& lc)
00248 {
00249
00250 const Vector<Scalar>& x = lc.vec(0);
00251 const Scalar& alpha = -lc.coeff(0);
00252 const Vector<Scalar>& y = lc.vec(1);
00253 const Scalar& beta = -lc.coeff(1);
00254 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00255
00256 TEUCHOS_TEST_FOR_EXCEPTION(!this->space().isCompatible(x.space()),
00257 std::runtime_error,
00258 "Spaces this=" << this->space() << " and other="
00259 << x.space() << " are not compatible in operator+=()");
00260
00261 this->update(alpha, x, beta, y, one);
00262
00263 return *this;
00264 }
00265
00266
00267 template <class Scalar>
00268 template <int N> inline
00269 Vector<Scalar>& Vector<Scalar>::operator-=(const LCN<Scalar, N>& lc)
00270 {
00271 Vector<Scalar> e = lc.eval();
00272 *this -= e;
00273
00274 return *this;
00275 }
00276
00277
00278
00279
00280
00281
00282
00283
00284 template <class Scalar>
00285 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00286 const Vector<Scalar>& x)
00287 {
00288 Vector<Scalar> rtn;
00289 A.apply(x, rtn);
00290 return rtn;
00291 }
00292
00293
00294
00295 template <class Scalar, int N>
00296 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00297 const LCN<Scalar, N>& x)
00298 {
00299 return A*x.eval();
00300 }
00301
00302
00303 template <class Scalar>
00304 Vector<Scalar> operator*(const LinearOperator<Scalar>& A,
00305 const LCN<Scalar, 1>& x)
00306 {
00307 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00308 Vector<Scalar> rtn = A*x.vec(0);
00309 if (x.coeff(0)!=one) rtn.scale(x.coeff(0));
00310 return rtn;
00311 }
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 template <class Scalar> inline
00322 LCN<Scalar, 1> operator*(const Scalar& alpha,
00323 const Vector<Scalar>& x)
00324 {
00325 return LCN<Scalar, 1>(alpha, x);
00326 }
00327
00328
00329 template <class Scalar> inline
00330 LCN<Scalar, 1> operator*(const Vector<Scalar>& x,
00331 const Scalar& alpha)
00332 {
00333 return LCN<Scalar, 1>(alpha, x);
00334 }
00335
00336
00337
00338 template <class Scalar> inline
00339 LCN<Scalar, 1> operator/(const Vector<Scalar>& x,
00340 const Scalar& alpha)
00341 {
00342 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00343 return LCN<Scalar, 1>(one/alpha, x);
00344 }
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 template <class Scalar, int N> inline
00356 LCN<Scalar, N> operator*(const LCN<Scalar, N>& lc, const Scalar& beta)
00357 {
00358 LCN<Scalar, N> rtn(lc);
00359 rtn.multiply(beta);
00360 return rtn;
00361 }
00362
00363
00364 template <class Scalar, int N> inline
00365 LCN<Scalar, N> operator*(const Scalar& beta, const LCN<Scalar, N>& lc)
00366 {
00367 LCN<Scalar, N> rtn(lc);
00368 rtn.multiply(beta);
00369 return rtn;
00370 }
00371
00372
00373 template <class Scalar, int N> inline
00374 LCN<Scalar, N> operator/(const LCN<Scalar, N>& lc, const Scalar& beta)
00375 {
00376 LCN<Scalar, N> rtn(lc);
00377 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00378 rtn.multiply(one/beta);
00379 return rtn;
00380 }
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392 template <class Scalar> inline
00393 LCN<Scalar, 2> operator+(const Vector<Scalar>& x,
00394 const Vector<Scalar>& y)
00395 {
00396 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00397 return LCN<Scalar, 2>(one, x, one, y);
00398 }
00399
00400
00401 template <class Scalar> inline
00402 LCN<Scalar, 2> operator-(const Vector<Scalar>& x,
00403 const Vector<Scalar>& y)
00404 {
00405 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00406 return LCN<Scalar, 2>(one, x, -one, y);
00407 }
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418 template <class Scalar, int N, int M> inline
00419 LCN<Scalar, N+M> operator+(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g)
00420 {
00421 LCN<Scalar, N+M> rtn;
00422 for (int i=0; i<N; i++) rtn.set(i, f.coeff(i), f.vec(i));
00423 for (int i=0; i<M; i++) rtn.set(i+N, g.coeff(i), g.vec(i));
00424 return rtn;
00425 }
00426
00427
00428 template <class Scalar, int N> inline
00429 LCN<Scalar, N+1> operator+(const LCN<Scalar, N>& f, const Vector<Scalar>& g)
00430 {
00431 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00432 LCN<Scalar, N+1> rtn;
00433 for (int i=0; i<N; i++) rtn.set(i, f.coeff(i), f.vec(i));
00434 rtn.set(N, one, g);
00435 return rtn;
00436 }
00437
00438
00439 template <class Scalar, int N> inline
00440 LCN<Scalar, N+1> operator+(const Vector<Scalar>& f, const LCN<Scalar, N>& g)
00441 {
00442 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00443 LCN<Scalar, N+1> rtn;
00444 rtn.set(0, one, f);
00445 for (int i=0; i<N; i++) rtn.set(i+1, g.coeff(i), g.vec(i));
00446
00447 return rtn;
00448 }
00449
00450
00451
00452 template <class Scalar> inline
00453 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x)
00454 {
00455 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00456 return LCN<Scalar, 2>(lc, one, x);
00457 }
00458
00459
00460
00461
00462 template <class Scalar> inline
00463 LCN<Scalar, 2> operator+(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc)
00464 {
00465 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00466 return LCN<Scalar, 2>(one, x, lc);
00467 }
00468
00469
00470
00471 template <class Scalar> inline
00472 LCN<Scalar, 2> operator+(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by)
00473 {
00474 return LCN<Scalar, 2>(ax, by);
00475 }
00476
00477
00478
00479
00480 template <class Scalar> inline
00481 LCN<Scalar, 3> operator+(const LCN<Scalar, 1>& ax, const LCN<Scalar, 2>& bycz)
00482 {
00483 return LCN<Scalar, 3>(ax, bycz);
00484 }
00485
00486
00487
00488 template <class Scalar> inline
00489 LCN<Scalar, 3> operator+(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz)
00490 {
00491 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00492 return LCN<Scalar, 3>(LCN<Scalar, 1>(one, x), bycz);
00493 }
00494
00495
00496
00497
00498 template <class Scalar> inline
00499 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, const LCN<Scalar, 1>& cz)
00500 {
00501 return LCN<Scalar, 3>(axby, cz);
00502 }
00503
00504
00505
00506
00507 template <class Scalar> inline
00508 LCN<Scalar, 3> operator+(const LCN<Scalar, 2>& axby, const Vector<Scalar>& z)
00509 {
00510 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00511 return LCN<Scalar, 3>(axby, LCN<Scalar,1>(one,z));
00512 }
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524 template <class Scalar, int N, int M> inline
00525 LCN<Scalar, N+M> operator-(const LCN<Scalar, N>& f, const LCN<Scalar, M>& g)
00526 {
00527 LCN<Scalar, N+M> rtn;
00528 for (int i=0; i<N; i++) rtn.set(i, f.coeff(i), f.vec(i));
00529 for (int i=0; i<M; i++) rtn.set(i+N, -g.coeff(i), g.vec(i));
00530 return rtn;
00531 }
00532
00533
00534 template <class Scalar, int N> inline
00535 LCN<Scalar, N+1> operator-(const LCN<Scalar, N>& f, const Vector<Scalar>& g)
00536 {
00537 LCN<Scalar, N+1> rtn;
00538 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00539 for (int i=0; i<N; i++) rtn.set(i, f.coeff(i), f.vec(i));
00540 rtn.set(N, -one, g);
00541 return rtn;
00542 }
00543
00544
00545 template <class Scalar, int N> inline
00546 LCN<Scalar, N+1> operator-(const Vector<Scalar>& f, const LCN<Scalar, N>& g)
00547 {
00548 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00549 LCN<Scalar, N+1> rtn;
00550 rtn.set(0, one, f);
00551 for (int i=0; i<N; i++) rtn.set(i+1, -g.coeff(i), g.vec(i));
00552
00553 return rtn;
00554 }
00555
00556
00557
00558 template <class Scalar> inline
00559 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& lc, const Vector<Scalar>& x)
00560 {
00561 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00562 return LCN<Scalar, 2>(lc, -one, x);
00563 }
00564
00565
00566
00567
00568 template <class Scalar> inline
00569 LCN<Scalar, 2> operator-(const Vector<Scalar>& x, const LCN<Scalar, 1>& lc)
00570 {
00571 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00572 return LCN<Scalar, 2>(one, x, -lc.coeff(0), lc.vec(0));
00573 }
00574
00575
00576
00577 template <class Scalar> inline
00578 LCN<Scalar, 2> operator-(const LCN<Scalar, 1>& ax, const LCN<Scalar, 1>& by)
00579 {
00580 return LCN<Scalar, 2>(ax, -by);
00581 }
00582
00583
00584
00585
00586 template <class Scalar> inline
00587 LCN<Scalar, 3> operator-(const LCN<Scalar, 1>& ax, const LCN<Scalar, 2>& bycz)
00588 {
00589 return LCN<Scalar, 3>(ax, -bycz);
00590 }
00591
00592
00593
00594 template <class Scalar> inline
00595 LCN<Scalar, 3> operator-(const Vector<Scalar>& x, const LCN<Scalar, 2>& bycz)
00596 {
00597 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00598 return LCN<Scalar, 3>(LCN<Scalar, 1>(one, x), -bycz);
00599 }
00600
00601
00602
00603
00604 template <class Scalar> inline
00605 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const LCN<Scalar, 1>& cz)
00606 {
00607 return LCN<Scalar, 3>(axby, -cz);
00608 }
00609
00610
00611
00612
00613 template <class Scalar> inline
00614 LCN<Scalar, 3> operator-(const LCN<Scalar, 2>& axby, const Vector<Scalar>& z)
00615 {
00616 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00617 return LCN<Scalar, 3>(axby, LCN<Scalar,1>(-one,z));
00618 }
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 template <class Scalar, int N> inline
00631 Scalar norm1(const LCN<Scalar, N>& lc)
00632 {
00633 return lc.norm1();
00634 }
00635
00636
00637 template <class Scalar, int N> inline
00638 Scalar norm2(const LCN<Scalar, N>& lc)
00639 {
00640 return lc.norm2();
00641 }
00642
00643
00644 template <class Scalar, int N> inline
00645 Scalar normInf(const LCN<Scalar, N>& lc)
00646 {
00647 return lc.normInf();
00648 }
00649
00650
00651 template <class Scalar, int N> inline
00652 Vector<Scalar> abs(const LCN<Scalar, N>& lc)
00653 {
00654 return lc.abs();
00655 }
00656
00657
00658 template <class Scalar, int N> inline
00659 Scalar min(const LCN<Scalar, N>& lc)
00660 {
00661 return lc.min();
00662 }
00663
00664
00665 template <class Scalar, int N> inline
00666 Scalar max(const LCN<Scalar, N>& lc)
00667 {
00668 return lc.max();
00669 }
00670
00671
00672 template <class Scalar, int N> inline
00673 Vector<Scalar> reciprocal(const LCN<Scalar, N>& lc)
00674 {
00675 return lc.reciprocal();
00676 }
00677
00678
00679
00680 }
00681
00682
00683
00684 #endif