Sundance 2.1

Go to the documentation of this file.
00001 /* @HEADER@ */
00002 /* ***********************************************************************
00003 // 
00004 //           TSFExtended: Trilinos Solver Framework Extended
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // **********************************************************************/
00027  /* @HEADER@ */
00028 
00029 #ifndef TSFLINEARCOMBINATIONIMPL_HPP
00030 #define TSFLINEARCOMBINATIONIMPL_HPP
00031 
00032 #include "SundanceDefs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceTabs.hpp"
00035 #include "TSFLinearCombinationDecl.hpp"
00036 
00037 
00038 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00039 #include "TSFVectorImpl.hpp"
00040 #include "TSFLinearOperatorImpl.hpp"
00041 #endif
00042 
00043 
00044 namespace TSFExtendedOps
00045 {
00046 
00047 using Sundance::Out;
00048 using Sundance::Tabs;
00049 
00050 
00051 /* -------- methods of OpTimesLC ------ */
00052 
00053 template <class Scalar, class Node> inline
00054 OpTimesLC<Scalar, Node>::OpTimesLC(const Scalar& alpha, 
00055   const Node& x)
00056   : alpha_(alpha), op_(), x_(x) 
00057 {;}
00058 
00059 template <class Scalar, class Node> inline
00060 OpTimesLC<Scalar, Node>
00061 ::OpTimesLC(const Scalar& alpha,
00062   const TSFExtended::LinearOperator<Scalar>& op, 
00063   const Node& x)
00064   : alpha_(alpha), op_(op), x_(x) 
00065 {;}
00066   
00067   
00068 template <class Scalar, class Node> inline
00069 void OpTimesLC<Scalar, Node>::evalInto(TSFExtended::Vector<Scalar>& result) const
00070 {
00071   if (op_.ptr().get() != 0)
00072   {
00073     op_.apply(x_.eval(), result, alpha_);
00074   }
00075   else
00076   {
00077     x_.evalInto(result);
00078     if (alpha_ != one()) result.scale(alpha_);
00079   }
00080 }
00081 
00082 template <class Scalar, class Node> inline
00083 void OpTimesLC<Scalar, Node>::addInto(TSFExtended::Vector<Scalar>& result,
00084   LCSign sign) const
00085 {
00086   if (op_.ptr().get() != 0)
00087   {
00088     Vector<Scalar> tmp;
00089     op_.apply(x_.eval(), tmp);
00090     result.update(sign*alpha_, tmp);
00091   }
00092   else
00093   {
00094     result.update(sign*alpha_, x_.eval());
00095   }
00096 } 
00097 
00098 template <class Scalar, class Node> inline
00099 TSFExtended::Vector<Scalar> OpTimesLC<Scalar, Node>::eval() const 
00100 {
00101   TSFExtended::Vector<Scalar> result;
00102   if (op_.ptr().get() != 0)
00103   {
00104     result = op_.range().createMember();
00105     op_.apply(x_.eval(), result, alpha_);
00106   }
00107   else
00108   {
00109     result = x_.eval();
00110     if (alpha_ != one()) result.scale(alpha_);    
00111   }
00112   return result;
00113 }
00114 
00115 
00116 template <class Scalar, class Node> inline
00117 bool OpTimesLC<Scalar, Node>::containsVector(const Thyra::VectorBase<Scalar>* vec) const 
00118 {return x_.containsVector(vec);}
00119 
00120 
00121  
00122 
00123 
00124   
00125 /* ------------------------ methods of LC2 --------------------------- */
00126   
00127 template <class Scalar, class Node1, class Node2> inline
00128 LC2<Scalar, Node1, Node2>::LC2(const Node1& x1, const Node2& x2, LCSign sign)
00129   : x1_(x1), x2_(x2), sign_(sign) 
00130 {;}
00131 
00132 template <class Scalar, class Node1, class Node2> inline
00133 bool LC2<Scalar, Node1, Node2>::containsVector(const Thyra::VectorBase<Scalar>* vec) const
00134 {return x1_.containsVector(vec) || x2_.containsVector(vec);}
00135 
00136 template <class Scalar, class Node1, class Node2> inline
00137 void LC2<Scalar, Node1, Node2>::evalInto(TSFExtended::Vector<Scalar>& result) const
00138 {
00139   Tabs tab;
00140   x1_.evalInto(result);
00141   x2_.addInto(result, sign_);
00142 } 
00143 
00144 template <class Scalar, class Node1, class Node2> inline
00145 void LC2<Scalar, Node1, Node2>::addInto(TSFExtended::Vector<Scalar>& result,
00146   TSFExtendedOps::LCSign sign) const
00147 {
00148   x1_.addInto(result, sign);
00149   if (sign_*sign < 0) x2_.addInto(result, LCSubtract);
00150   else x2_.addInto(result, LCAdd);
00151 }
00152 
00153 template <class Scalar, class Node1, class Node2> inline
00154 TSFExtended::Vector<Scalar> LC2<Scalar, Node1, Node2>::eval() const
00155 {
00156   TSFExtended::Vector<Scalar> result = x1_.eval();
00157   x2_.addInto(result, sign_);
00158   return result;
00159 }
00160 }
00161 
00162 namespace TSFExtended
00163 {
00164 using TSFExtendedOps::OpTimesLC;
00165 using TSFExtendedOps::LC2;
00166 using TSFExtendedOps::LCAdd;
00167 using TSFExtendedOps::LCSubtract;
00168 
00169 /* ------------------------ global methods ----------------------- */
00170 
00171 
00172 /*======================================================================
00173  *
00174  *    scalar times vector
00175  *
00176  *======================================================================*/
00177 
00178 /* scalar * vec */
00179 template <class Scalar> inline
00180 OpTimesLC<Scalar, Vector<Scalar> > operator*(const Scalar& alpha, 
00181   const Vector<Scalar>& x)
00182 {
00183   return OpTimesLC<Scalar, Vector<Scalar> >(alpha, x);
00184 } 
00185 
00186 /* vec * scalar */
00187 template <class Scalar> inline
00188 OpTimesLC<Scalar, Vector<Scalar> > operator*(const Vector<Scalar>& x, 
00189   const Scalar& alpha)
00190 {
00191   return OpTimesLC<Scalar, Vector<Scalar> >(alpha, x);
00192 }
00193 
00194 
00195 /*======================================================================
00196  *
00197  *    scalar times OpTimesLC
00198  *
00199  *======================================================================*/
00200 
00201 /* scalar * OpTimesLC */
00202 template <class Scalar, class Node> inline
00203 OpTimesLC<Scalar, Node> 
00204 operator*(const Scalar& alpha, 
00205   const OpTimesLC<Scalar, Node>& x)
00206 {
00207   return OpTimesLC<Scalar, Node>(alpha * x.alpha(), x.op(), x.node());
00208 }
00209 
00210 /* OpTimesLC * scalar */
00211 template <class Scalar, class Node> inline
00212 OpTimesLC<Scalar, Node> 
00213 operator*(const OpTimesLC<Scalar, Node>& x, const Scalar& alpha)
00214 {
00215   return alpha * x;
00216 }
00217 
00218 
00219 /*======================================================================
00220  *
00221  *    scalar times LC2
00222  *
00223  *======================================================================*/
00224 
00225 /* scalar * LC2 */
00226 template <class Scalar, class Node1, class Node2> inline
00227 OpTimesLC<Scalar, LC2<Scalar, Node1, Node2> > 
00228 operator*(const Scalar& alpha, 
00229   const LC2<Scalar, Node1, Node2>& x)
00230 {
00231   return OpTimesLC<Scalar, LC2<Scalar, Node1, Node2> >(alpha, x);
00232 }
00233 
00234 /* LC2 * scalar */
00235 template <class Scalar, class Node1, class Node2> inline
00236 OpTimesLC<Scalar, LC2<Scalar, Node1, Node2> > 
00237 operator*(const LC2<Scalar, Node1, Node2>& x, const Scalar& alpha)
00238 {
00239   return alpha * x;
00240 }
00241   
00242 
00243 
00244 /*======================================================================
00245  *
00246  *    operator times [vectors, OpTimesLC, LC2]
00247  *
00248  *======================================================================*/
00249 
00250 /* op * vec */
00251 template <class Scalar> inline
00252 OpTimesLC<Scalar, Vector<Scalar> > 
00253 operator*(const LinearOperator<Scalar>& op, 
00254   const Vector<Scalar>& x)
00255 {
00256   return OpTimesLC<Scalar, Vector<Scalar> >(Teuchos::ScalarTraits<Scalar>::one(), op, x);
00257 }
00258 
00259 
00260 /* op * OpTimesLC */
00261 template <class Scalar, class Node> inline
00262 OpTimesLC<Scalar, Node> 
00263 operator*(const LinearOperator<Scalar>& op, 
00264   const OpTimesLC<Scalar, Node>& x)
00265 {
00266   TEST_FOR_EXCEPTION(op.ptr().get()==0, std::runtime_error,
00267     "null operator in LinearOperator * ( OpTimesLC )");
00268   if (x.op().ptr().get()==0)
00269   {
00270     return OpTimesLC<Scalar, Node>(x.alpha(), op, x.node());
00271   }
00272   else
00273   {
00274     return OpTimesLC<Scalar, Node>(x.alpha(), op * x.op(), x.node());
00275   }
00276 }
00277 
00278 
00279 /* op * LC2 */
00280 template <class Scalar, class Node1, class Node2> inline
00281 OpTimesLC<Scalar, LC2<Scalar, Node1, Node2> > 
00282 operator*(const LinearOperator<Scalar>& op, 
00283   const LC2<Scalar, Node1, Node2>& x)
00284 {
00285   return OpTimesLC<Scalar, LC2<Scalar, Node1, Node2> >(Teuchos::ScalarTraits<Scalar>::one(), op, x);
00286 }
00287 
00288 
00289 /*======================================================================
00290  *
00291  *    add/subtract vector, vector
00292  *
00293  *======================================================================*/
00294   
00295 /* vec + vec */
00296 template <class Scalar> inline
00297 LC2<Scalar, Vector<Scalar>, Vector<Scalar> >
00298 operator+(const Vector<Scalar>& x1, 
00299   const Vector<Scalar>& x2)
00300 {
00301   return LC2<Scalar, Vector<Scalar>, Vector<Scalar> >(x1, x2);
00302 }
00303   
00304 /* vec - vec */
00305 template <class Scalar> inline
00306 LC2<Scalar, Vector<Scalar>, Vector<Scalar> >
00307 operator-(const Vector<Scalar>& x1, 
00308   const Vector<Scalar>& x2)
00309 {
00310   return LC2<Scalar, Vector<Scalar>, Vector<Scalar> >(x1, x2, LCSubtract);
00311 }
00312 
00313 /*======================================================================
00314  *
00315  *    add/subtract vector, OpTimesLC
00316  *
00317  *======================================================================*/
00318 
00319 
00320 /* vec + OpTimesLC */
00321 template <class Scalar, class Node> inline
00322 LC2<Scalar, Vector<Scalar>, OpTimesLC<Scalar, Node> >
00323 operator+(const Vector<Scalar>& x1, 
00324   const OpTimesLC<Scalar, Node>& x2)
00325 {
00326   return LC2<Scalar, Vector<Scalar>, OpTimesLC<Scalar, Node> >(x1, x2);
00327 }
00328 
00329 
00330 /* vec - OpTimesLC */
00331 template <class Scalar, class Node> inline
00332 LC2<Scalar, Vector<Scalar>, OpTimesLC<Scalar, Node> >
00333 operator-(const Vector<Scalar>& x1, 
00334   const OpTimesLC<Scalar, Node>& x2)
00335 {
00336   return LC2<Scalar, Vector<Scalar>, OpTimesLC<Scalar, Node> >(x1, x2, 
00337     LCSubtract);
00338 }
00339 
00340 /* OpTimesLC + vec */
00341 template <class Scalar, class Node> inline
00342 LC2<Scalar, OpTimesLC<Scalar, Node>, Vector<Scalar> >
00343 operator+(const OpTimesLC<Scalar, Node>& x1, 
00344   const Vector<Scalar>& x2)
00345 {
00346   return LC2<Scalar, OpTimesLC<Scalar, Node>, Vector<Scalar> >(x1, x2);
00347 }
00348   
00349 /* OpTimesLC - vec */
00350 template <class Scalar, class Node> inline
00351 LC2<Scalar, OpTimesLC<Scalar, Node>, Vector<Scalar> >
00352 operator-(const OpTimesLC<Scalar, Node>& x1, 
00353   const Vector<Scalar>& x2)
00354 {
00355   return LC2<Scalar, OpTimesLC<Scalar, Node>, Vector<Scalar> >(x1, x2,
00356     LCSubtract);
00357 }
00358 
00359   
00360 /*======================================================================
00361  *
00362  *    add/subtract OpTimesLC, OpTimesLC
00363  *
00364  *======================================================================*/
00365   
00366 /* OpTimesLC + OpTimesLC */
00367 template <class Scalar, class Node1, class Node2> inline
00368 LC2<Scalar, OpTimesLC<Scalar, Node1>, OpTimesLC<Scalar, Node2> >
00369 operator+(const OpTimesLC<Scalar, Node1>& x1, 
00370   const OpTimesLC<Scalar, Node2>& x2)
00371 {
00372   return LC2<Scalar, OpTimesLC<Scalar, Node1>, 
00373     OpTimesLC<Scalar, Node2> >(x1, x2);
00374 }
00375   
00376 /* OpTimesLC - OpTimesLC */
00377 template <class Scalar, class Node1, class Node2> inline
00378 LC2<Scalar, OpTimesLC<Scalar, Node1>, OpTimesLC<Scalar, Node2> >
00379 operator-(const OpTimesLC<Scalar, Node1>& x1, 
00380   const OpTimesLC<Scalar, Node2>& x2)
00381 {
00382   return LC2<Scalar, OpTimesLC<Scalar, Node1>, 
00383     OpTimesLC<Scalar, Node2> >(x1, x2, LCSubtract);
00384 }
00385   
00386 
00387   
00388 /*======================================================================
00389  *
00390  *    add/subtract Vector, LC2
00391  *
00392  *======================================================================*/
00393 
00394   
00395 /* vec + LC2 */
00396 template <class Scalar, class Node1, class Node2> inline
00397 LC2<Scalar, Vector<Scalar>, LC2<Scalar, Node1, Node2> >
00398 operator+(const Vector<Scalar>& x1, 
00399   const LC2<Scalar, Node1, Node2>& x2)
00400 {
00401   return LC2<Scalar, Vector<Scalar>, LC2<Scalar, Node1, Node2> >(x1, x2);
00402 }
00403 
00404 /* vec - LC2 */
00405 template <class Scalar, class Node1, class Node2> inline
00406 LC2<Scalar, Vector<Scalar>, LC2<Scalar, Node1, Node2> >
00407 operator-(const Vector<Scalar>& x1, 
00408   const LC2<Scalar, Node1, Node2>& x2)
00409 {
00410   return LC2<Scalar, Vector<Scalar>, LC2<Scalar, Node1, Node2> >(x1, x2,
00411     LCSubtract);
00412 }
00413 
00414 
00415 /* LC2 + vec */
00416 template <class Scalar, class Node1, class Node2> inline
00417 LC2<Scalar, LC2<Scalar, Node1, Node2>, Vector<Scalar> >
00418 operator+(const LC2<Scalar, Node1, Node2>& x1, 
00419   const Vector<Scalar>& x2)
00420 {
00421   return LC2<Scalar, LC2<Scalar, Node1, Node2>, Vector<Scalar> >(x1, x2);
00422 }
00423 
00424 /* LC2 - vec */
00425 template <class Scalar, class Node1, class Node2> inline
00426 LC2<Scalar, LC2<Scalar, Node1, Node2>, Vector<Scalar> >
00427 operator-(const LC2<Scalar, Node1, Node2>& x1, 
00428   const Vector<Scalar>& x2)
00429 {
00430   return LC2<Scalar, LC2<Scalar, Node1, Node2>, Vector<Scalar> >(x1, x2,
00431     LCSubtract);
00432 }
00433 
00434 
00435 /*======================================================================
00436  *
00437  *    add/subtract OpTimesLC, LC2
00438  *
00439  *======================================================================*/
00440 
00441 
00442 /* OpTimesLC + LC2 */
00443 template <class Scalar, class Node0, class Node1, class Node2> inline
00444 LC2<Scalar, OpTimesLC<Scalar, Node0>, LC2<Scalar, Node1, Node2> > 
00445 operator+(const OpTimesLC<Scalar, Node0>& x1, 
00446   const LC2<Scalar, Node1, Node2>& x2)
00447 {
00448   return LC2<Scalar, OpTimesLC<Scalar, Node0>,
00449     LC2<Scalar, Node1, Node2> >(x1, x2);
00450 }
00451 
00452 /* OpTimesLC - LC2 */
00453 template <class Scalar, class Node0, class Node1, class Node2> inline
00454 LC2<Scalar, OpTimesLC<Scalar, Node0>, LC2<Scalar, Node1, Node2> > 
00455 operator-(const OpTimesLC<Scalar, Node0>& x1, 
00456   const LC2<Scalar, Node1, Node2>& x2)
00457 {
00458   return LC2<Scalar, OpTimesLC<Scalar, Node0>,
00459     LC2<Scalar, Node1, Node2> >(x1, x2, LCSubtract);
00460 }
00461 
00462 
00463 /* LC2 + OpTimesLC */
00464 template <class Scalar, class Node1, class Node2, class Node3> inline
00465 LC2<Scalar, LC2<Scalar, Node1, Node2>, OpTimesLC<Scalar, Node3> > 
00466 operator+(const LC2<Scalar, Node1, Node2>& x1, 
00467   const OpTimesLC<Scalar, Node3>& x2)
00468 {
00469   return LC2<Scalar, LC2<Scalar, Node1, Node2>, 
00470     OpTimesLC<Scalar, Node3> >(x1, x2);
00471 }
00472 
00473 /* LC2 - OpTimesLC */
00474 template <class Scalar, class Node1, class Node2, class Node3> inline
00475 LC2<Scalar, LC2<Scalar, Node1, Node2>, OpTimesLC<Scalar, Node3> > 
00476 operator-(const LC2<Scalar, Node1, Node2>& x1, 
00477   const OpTimesLC<Scalar, Node3>& x2)
00478 {
00479   return LC2<Scalar, LC2<Scalar, Node1, Node2>, 
00480     OpTimesLC<Scalar, Node3> >(x1, x2, LCSubtract);
00481 }
00482 
00483 
00484 /*======================================================================
00485  *
00486  *    add/subtract LC2, LC2
00487  *
00488  *======================================================================*/
00489   
00490 /* LC2 + LC2 */
00491 template <class Scalar, class Node1, class Node2, 
00492           class Node3, class Node4> inline
00493 LC2<Scalar, LC2<Scalar, Node1, Node2>, LC2<Scalar, Node3, Node4> >
00494 operator+(const LC2<Scalar, Node1, Node2>& x1, 
00495   const LC2<Scalar, Node3, Node4>& x2)
00496 {
00497   return LC2<Scalar, LC2<Scalar, Node1, Node2>, 
00498     LC2<Scalar, Node3, Node4> >(x1, x2);
00499 }
00500 
00501 /* LC2 - LC2 */
00502 template <class Scalar, class Node1, class Node2, 
00503           class Node3, class Node4> inline
00504 LC2<Scalar, LC2<Scalar, Node1, Node2>, LC2<Scalar, Node3, Node4> >
00505 operator-(const LC2<Scalar, Node1, Node2>& x1, 
00506   const LC2<Scalar, Node3, Node4>& x2)
00507 {
00508   return LC2<Scalar, LC2<Scalar, Node1, Node2>, 
00509     LC2<Scalar, Node3, Node4> >(x1, x2, LCSubtract);
00510 }
00511 
00512 
00513 /*======================================================================
00514  *
00515  *    assignment of [OpTimesLC, LC2] to vector
00516  *
00517  *======================================================================*/
00518   
00519   
00520 /* definition of assignment from 1-term linear combination to a vector */
00521 template <class Scalar> 
00522 template <class Node> inline
00523 Vector<Scalar>& Vector<Scalar>::operator=(const TSFExtendedOps::OpTimesLC<Scalar, Node>& x)
00524 {
00525   if (this->ptr().get()==0)
00526   {
00527     *this = x.eval();
00528   }
00529   else if (x.containsVector(this->ptr().get()))
00530   {
00531     Vector<Scalar> rtn = x.eval();
00532     acceptCopyOf(rtn);
00533   }
00534   else
00535   {
00536     x.evalInto(*this);
00537   }
00538   return *this;
00539 }
00540 
00541  
00542 /* definition of assignment from N-term linear combination to a vector */
00543 template <class Scalar>
00544 template <class Node1, class Node2> inline
00545 Vector<Scalar>& Vector<Scalar>::operator=(const TSFExtendedOps::LC2<Scalar, Node1, Node2>& x)
00546 {
00547   if (this->ptr().get()==0)
00548   {
00549     *this = x.eval();
00550   }
00551   else if (x.containsVector(this->ptr().get()))
00552   {
00553     Vector<Scalar> rtn = x.eval();
00554     acceptCopyOf(rtn);
00555   }
00556   else
00557   {
00558     x.evalInto(*this);
00559   }
00560   return *this;
00561 }
00562 
00563  
00564 
00565 /*======================================================================
00566  *
00567  *    construction of vectors from [OpTimesLC, LC2]
00568  *
00569  *======================================================================*/
00570    
00571 
00572 template <class Scalar>
00573 template <class Node1, class Node2> inline
00574 Vector<Scalar>::Vector(const TSFExtendedOps::LC2<Scalar, Node1, Node2>& x)
00575   : Handle<Thyra::VectorBase<Scalar> >(x.eval().ptr())
00576 {;}
00577 
00578 template <class Scalar> 
00579 template <class Node> inline
00580 Vector<Scalar>::Vector(const TSFExtendedOps::OpTimesLC<Scalar, Node>& x)
00581   : Handle<Thyra::VectorBase<Scalar> >(x.eval().ptr())
00582 {;}
00583 
00584 
00585 
00586 
00587   
00588 
00589 }
00590 
00591 
00592 
00593 #endif