Thyra Version of the Day
Thyra_DefaultMultipliedLinearOp_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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_DEF_HPP
00043 #define THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_DEF_HPP
00044 
00045 #include "Thyra_DefaultMultipliedLinearOp_decl.hpp"
00046 #include "Thyra_AssertOp.hpp"
00047 
00048 
00049 namespace Thyra {
00050 
00051 
00052 // Constructors/initializers/accessors
00053 
00054 
00055 template<class Scalar>
00056 DefaultMultipliedLinearOp<Scalar>::DefaultMultipliedLinearOp()
00057 {}
00058 
00059 
00060 template<class Scalar>
00061 void DefaultMultipliedLinearOp<Scalar>::initialize(
00062   const ArrayView<const RCP<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>::initialize(
00075   const ArrayView<const RCP<const LinearOpBase<Scalar> > > &Ops )
00076 {
00077   const int nOps = Ops.size();
00078   Ops_.resize(nOps);
00079   for( int k = 0; k < nOps; ++k )
00080     Ops_[k].initialize(Ops[k]);
00081   validateOps();
00082   setupDefaultObjectLabel();
00083 }
00084 
00085 
00086 template<class Scalar>
00087 void DefaultMultipliedLinearOp<Scalar>::uninitialize()
00088 {
00089   Ops_.resize(0);
00090   setupDefaultObjectLabel();
00091 }
00092 
00093 
00094 // Overridden form MultipliedLinearOpBase
00095 
00096 
00097 template<class Scalar>
00098 int DefaultMultipliedLinearOp<Scalar>::numOps() const
00099 {
00100   return Ops_.size();
00101 }
00102 
00103 
00104 template<class Scalar>
00105 bool DefaultMultipliedLinearOp<Scalar>::opIsConst(const int k) const
00106 {
00107 #ifdef TEUCHOS_DEBUG
00108   TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00109 #endif
00110   return Ops_[k].isConst();
00111 }
00112 
00113 
00114 template<class Scalar>
00115 RCP<LinearOpBase<Scalar> >
00116 DefaultMultipliedLinearOp<Scalar>::getNonconstOp(const int k)
00117 {
00118 #ifdef TEUCHOS_DEBUG
00119   TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00120 #endif
00121   return Ops_[k].getNonconstObj();
00122 }
00123 
00124 
00125 template<class Scalar>
00126 RCP<const LinearOpBase<Scalar> >
00127 DefaultMultipliedLinearOp<Scalar>::getOp(const int k) const
00128 {
00129 #ifdef TEUCHOS_DEBUG
00130   TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= k && k < numOps() ) );
00131 #endif
00132   return Ops_[k];
00133 }
00134 
00135 
00136 // Overridden from LinearOpBase
00137 
00138 
00139 template<class Scalar>
00140 RCP< const VectorSpaceBase<Scalar> >
00141 DefaultMultipliedLinearOp<Scalar>::range() const
00142 {
00143   assertInitialized();
00144   return getOp(0)->range();
00145 }
00146 
00147 
00148 template<class Scalar>
00149 RCP< const VectorSpaceBase<Scalar> >
00150 DefaultMultipliedLinearOp<Scalar>::domain() const
00151 {
00152   assertInitialized();
00153   return getOp(numOps()-1)->domain();
00154 }
00155 
00156 
00157 template<class Scalar>
00158 RCP<const LinearOpBase<Scalar> >
00159 DefaultMultipliedLinearOp<Scalar>::clone() const
00160 {
00161   return Teuchos::null; // Not supported yet but could be!
00162 }
00163 
00164 
00165 // Overridden from Teuchos::Describable
00166 
00167                                                 
00168 template<class Scalar>
00169 std::string DefaultMultipliedLinearOp<Scalar>::description() const
00170 {
00171   assertInitialized();
00172   typedef Teuchos::ScalarTraits<Scalar>  ST;
00173   std::ostringstream oss;
00174   oss << Teuchos::Describable::description() << "{numOps = "<<numOps()<<"}";
00175   return oss.str();
00176 }
00177 
00178 template<class Scalar>
00179 void DefaultMultipliedLinearOp<Scalar>::describe(
00180   Teuchos::FancyOStream &out_arg,
00181   const Teuchos::EVerbosityLevel verbLevel
00182   ) const
00183 {
00184   typedef Teuchos::ScalarTraits<Scalar>  ST;
00185   using Teuchos::FancyOStream;
00186   using Teuchos::OSTab;
00187   assertInitialized();
00188   RCP<FancyOStream> out = rcp(&out_arg,false);
00189   OSTab tab(out);
00190   const int nOps = Ops_.size();
00191   switch(verbLevel) {
00192     case Teuchos::VERB_DEFAULT:
00193     case Teuchos::VERB_LOW:
00194       *out << this->description() << std::endl;
00195       break;
00196     case Teuchos::VERB_MEDIUM:
00197     case Teuchos::VERB_HIGH:
00198     case Teuchos::VERB_EXTREME:
00199     {
00200       *out
00201         << Teuchos::Describable::description() << "{"
00202         << "rangeDim=" << this->range()->dim()
00203         << ",domainDim="<< this->domain()->dim() << "}\n";
00204       OSTab tab2(out);
00205       *out
00206         <<  "numOps = "<< nOps << std::endl
00207         <<  "Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:\n";
00208       OSTab tab3(out);
00209       for( int k = 0; k < nOps; ++k ) {
00210         *out << "Op["<<k<<"] = " << Teuchos::describe(*getOp(k),verbLevel);
00211       }
00212       break;
00213     }
00214     default:
00215       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00216   }
00217 }
00218 
00219 
00220 // protected
00221 
00222 
00223 // Overridden from LinearOpBase
00224 
00225 
00226 template<class Scalar>
00227 bool DefaultMultipliedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00228 {
00229   bool overallOpSupported = true;
00230   for( int k = 0; k < static_cast<int>(Ops_.size()); ++k )
00231     if(!Thyra::opSupported(*getOp(k),M_trans)) overallOpSupported = false;
00232   return overallOpSupported;
00233   // ToDo: Cache these?
00234 }
00235 
00236 
00237 template<class Scalar>
00238 void DefaultMultipliedLinearOp<Scalar>::applyImpl(
00239   const EOpTransp M_trans,
00240   const MultiVectorBase<Scalar> &X,
00241   const Ptr<MultiVectorBase<Scalar> > &Y,
00242   const Scalar alpha,
00243   const Scalar beta
00244   ) const
00245 {
00246   using Teuchos::rcpFromPtr;
00247   using Teuchos::rcpFromRef;
00248 #ifdef TEUCHOS_DEBUG
00249   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00250     "DefaultMultipliedLinearOp<Scalar>::apply(...)", *this, M_trans, X, &*Y
00251     );
00252 #endif // TEUCHOS_DEBUG  
00253   const int nOps = Ops_.size();
00254   const Ordinal m = X.domain()->dim();
00255   if( real_trans(M_trans)==NOTRANS ) {
00256     //
00257     // Y = alpha * M * X + beta*Y
00258     // =>
00259     // Y = alpha * op(Op[0]) * op(Op[1]) * ... * op(Op[numOps-1]) * X + beta*Y
00260     //
00261     RCP<MultiVectorBase<Scalar> > T_kp1, T_k; // Temporary propagated between loops 
00262     for( int k = nOps-1; k >= 0; --k ) {
00263       RCP<MultiVectorBase<Scalar> > Y_k;
00264       RCP<const MultiVectorBase<Scalar> > X_k;
00265       if(k==0) Y_k = rcpFromPtr(Y);  else Y_k = T_k = createMembers(getOp(k)->range(), m);
00266       if(k==nOps-1) X_k = rcpFromRef(X); else X_k = T_kp1;
00267       if( k > 0 )
00268         Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
00269       else
00270         Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
00271       T_kp1 = T_k;
00272     }
00273   }
00274   else {
00275     //
00276     // Y = alpha * M' * X + beta*Y
00277     // =>
00278     // Y = alpha * Op[numOps-1]' * Op[1]' * ... * Op[0]' * X + beta * Y
00279     //
00280     RCP<MultiVectorBase<Scalar> > T_km1, T_k; // Temporary propagated between loops 
00281     for( int k = 0; k <= nOps-1; ++k ) {
00282       RCP<MultiVectorBase<Scalar> >         Y_k;
00283       RCP<const MultiVectorBase<Scalar> >   X_k;
00284       if(k==nOps-1) Y_k = rcpFromPtr(Y);  else Y_k = T_k = createMembers(getOp(k)->domain(), m);
00285       if(k==0) X_k = rcpFromRef(X); else X_k = T_km1;
00286       if( k < nOps-1 )
00287         Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
00288       else
00289         Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
00290       T_km1 = T_k;
00291     }
00292   }
00293 }
00294 
00295 
00296 // private
00297 
00298 
00299 template<class Scalar>
00300 void DefaultMultipliedLinearOp<Scalar>::validateOps()
00301 {
00302   using Teuchos::toString;
00303 #ifdef TEUCHOS_DEBUG
00304   try {
00305     const int nOps = Ops_.size();
00306     for( int k = 0; k < nOps; ++k ) {
00307       TEUCHOS_TEST_FOR_EXCEPT( Ops_[k]().get() == NULL );
00308       if( k < nOps-1 ) {
00309         THYRA_ASSERT_LINEAR_OP_TIMES_LINEAR_OP_SPACES_NAMES(
00310           "DefaultMultipliedLinearOp<Scalar>::initialize(...)"
00311           ,*Ops_[k],NOTRANS,("Ops["+toString(k)+"]")
00312           ,*Ops_[k+1],NOTRANS,("Ops["+toString(k+1)+"]")
00313           );
00314       }
00315     }
00316   }
00317   catch(...) {
00318     uninitialize();
00319     throw;
00320   }
00321 #endif
00322 }
00323 
00324 
00325 template<class Scalar>
00326 void DefaultMultipliedLinearOp<Scalar>::setupDefaultObjectLabel()
00327 {
00328   std::ostringstream label;
00329   const int nOps = Ops_.size();
00330   for( int k = 0; k < nOps; ++k ) {
00331     std::string Op_k_label = Ops_[k]->getObjectLabel();
00332     if (Op_k_label.length() == 0)
00333       Op_k_label = "ANYM";
00334     if (k > 0)
00335       label << "*";
00336     label << "("<<Op_k_label<<")";
00337   }
00338   this->setObjectLabel(label.str());
00339   validateOps();
00340 }
00341 
00342 
00343 } // end namespace Thyra
00344 
00345 
00346 //
00347 // Nonmember implementations
00348 //
00349 
00350 template<class Scalar>
00351 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00352 Thyra::nonconstMultiply(
00353   const RCP<LinearOpBase<Scalar> > &A,
00354   const RCP<LinearOpBase<Scalar> > &B,
00355   const std::string &M_label
00356   )
00357 {
00358   using Teuchos::tuple;
00359   RCP<DefaultMultipliedLinearOp<Scalar> > multOp = 
00360     defaultMultipliedLinearOp<Scalar>(tuple(A, B)());
00361   if(M_label.length())
00362     multOp->setObjectLabel(M_label);
00363   return multOp;
00364 }
00365 
00366 
00367 template<class Scalar>
00368 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00369 Thyra::multiply(
00370   const RCP<const LinearOpBase<Scalar> > &A,
00371   const RCP<const LinearOpBase<Scalar> > &B,
00372   const std::string &M_label
00373   )
00374 {
00375   using Teuchos::tuple;
00376   RCP<DefaultMultipliedLinearOp<Scalar> > multOp =
00377     defaultMultipliedLinearOp<Scalar>(tuple(A, B)());
00378   if(M_label.length())
00379     multOp->setObjectLabel(M_label);
00380   return multOp;
00381 }
00382 
00383 
00384 template<class Scalar>
00385 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
00386 Thyra::multiply(
00387   const RCP<const LinearOpBase<Scalar> > &A,
00388   const RCP<const LinearOpBase<Scalar> > &B,
00389   const RCP<const LinearOpBase<Scalar> > &C,
00390   const std::string &M_label
00391   )
00392 {
00393   using Teuchos::tuple;
00394   RCP<DefaultMultipliedLinearOp<Scalar> > multOp =
00395     defaultMultipliedLinearOp<Scalar>(tuple(A, B, C)());
00396   if(M_label.length())
00397     multOp->setObjectLabel(M_label);
00398   return multOp;
00399 }
00400 
00401 
00402 //
00403 // Explicit instantiation macro
00404 //
00405 // Must be expanded from within the Thyra namespace!
00406 //
00407 
00408 
00409 #define THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_INSTANT(SCALAR) \
00410   \
00411   template class DefaultMultipliedLinearOp<SCALAR >; \
00412  \
00413   template RCP<LinearOpBase<SCALAR > > \
00414   nonconstMultiply( \
00415     const RCP<LinearOpBase<SCALAR > > &A, \
00416     const RCP<LinearOpBase<SCALAR > > &B, \
00417     const std::string &M_label \
00418     ); \
00419    \
00420   template RCP<const LinearOpBase<SCALAR > > \
00421   multiply( \
00422     const RCP<const LinearOpBase<SCALAR > > &A, \
00423     const RCP<const LinearOpBase<SCALAR > > &B, \
00424     const std::string &M_label \
00425     ); \
00426    \
00427   template RCP<const LinearOpBase<SCALAR > > \
00428   multiply( \
00429     const RCP<const LinearOpBase<SCALAR > > &A, \
00430     const RCP<const LinearOpBase<SCALAR > > &B, \
00431     const RCP<const LinearOpBase<SCALAR > > &C, \
00432     const std::string &M_label \
00433     ); \
00434 
00435 
00436 #endif  // THYRA_DEFAULT_MULTIPLIED_LINEAR_OP_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines