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

Generated on Tue Oct 20 12:47:25 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7