Thyra Version of the Day
Thyra_DefaultAddedLinearOp_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
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 THYRA_DEFAULT_ADDED_LINEAR_OP_DEF_HPP
00030 #define THYRA_DEFAULT_ADDED_LINEAR_OP_DEF_HPP
00031 
00032 #include "Thyra_DefaultAddedLinearOp_decl.hpp"
00033 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00034 #include "Thyra_AssertOp.hpp"
00035 #include "Teuchos_Utils.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00041 // Constructors/initializers/accessors
00042 
00043 
00044 template<class Scalar>
00045 DefaultAddedLinearOp<Scalar>::DefaultAddedLinearOp()
00046 {}
00047 
00048 
00049 template<class Scalar>
00050 DefaultAddedLinearOp<Scalar>::DefaultAddedLinearOp(
00051   const ArrayView<const RCP<LinearOpBase<Scalar> > > &Ops )
00052 {
00053   initialize(Ops);
00054 }
00055 
00056 
00057 template<class Scalar>
00058 DefaultAddedLinearOp<Scalar>::DefaultAddedLinearOp(
00059   const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
00060 {
00061   initialize(Ops);
00062 }
00063 
00064 
00065 template<class Scalar>
00066 void DefaultAddedLinearOp<Scalar>::initialize(
00067   const ArrayView<const RCP<LinearOpBase<Scalar> > > &Ops )
00068 {
00069   const int l_numOps = Ops.size();
00070   Ops_.resize(l_numOps);
00071   for( int k = 0; k < l_numOps; ++k )
00072     Ops_[k].initialize(Ops[k]);
00073   validateOps();
00074   setupDefaultObjectLabel();
00075 }
00076 
00077 
00078 template<class Scalar>
00079 void DefaultAddedLinearOp<Scalar>::initialize(
00080   const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
00081 {
00082   const int l_numOps = Ops.size();
00083   Ops_.resize(l_numOps);
00084   for( int k = 0; k < l_numOps; ++k )
00085     Ops_[k].initialize(Ops[k]);
00086   validateOps();
00087   setupDefaultObjectLabel();
00088 }
00089 
00090 
00091 template<class Scalar>
00092 void DefaultAddedLinearOp<Scalar>::uninitialize()
00093 {
00094   Ops_.resize(0);
00095   setupDefaultObjectLabel();
00096 }
00097 
00098 
00099 // Overridden form AddedLinearOpBase
00100 
00101 
00102 template<class Scalar>
00103 int DefaultAddedLinearOp<Scalar>::numOps() const
00104 {
00105   return Ops_.size();
00106 }
00107 
00108 
00109 template<class Scalar>
00110 bool DefaultAddedLinearOp<Scalar>::opIsConst(const int k) const
00111 {
00112 #ifdef TEUCHOS_DEBUG
00113   TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00114 #endif
00115   return Ops_[k].isConst();
00116 }
00117 
00118 
00119 template<class Scalar>
00120 Teuchos::RCP<LinearOpBase<Scalar> >
00121 DefaultAddedLinearOp<Scalar>::getNonconstOp(const int k)
00122 {
00123 #ifdef TEUCHOS_DEBUG
00124   TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00125 #endif
00126   return Ops_[k].getNonconstObj();
00127 }
00128 
00129 
00130 template<class Scalar>
00131 Teuchos::RCP<const LinearOpBase<Scalar> >
00132 DefaultAddedLinearOp<Scalar>::getOp(const int k) const
00133 {
00134 #ifdef TEUCHOS_DEBUG
00135   TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00136 #endif
00137   return Ops_[k];
00138 }
00139 
00140 
00141 // Overridden from LinearOpBase
00142 
00143 
00144 template<class Scalar>
00145 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00146 DefaultAddedLinearOp<Scalar>::range() const
00147 {
00148   assertInitialized();
00149   return getOp(0)->range();
00150 }
00151 
00152 
00153 template<class Scalar>
00154 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00155 DefaultAddedLinearOp<Scalar>::domain() const
00156 {
00157   assertInitialized();
00158   return getOp(numOps()-1)->domain();
00159 }
00160 
00161 
00162 template<class Scalar>
00163 Teuchos::RCP<const LinearOpBase<Scalar> >
00164 DefaultAddedLinearOp<Scalar>::clone() const
00165 {
00166   return Teuchos::null; // Not supported yet but could be!
00167 }
00168 
00169 
00170 // Overridden from Teuchos::Describable
00171 
00172                                                 
00173 template<class Scalar>
00174 std::string DefaultAddedLinearOp<Scalar>::description() const
00175 {
00176   assertInitialized();
00177   std::ostringstream oss;
00178   oss << Teuchos::Describable::description() << "{numOps = "<<numOps()<<"}";
00179   return oss.str();
00180 }
00181 
00182 
00183 template<class Scalar>
00184 void DefaultAddedLinearOp<Scalar>::describe(
00185   Teuchos::FancyOStream                &out_arg
00186   ,const Teuchos::EVerbosityLevel      verbLevel
00187   ) const
00188 {
00189   typedef Teuchos::ScalarTraits<Scalar>  ST;
00190   using Teuchos::RCP;
00191   using Teuchos::FancyOStream;
00192   using Teuchos::OSTab;
00193   assertInitialized();
00194   RCP<FancyOStream> out = rcp(&out_arg,false);
00195   OSTab tab(out);
00196   const int l_numOps = Ops_.size();
00197   switch(verbLevel) {
00198     case Teuchos::VERB_DEFAULT:
00199     case Teuchos::VERB_LOW:
00200       *out << this->description() << std::endl;
00201       break;
00202     case Teuchos::VERB_MEDIUM:
00203     case Teuchos::VERB_HIGH:
00204     case Teuchos::VERB_EXTREME:
00205     {
00206       *out
00207         << Teuchos::Describable::description() << "{"
00208         << "rangeDim=" << this->range()->dim()
00209         << ",domainDim="<< this->domain()->dim() << "}\n";
00210       OSTab tab2(out);
00211       *out
00212         <<  "numOps="<< l_numOps << std::endl
00213         <<  "Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:\n";
00214       tab.incrTab();
00215       for( int k = 0; k < l_numOps; ++k ) {
00216         *out << "Op["<<k<<"] = " << Teuchos::describe(*getOp(k),verbLevel);
00217       }
00218       break;
00219     }
00220     default:
00221       TEST_FOR_EXCEPT(true); // Should never get here!
00222   }
00223 }
00224 
00225 
00226 // Deprecated
00227 
00228 
00229 template<class Scalar>
00230 DefaultAddedLinearOp<Scalar>::DefaultAddedLinearOp(
00231   const int numOps_in,
00232   const RCP<LinearOpBase<Scalar> > Ops[]
00233   )
00234 {
00235   initialize(Teuchos::arrayView(Ops, numOps_in));
00236 }
00237 
00238 
00239 template<class Scalar>
00240 DefaultAddedLinearOp<Scalar>::DefaultAddedLinearOp(
00241   const int numOps_in,
00242   const RCP<const LinearOpBase<Scalar> > Ops[]
00243   )
00244 {
00245   initialize(Teuchos::arrayView(Ops, numOps_in));
00246 }
00247 
00248 
00249 // protected
00250 
00251 
00252 // Overridden from LinearOpBase
00253 
00254 
00255 template<class Scalar>
00256 bool DefaultAddedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00257 {
00258   bool isOpSupported = true;
00259   for( int k = 0; k < static_cast<int>(Ops_.size()); ++k )
00260     if(!Thyra::opSupported(*getOp(k),M_trans)) isOpSupported = false;
00261   return isOpSupported;
00262   // ToDo: Cache these?
00263 }
00264 
00265 
00266 template<class Scalar>
00267 void DefaultAddedLinearOp<Scalar>::applyImpl(
00268   const EOpTransp M_trans,
00269   const MultiVectorBase<Scalar> &X,
00270   const Ptr<MultiVectorBase<Scalar> > &Y,
00271   const Scalar alpha,
00272   const Scalar beta
00273   ) const
00274 {
00275   typedef Teuchos::ScalarTraits<Scalar> ST;
00276 #ifdef TEUCHOS_DEBUG
00277   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00278     "DefaultAddedLinearOp<Scalar>::apply(...)", *this, M_trans, X, &*Y
00279     );
00280 #endif // TEUCHOS_DEBUG  
00281   //
00282   // Y = alpha * op(M) * X + beta*Y
00283   //
00284   // =>
00285   //
00286   // Y = beta*Y + sum(alpha*op(Op[j])*X),j=0...numOps-1)
00287   //
00288   const int l_numOps = Ops_.size();
00289   for( int j = 0; j < l_numOps; ++j )
00290     Thyra::apply(*getOp(j), M_trans, X, Y, alpha, j==0 ? beta : ST::one());
00291 }
00292 
00293 
00294 // private
00295 
00296 
00297 template<class Scalar>
00298 void DefaultAddedLinearOp<Scalar>::validateOps()
00299 {
00300   typedef std::string s;
00301   using Teuchos::toString;
00302 #ifdef TEUCHOS_DEBUG
00303   try {
00304     const int l_numOps = Ops_.size();
00305     for( int k = 0; k < l_numOps; ++k ) {
00306       TEST_FOR_EXCEPT( Ops_[k]().get() == NULL );
00307       if( k > 0 ) {
00308         THYRA_ASSERT_LINEAR_OP_PLUS_LINEAR_OP_SPACES_NAMES(
00309           "DefaultMultipliedLinearOp<Scalar>::initialize(...)"
00310           ,*Ops_[0], NOTRANS, ("Ops[0]")
00311           ,*Ops_[k], NOTRANS, ("Ops["+toString(k)+"]")
00312           );
00313       }
00314     }
00315   }
00316   catch(...) {
00317     uninitialize();
00318     throw;
00319   }
00320 #endif
00321 }
00322 
00323 
00324 template<class Scalar>
00325 void DefaultAddedLinearOp<Scalar>::setupDefaultObjectLabel()
00326 {
00327   std::ostringstream label;
00328   const int l_numOps = Ops_.size();
00329   for( int k = 0; k < l_numOps; ++k ) {
00330     std::string Op_k_label = Ops_[k]->getObjectLabel();
00331     if (Op_k_label.length() == 0)
00332       Op_k_label = "ANYM";
00333     if (k > 0)
00334       label << "+";
00335     label << "("<<Op_k_label<<")";
00336   }
00337   this->setObjectLabel(label.str());
00338   validateOps();
00339 }
00340 
00341 
00342 } // end namespace Thyra
00343 
00344 
00345 template<class Scalar>
00346 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00347 Thyra::nonconstAdd(
00348   const RCP<LinearOpBase<Scalar> > &A,
00349   const RCP<LinearOpBase<Scalar> > &B,
00350   const std::string &label
00351   )
00352 {
00353   using Teuchos::tuple;
00354   RCP<LinearOpBase<Scalar> > alo =
00355     defaultAddedLinearOp<Scalar>(
00356       tuple<RCP<LinearOpBase<Scalar> > >(A, B)()
00357       );
00358   if (label.length())
00359     alo->setObjectLabel(label);
00360   return alo;
00361 }
00362 
00363 
00364 template<class Scalar>
00365 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00366 Thyra::add(
00367   const RCP<const LinearOpBase<Scalar> > &A,
00368   const RCP<const LinearOpBase<Scalar> > &B,
00369   const std::string &label
00370   )
00371 {
00372   using Teuchos::tuple;
00373   RCP<LinearOpBase<Scalar> > alo =
00374     defaultAddedLinearOp<Scalar>(
00375       tuple<RCP<const LinearOpBase<Scalar> > >(A, B)()
00376       );
00377   if (label.length())
00378     alo->setObjectLabel(label);
00379   return alo;
00380 }
00381 
00382 
00383 template<class Scalar>
00384 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00385 Thyra::nonconstSubtract(
00386   const RCP<LinearOpBase<Scalar> > &A,
00387   const RCP<LinearOpBase<Scalar> > &B,
00388   const std::string &label
00389   )
00390 {
00391   typedef ScalarTraits<Scalar> ST;
00392   using Teuchos::tuple;
00393   RCP<LinearOpBase<Scalar> > alo =
00394     defaultAddedLinearOp<Scalar>(
00395       tuple<RCP<LinearOpBase<Scalar> > >(
00396         A, nonconstScale<Scalar>(-ST::one(),B) )()
00397       );
00398   if (label.length())
00399     alo->setObjectLabel(label);
00400   return alo;
00401 }
00402 
00403 
00404 template<class Scalar>
00405 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00406 Thyra::subtract(
00407   const RCP<const LinearOpBase<Scalar> > &A,
00408   const RCP<const LinearOpBase<Scalar> > &B,
00409   const std::string &label
00410   )
00411 {
00412   typedef ScalarTraits<Scalar> ST;
00413   using Teuchos::tuple;
00414   RCP<LinearOpBase<Scalar> > alo =
00415     defaultAddedLinearOp<Scalar>(
00416       tuple<RCP<const LinearOpBase<Scalar> > >(
00417           A, scale<Scalar>(-ST::one(),B)
00418         )()
00419       );
00420   if (label.length())
00421     alo->setObjectLabel(label);
00422   return alo;
00423 }
00424 
00425 
00426 //
00427 // Explicit instantiation macro
00428 //
00429 
00430 #define THYRA_DEFAULT_ADDED_LINEAR_OP_INSTANT(SCALAR) \
00431   template class DefaultAddedLinearOp<SCALAR >; \
00432   \
00433   template RCP<LinearOpBase<SCALAR > > \
00434   nonconstAdd( \
00435     const RCP<LinearOpBase<SCALAR > > &A, \
00436     const RCP<LinearOpBase<SCALAR > > &B, \
00437     const std::string &label \
00438     ); \
00439    \
00440   template RCP<const LinearOpBase<SCALAR > > \
00441   add( \
00442     const RCP<const LinearOpBase<SCALAR > > &A, \
00443     const RCP<const LinearOpBase<SCALAR > > &B, \
00444     const std::string &label \
00445     ); \
00446    \
00447   template RCP<LinearOpBase<SCALAR > > \
00448   nonconstSubtract( \
00449     const RCP<LinearOpBase<SCALAR > > &A, \
00450     const RCP<LinearOpBase<SCALAR > > &B, \
00451     const std::string &label \
00452     ); \
00453    \
00454   template RCP<const LinearOpBase<SCALAR > > \
00455   subtract( \
00456     const RCP<const LinearOpBase<SCALAR > > &A, \
00457     const RCP<const LinearOpBase<SCALAR > > &B, \
00458     const std::string &label \
00459     );
00460 
00461 
00462 #endif  // THYRA_DEFAULT_ADDED_LINEAR_OP_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines