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

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1