Thyra_EpetraLinearOp.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //               Thyra: Trilinos Solver Framework Core
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 #include "Thyra_EpetraLinearOp.hpp"
00030 #include "Thyra_EpetraThyraWrappers.hpp"
00031 #include "Teuchos_dyn_cast.hpp"
00032 #include "Teuchos_TestForException.hpp"
00033 #include "Teuchos_getConst.hpp"
00034 
00035 #include "Epetra_Map.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_Operator.h"
00038 
00039 namespace Thyra {
00040 
00041 // Constructors / initializers / accessors
00042 
00043 EpetraLinearOp::EpetraLinearOp()
00044 {}
00045 
00046 EpetraLinearOp::EpetraLinearOp(
00047   const Teuchos::RefCountPtr<Epetra_Operator>                        &op
00048   ,ETransp                                                           opTrans
00049   ,EApplyEpetraOpAs                                                  applyAs
00050   ,EAdjointEpetraOp                                                  adjointSupport
00051   ,const Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >   &spmdRange
00052   ,const Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >   &spmdDomain
00053   )
00054 {
00055   initialize(op,opTrans,applyAs,adjointSupport,spmdRange,spmdDomain);
00056 }
00057 
00058 void EpetraLinearOp::initialize(
00059   const Teuchos::RefCountPtr<Epetra_Operator>                        &op
00060   ,ETransp                                                           opTrans
00061   ,EApplyEpetraOpAs                                                  applyAs
00062   ,EAdjointEpetraOp                                                  adjointSupport
00063   ,const Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >   &spmdRange
00064   ,const Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >   &spmdDomain
00065   )
00066 {
00067   // Validate input, allocate spaces, validate ...
00068 #ifdef TEUCHOS_DEBUG
00069   TEST_FOR_EXCEPTION( op.get()==NULL, std::invalid_argument, "EpetraLinearOp::initialize(...): Error!" );
00070 #endif
00071   Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> > range, domain;
00072   if(spmdRange.get())
00073     range = spmdRange;
00074   else
00075     range = ( applyAs==EPETRA_OP_APPLY_APPLY ? allocateRange(op,opTrans)  : allocateDomain(op,opTrans) );
00076   if(spmdDomain.get())
00077     domain = spmdDomain;
00078   else
00079     domain  = ( applyAs==EPETRA_OP_APPLY_APPLY ? allocateDomain(op,opTrans) : allocateRange(op,opTrans)  );
00080   Teuchos::RefCountPtr<const ScalarProdVectorSpaceBase<EpetraLinearOp::Scalar> >
00081     sp_range = Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<EpetraLinearOp::Scalar> >(range),
00082     sp_domain = Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<EpetraLinearOp::Scalar> >(domain);
00083   // Set data (no exceptions should be thrown now)
00084   op_      = op;
00085   opTrans_ = opTrans;
00086   applyAs_ = applyAs;
00087   adjointSupport_ = adjointSupport;
00088   range_ = range;
00089   domain_ = domain;
00090   sp_range_ = sp_range;
00091   sp_domain_ = sp_domain;
00092 }
00093 
00094 void EpetraLinearOp::uninitialize(
00095   Teuchos::RefCountPtr<Epetra_Operator>                       *op
00096   ,ETransp                                                    *opTrans
00097   ,EApplyEpetraOpAs                                           *applyAs
00098   ,EAdjointEpetraOp                                           *adjointSupport
00099   ,Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >  *spmdRange
00100   ,Teuchos::RefCountPtr< const SpmdVectorSpaceBase<Scalar> >  *spmdDomain
00101   )
00102 {
00103 
00104   if(op)      *op      = op_;
00105   if(opTrans) *opTrans = opTrans_;
00106   if(applyAs) *applyAs = applyAs_;
00107   if(adjointSupport) *adjointSupport = adjointSupport_;
00108   if(spmdRange) *spmdRange = range_;
00109   if(spmdDomain) *spmdDomain = domain_;
00110 
00111   op_      = Teuchos::null;
00112   opTrans_ = NOTRANS;
00113   applyAs_ = EPETRA_OP_APPLY_APPLY;
00114   adjointSupport_ = EPETRA_OP_ADJOINT_SUPPORTED;
00115   range_   = Teuchos::null;
00116   domain_  = Teuchos::null;
00117   sp_range_   = Teuchos::null;
00118   sp_domain_  = Teuchos::null;
00119 
00120 }
00121 
00122 Teuchos::RefCountPtr< const SpmdVectorSpaceBase<EpetraLinearOp::Scalar> >
00123 EpetraLinearOp::spmdRange() const
00124 {
00125   return range_;
00126 }
00127 
00128 Teuchos::RefCountPtr< const SpmdVectorSpaceBase<EpetraLinearOp::Scalar> >
00129 EpetraLinearOp::spmdDomain() const
00130 {
00131   return domain_;
00132 }
00133 
00134 Teuchos::RefCountPtr<Epetra_Operator>
00135 EpetraLinearOp::epetra_op() 
00136 {
00137   return op_;
00138 }
00139 
00140 Teuchos::RefCountPtr<const Epetra_Operator>
00141 EpetraLinearOp::epetra_op() const 
00142 {
00143   return op_;
00144 }
00145 
00146 // Overridden from EpetraLinearOpBase
00147 
00148 void EpetraLinearOp::getEpetraOpView(
00149   Teuchos::RefCountPtr<Epetra_Operator>   *epetraOp
00150   ,ETransp                                *epetraOpTransp
00151   ,EApplyEpetraOpAs                       *epetraOpApplyAs
00152   ,EAdjointEpetraOp                       *epetraOpAdjointSupport
00153   )
00154 {
00155 #ifdef TEUCHOS_DEBUG
00156   TEST_FOR_EXCEPT(epetraOp==NULL);
00157   TEST_FOR_EXCEPT(epetraOpTransp==NULL);
00158   TEST_FOR_EXCEPT(epetraOpApplyAs==NULL);
00159   TEST_FOR_EXCEPT(epetraOpAdjointSupport==NULL);
00160 #endif
00161   *epetraOp = op_;
00162   *epetraOpTransp = opTrans_;
00163   *epetraOpApplyAs = applyAs_;
00164   *epetraOpAdjointSupport = adjointSupport_;
00165 }
00166 
00167 void EpetraLinearOp::getEpetraOpView(
00168   Teuchos::RefCountPtr<const Epetra_Operator>   *epetraOp
00169   ,ETransp                                      *epetraOpTransp
00170   ,EApplyEpetraOpAs                             *epetraOpApplyAs
00171   ,EAdjointEpetraOp                             *epetraOpAdjointSupport
00172   ) const
00173 {
00174 #ifdef TEUCHOS_DEBUG
00175   TEST_FOR_EXCEPT(epetraOp==NULL);
00176   TEST_FOR_EXCEPT(epetraOpTransp==NULL);
00177   TEST_FOR_EXCEPT(epetraOpApplyAs==NULL);
00178   TEST_FOR_EXCEPT(epetraOpAdjointSupport==NULL);
00179 #endif
00180   *epetraOp = op_;
00181   *epetraOpTransp = opTrans_;
00182   *epetraOpApplyAs = applyAs_;
00183   *epetraOpAdjointSupport = adjointSupport_;
00184 }
00185 
00186 // Overridden from SingleScalarLinearOpBase
00187 
00188 bool EpetraLinearOp::opSupported(ETransp M_trans) const
00189 {
00190   return ( M_trans == NOTRANS ? true : adjointSupport_==EPETRA_OP_ADJOINT_SUPPORTED );
00191 }
00192 
00193 // Overridden from EuclideanLinearOpBase
00194 
00195 Teuchos::RefCountPtr<const ScalarProdVectorSpaceBase<EpetraLinearOp::Scalar> >
00196 EpetraLinearOp::rangeScalarProdVecSpc() const
00197 {
00198   return sp_range_;
00199 }
00200 
00201 Teuchos::RefCountPtr<const ScalarProdVectorSpaceBase<EpetraLinearOp::Scalar> >
00202 EpetraLinearOp::domainScalarProdVecSpc() const
00203 {
00204   return sp_domain_;
00205 }
00206 
00207 void EpetraLinearOp::euclideanApply(
00208   const ETransp                     M_trans
00209   ,const MultiVectorBase<Scalar>    &X_in
00210   ,MultiVectorBase<Scalar>          *Y_inout
00211   ,const Scalar                     alpha
00212   ,const Scalar                     beta
00213   ) const
00214 {
00215   const ETransp real_M_trans = real_trans(M_trans);
00216 #ifdef TEUCHOS_DEBUG
00217   // ToDo: Assert vector spaces!
00218   TEST_FOR_EXCEPTION(
00219     real_M_trans==TRANS && adjointSupport_==EPETRA_OP_ADJOINT_UNSUPPORTED
00220     ,Exceptions::OpNotSupported
00221     ,"EpetraLinearOp::apply(...): *this was informed that adjoints are not supported when initialized." 
00222     );
00223 #endif
00224   //
00225   // Get Epetra_MultiVector objects for the arguments
00226   //
00227   Teuchos::RefCountPtr<const Epetra_MultiVector>
00228     X = get_Epetra_MultiVector(
00229       real_M_trans==NOTRANS ? getDomainMap() : getRangeMap()
00230       ,Teuchos::rcp(&X_in,false)
00231       );
00232   Teuchos::RefCountPtr<Epetra_MultiVector>
00233     Y;
00234   if( beta == 0 ) {
00235     Y = get_Epetra_MultiVector(
00236       real_M_trans==NOTRANS ? getRangeMap() : getDomainMap()
00237       ,Teuchos::rcp(Y_inout,false)
00238       );
00239   }
00240   //
00241   // Set the operator mode
00242   //
00243   /* We need to save the transpose state here, and then reset it after 
00244    * application. The reason for this is that if we later apply the 
00245    * operator outside Thyra (in Aztec, for instance), it will remember
00246    * the transpose flag set here. */
00247   bool oldState = op_->UseTranspose();
00248   op_->SetUseTranspose( real_trans(trans_trans(opTrans_,M_trans)) == NOTRANS ? false : true );
00249   //
00250   // Perform the operation
00251   //
00252   if( beta == 0.0 ) {
00253     // Y = M * X
00254     if( applyAs_ == EPETRA_OP_APPLY_APPLY )
00255       op_->Apply( *X, *Y );
00256     else if( applyAs_ == EPETRA_OP_APPLY_APPLY_INVERSE )
00257       op_->ApplyInverse( *X, *Y );
00258     else
00259       TEST_FOR_EXCEPT(true);
00260     // Y = alpha * Y
00261     if( alpha != 1.0 ) Y->Scale(alpha);
00262   }
00263   else {
00264     // Y_inout = beta * Y_inout
00265     if(beta != 0.0) scale( beta, Y_inout );
00266     else assign( Y_inout, 0.0 );
00267     // T = M * X
00268     Epetra_MultiVector T(
00269       real_M_trans == NOTRANS ? op_->OperatorRangeMap() : op_->OperatorDomainMap()
00270       ,X_in.domain()->dim()
00271       ,false
00272       );
00273     if( applyAs_ == EPETRA_OP_APPLY_APPLY )
00274       op_->Apply( *X, T );
00275     else if( applyAs_ == EPETRA_OP_APPLY_APPLY_INVERSE )
00276       op_->ApplyInverse( *X, T );
00277     else
00278       TEST_FOR_EXCEPT(true);
00279     // Y_inout += alpha * T
00280     update(
00281       alpha
00282       ,*create_MultiVector(
00283         Teuchos::rcp(&Teuchos::getConst(T),false)
00284         ,Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<Scalar> >(Y_inout->range(),true)
00285         ,Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<Scalar> >(Y_inout->domain(),true)
00286         )
00287       ,Y_inout
00288       );
00289   }
00290   // Reset the transpose state
00291   op_->SetUseTranspose(oldState);
00292 }
00293 
00294 // Overridden from LinearOpBase
00295 
00296 Teuchos::RefCountPtr<const LinearOpBase<EpetraLinearOp::Scalar> >
00297 EpetraLinearOp::clone() const
00298 {
00299   assert(0); // ToDo: Implement when needed
00300   return Teuchos::null;
00301 }
00302 
00303 // Overridden from Teuchos::Describable
00304 
00305 std::string EpetraLinearOp::description() const
00306 {
00307   std::ostringstream oss;
00308   oss << "Thyra::EpetraLinearOp";
00309   oss << "{";
00310   if(op_.get()) {
00311     oss << "op=\'"<<typeid(*op_).name()<<"\'";
00312   }
00313   else {
00314     oss << "op=NULL";
00315   }
00316   oss << "}";
00317   return oss.str();
00318 }
00319 
00320 void EpetraLinearOp::describe(
00321   Teuchos::FancyOStream                &out_arg
00322   ,const Teuchos::EVerbosityLevel      verbLevel
00323   ) const
00324 {
00325   typedef Teuchos::ScalarTraits<Scalar>  ST;
00326   using Teuchos::RefCountPtr;
00327   using Teuchos::FancyOStream;
00328   using Teuchos::OSTab;
00329   using Teuchos::describe;
00330   RefCountPtr<FancyOStream> out = rcp(&out_arg,false);
00331   OSTab tab(out);
00332   switch(verbLevel) {
00333     case Teuchos::VERB_DEFAULT:
00334     case Teuchos::VERB_LOW:
00335       *out << this->description() << std::endl;
00336       break;
00337     case Teuchos::VERB_MEDIUM:
00338     case Teuchos::VERB_HIGH:
00339     case Teuchos::VERB_EXTREME:
00340     {
00341       *out
00342         << "type = \'Thyra::EpetraLinearOp\', "
00343         << "rangeDim = " << this->range()->dim() << ", domainDim = " << this->domain()->dim() << std::endl;
00344       OSTab tab(out);
00345       if(op_.get()) {
00346         *out << "op=\'"<<typeid(*op_).name()<<"\'\n";
00347         *out << "opTrans="<<toString(opTrans_)<<"\n";
00348         *out << "applyAs="<<toString(applyAs_)<<"\n";
00349         *out << "adjointSupport="<<toString(adjointSupport_)<<"\n";
00350       }
00351       else {
00352         *out << "op=NULL"<<"\n";
00353       }
00354       break;
00355     }
00356     default:
00357       TEST_FOR_EXCEPT(true); // Should never get here!
00358   }
00359 }
00360 
00361 // protected
00362 
00363 // Allocators for domain and range spaces
00364 
00365 Teuchos::RefCountPtr< const SpmdVectorSpaceBase<EpetraLinearOp::Scalar> > 
00366 EpetraLinearOp::allocateDomain(
00367   const Teuchos::RefCountPtr<Epetra_Operator>  &op 
00368   ,ETransp                                     op_trans 
00369   )  const
00370 {
00371   return create_VectorSpace(Teuchos::rcp(&op->OperatorDomainMap(),false));
00372   // ToDo: What about the transpose argument???, test this!!!
00373 }
00374 
00375 Teuchos::RefCountPtr< const SpmdVectorSpaceBase<EpetraLinearOp::Scalar> > 
00376 EpetraLinearOp::allocateRange(
00377   const Teuchos::RefCountPtr<Epetra_Operator>  &op 
00378   ,ETransp                                     op_trans 
00379   )  const
00380 {
00381   return create_VectorSpace(Teuchos::rcp(&op->OperatorRangeMap(),false));
00382   // ToDo: What about the transpose argument???, test this!!!
00383 }
00384 
00385 // private
00386 
00387 const Epetra_Map& EpetraLinearOp::getRangeMap() const
00388 {
00389   return ( applyAs_ == EPETRA_OP_APPLY_APPLY ? op_->OperatorRangeMap() : op_->OperatorDomainMap() );
00390   // ToDo: What about the transpose argument???, test this!!!
00391 }
00392 
00393 const Epetra_Map& EpetraLinearOp::getDomainMap() const
00394 {
00395   return ( applyAs_ == EPETRA_OP_APPLY_APPLY ? op_->OperatorDomainMap() : op_->OperatorRangeMap() );
00396   // ToDo: What about the transpose argument???, test this!!!
00397 }
00398 
00399 } // end namespace Thyra

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1