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