Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_EpetraLinearOp.cpp
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 // 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 #include "Thyra_EpetraLinearOp.hpp"
00043 #include "Thyra_EpetraThyraWrappers.hpp"
00044 #include "Thyra_SpmdMultiVectorBase.hpp"
00045 #include "Thyra_MultiVectorStdOps.hpp"
00046 #include "Thyra_AssertOp.hpp"
00047 #include "Teuchos_dyn_cast.hpp"
00048 #include "Teuchos_Assert.hpp"
00049 #include "Teuchos_getConst.hpp"
00050 #include "Teuchos_as.hpp"
00051 #include "Teuchos_TimeMonitor.hpp"
00052 
00053 #include "Epetra_Map.h"
00054 #include "Epetra_Vector.h"
00055 #include "Epetra_Operator.h"
00056 #include "Epetra_CrsMatrix.h" // Printing only!
00057 
00058 
00059 namespace Thyra {
00060 
00061 
00062 // Constructors / initializers / accessors
00063 
00064 
00065 EpetraLinearOp::EpetraLinearOp()
00066   :isFullyInitialized_(false),
00067    opTrans_(NOTRANS),
00068    applyAs_(EPETRA_OP_APPLY_APPLY),
00069    adjointSupport_(EPETRA_OP_ADJOINT_UNSUPPORTED)
00070 {}
00071 
00072 
00073 void EpetraLinearOp::initialize(
00074   const RCP<Epetra_Operator> &op,
00075   EOpTransp opTrans,
00076   EApplyEpetraOpAs applyAs,
00077   EAdjointEpetraOp adjointSupport,
00078   const RCP< const VectorSpaceBase<double> > &range_in,
00079   const RCP< const VectorSpaceBase<double> > &domain_in
00080   )
00081 {
00082   
00083   using Teuchos::rcp_dynamic_cast;
00084   typedef SpmdVectorSpaceBase<double> SPMDVSB;
00085   typedef ScalarProdVectorSpaceBase<double> SPVSB;
00086 
00087   // Validate input, allocate spaces, validate ...
00088 #ifdef TEUCHOS_DEBUG
00089   TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
00090     "Thyra::EpetraLinearOp::initialize(...): Error!" );
00091   // ToDo: Validate spmdRange, spmdDomain against op maps!
00092 #endif
00093 
00094   RCP<const SPMDVSB> l_spmdRange;
00095   if(!is_null(range_in))
00096     l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,true);
00097   else
00098     l_spmdRange = ( applyAs==EPETRA_OP_APPLY_APPLY
00099       ? allocateRange(op,opTrans) : allocateDomain(op,opTrans) );
00100 
00101   RCP<const SPMDVSB> l_spmdDomain;
00102   if(!is_null(domain_in))
00103     l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,true);
00104   else
00105     l_spmdDomain = ( applyAs==EPETRA_OP_APPLY_APPLY
00106       ? allocateDomain(op,opTrans) : allocateRange(op,opTrans) );
00107   
00108   // Set data (no exceptions should be thrown now)
00109   isFullyInitialized_ = true;
00110   op_ = op;
00111   rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
00112   opTrans_ = opTrans;
00113   applyAs_ = applyAs;
00114   adjointSupport_ = adjointSupport;
00115   range_ = l_spmdRange;
00116   domain_ = l_spmdDomain;
00117 
00118 }
00119 
00120 
00121 void EpetraLinearOp::partiallyInitialize(
00122   const RCP<const VectorSpaceBase<double> > &range_in,
00123   const RCP<const VectorSpaceBase<double> > &domain_in,
00124   const RCP<Epetra_Operator> &op,
00125   EOpTransp opTrans,
00126   EApplyEpetraOpAs applyAs,
00127   EAdjointEpetraOp adjointSupport
00128   )
00129 {
00130   
00131   using Teuchos::rcp_dynamic_cast;
00132   typedef SpmdVectorSpaceBase<double> SPMDVSB;
00133   typedef ScalarProdVectorSpaceBase<double> SPVSB;
00134 
00135   // Validate input, allocate spaces, validate ...
00136 #ifdef TEUCHOS_DEBUG
00137   TEUCHOS_TEST_FOR_EXCEPTION( is_null(range_in), std::invalid_argument,
00138     "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
00139   TEUCHOS_TEST_FOR_EXCEPTION( is_null(domain_in), std::invalid_argument,
00140     "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
00141   TEUCHOS_TEST_FOR_EXCEPTION( is_null(op), std::invalid_argument,
00142     "Thyra::EpetraLinearOp::partiallyInitialize(...): Error!" );
00143 #endif
00144 
00145   RCP<const SPMDVSB>
00146     l_spmdRange = rcp_dynamic_cast<const SPMDVSB>(range_in,true);
00147   RCP<const SPMDVSB>
00148     l_spmdDomain = rcp_dynamic_cast<const SPMDVSB>(domain_in,true);
00149   
00150   // Set data (no exceptions should be thrown now)
00151   isFullyInitialized_ = false;
00152   op_ = op;
00153   rowMatrix_ = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(op_);
00154   opTrans_ = opTrans;
00155   applyAs_ = applyAs;
00156   adjointSupport_ = adjointSupport;
00157   range_ = l_spmdRange;
00158   domain_ = l_spmdDomain;
00159   
00160 }
00161 
00162 
00163 void EpetraLinearOp::setFullyInitialized(bool isFullyInitialized)
00164 {
00165   // ToDo: Validate that everything matches up!
00166   isFullyInitialized_ = true;
00167 }
00168 
00169 
00170 void EpetraLinearOp::uninitialize(
00171   RCP<Epetra_Operator> *op,
00172   EOpTransp *opTrans,
00173   EApplyEpetraOpAs *applyAs,
00174   EAdjointEpetraOp *adjointSupport,
00175   RCP<const VectorSpaceBase<double> > *range_out,
00176   RCP<const VectorSpaceBase<double> > *domain_out
00177   )
00178 {
00179 
00180   if(op) *op = op_;
00181   if(opTrans) *opTrans = opTrans_;
00182   if(applyAs) *applyAs = applyAs_;
00183   if(adjointSupport) *adjointSupport = adjointSupport_;
00184   if(range_out) *range_out = range_;
00185   if(domain_out) *domain_out = domain_;
00186 
00187   isFullyInitialized_ = false;
00188   op_ = Teuchos::null;
00189   rowMatrix_ = Teuchos::null;
00190   opTrans_ = NOTRANS;
00191   applyAs_ = EPETRA_OP_APPLY_APPLY;
00192   adjointSupport_ = EPETRA_OP_ADJOINT_SUPPORTED;
00193   range_ = Teuchos::null;
00194   domain_ = Teuchos::null;
00195 
00196 }
00197 
00198 
00199 RCP< const SpmdVectorSpaceBase<double> >
00200 EpetraLinearOp::spmdRange() const
00201 {
00202   return range_;
00203 }
00204 
00205 
00206 RCP< const SpmdVectorSpaceBase<double> >
00207 EpetraLinearOp::spmdDomain() const
00208 {
00209   return domain_;
00210 }
00211 
00212 
00213 RCP<Epetra_Operator>
00214 EpetraLinearOp::epetra_op() 
00215 {
00216   return op_;
00217 }
00218 
00219 
00220 RCP<const Epetra_Operator>
00221 EpetraLinearOp::epetra_op() const 
00222 {
00223   return op_;
00224 }
00225 
00226 
00227 // Overridden from EpetraLinearOpBase
00228 
00229 
00230 void EpetraLinearOp::getNonconstEpetraOpView(
00231   const Ptr<RCP<Epetra_Operator> > &epetraOp,
00232   const Ptr<EOpTransp> &epetraOpTransp,
00233   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00234   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00235   )
00236 {
00237   *epetraOp = op_;
00238   *epetraOpTransp = opTrans_;
00239   *epetraOpApplyAs = applyAs_;
00240   *epetraOpAdjointSupport = adjointSupport_;
00241 }
00242 
00243 
00244 void EpetraLinearOp::getEpetraOpView(
00245   const Ptr<RCP<const Epetra_Operator> > &epetraOp,
00246   const Ptr<EOpTransp> &epetraOpTransp,
00247   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00248   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00249   ) const
00250 {
00251   *epetraOp = op_;
00252   *epetraOpTransp = opTrans_;
00253   *epetraOpApplyAs = applyAs_;
00254   *epetraOpAdjointSupport = adjointSupport_;
00255 }
00256 
00257 
00258 // Overridden from LinearOpBase
00259 
00260 
00261 RCP<const VectorSpaceBase<double> >
00262 EpetraLinearOp::range() const
00263 {
00264   return range_;
00265 }
00266 
00267 
00268 RCP<const VectorSpaceBase<double> >
00269 EpetraLinearOp::domain() const
00270 {
00271   return domain_;
00272 }
00273 
00274 
00275 RCP<const LinearOpBase<double> >
00276 EpetraLinearOp::clone() const
00277 {
00278   assert(0); // ToDo: Implement when needed
00279   return Teuchos::null;
00280 }
00281 
00282 
00283 // Overridden from Teuchos::Describable
00284 
00285 
00286 std::string EpetraLinearOp::description() const
00287 {
00288   std::ostringstream oss;
00289   oss << Teuchos::Describable::description() << "{";
00290   if(op_.get()) {
00291     oss << "op=\'"<<typeName(*op_)<<"\'";
00292     oss << ",rangeDim="<<this->range()->dim();
00293     oss << ",domainDim="<<this->domain()->dim();
00294   }
00295   else {
00296     oss << "op=NULL";
00297   }
00298   oss << "}";
00299   return oss.str();
00300 }
00301 
00302 
00303 void EpetraLinearOp::describe(
00304   FancyOStream &out,
00305   const Teuchos::EVerbosityLevel verbLevel
00306   ) const
00307 {
00308   typedef Teuchos::ScalarTraits<double> ST;
00309   using Teuchos::includesVerbLevel;
00310   using Teuchos::as;
00311   using Teuchos::rcp_dynamic_cast;
00312   using Teuchos::OSTab;
00313   using Teuchos::describe;
00314   OSTab tab(out);
00315   if ( as<int>(verbLevel) == as<int>(Teuchos::VERB_LOW) || is_null(op_)) {
00316     out << this->description() << std::endl;
00317   }
00318   else if (includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
00319     out
00320       << Teuchos::Describable::description()
00321       << "{"
00322       << "rangeDim=" << this->range()->dim()
00323       << ",domainDim=" << this->domain()->dim()
00324       << "}\n";
00325     OSTab tab2(out);
00326     if (op_.get()) {
00327       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00328         out << "opTrans="<<toString(opTrans_)<<"\n";
00329         out << "applyAs="<<toString(applyAs_)<<"\n";
00330         out << "adjointSupport="<<toString(adjointSupport_)<<"\n";
00331         out << "op="<<typeName(*op_)<<"\n";
00332       }
00333       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00334         OSTab tab3(out);
00335         RCP<const Epetra_CrsMatrix>
00336           csr_op = rcp_dynamic_cast<const Epetra_CrsMatrix>(op_);
00337         if (!is_null(csr_op)) {
00338           csr_op->Print(out);
00339         }
00340       }
00341     }
00342     else {
00343       out << "op=NULL"<<"\n";
00344     }
00345   }
00346 }
00347 
00348 
00349 // protected
00350 
00351 
00352 // Protected member functions overridden from LinearOpBase
00353 
00354 
00355 bool EpetraLinearOp::opSupportedImpl(EOpTransp M_trans) const
00356 {
00357   if (!isFullyInitialized_)
00358     return false;
00359   return ( M_trans == NOTRANS
00360     ? true : adjointSupport_==EPETRA_OP_ADJOINT_SUPPORTED );
00361 }
00362 
00363 
00364 void EpetraLinearOp::applyImpl(
00365   const EOpTransp M_trans,
00366   const MultiVectorBase<double> &X_in,
00367   const Ptr<MultiVectorBase<double> > &Y_inout,
00368   const double alpha,
00369   const double beta
00370   ) const
00371 {
00372 
00373   THYRA_FUNC_TIME_MONITOR("Thyra::EpetraLinearOp::euclideanApply");
00374 
00375   const EOpTransp real_M_trans = real_trans(M_trans);
00376 
00377 #ifdef TEUCHOS_DEBUG
00378   TEUCHOS_TEST_FOR_EXCEPT(!isFullyInitialized_);
00379   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00380     "EpetraLinearOp::euclideanApply(...)", *this, M_trans, X_in, Y_inout
00381     );
00382   TEUCHOS_TEST_FOR_EXCEPTION(
00383     real_M_trans==TRANS && adjointSupport_==EPETRA_OP_ADJOINT_UNSUPPORTED,
00384     Exceptions::OpNotSupported,
00385     "EpetraLinearOp::apply(...): *this was informed that adjoints "
00386     "are not supported when initialized." 
00387     );
00388 #endif
00389 
00390   const RCP<const VectorSpaceBase<double> > XY_domain = X_in.domain();
00391   const int numCols = XY_domain->dim();
00392  
00393   //
00394   // Get Epetra_MultiVector objects for the arguments
00395   //
00396   // 2007/08/18: rabartl: Note: After profiling, I found that calling the more
00397   // general functions get_Epetra_MultiVector(...) was too slow. These
00398   // functions must ensure that memory is being remembered efficiently and the
00399   // use of extra data with the RCP and other things is slow.
00400   //
00401   RCP<const Epetra_MultiVector> X;
00402   RCP<Epetra_MultiVector> Y;
00403   {
00404     THYRA_FUNC_TIME_MONITOR_DIFF(
00405       "Thyra::EpetraLinearOp::euclideanApply: Convert MultiVectors", MultiVectors);
00406     // X
00407     X = get_Epetra_MultiVector(
00408       real_M_trans==NOTRANS ? getDomainMap() : getRangeMap(), X_in );
00409     // Y
00410     if( beta == 0 ) {
00411       Y = get_Epetra_MultiVector(
00412         real_M_trans==NOTRANS ? getRangeMap() : getDomainMap(), *Y_inout );
00413     }
00414   }
00415 
00416   //
00417   // Set the operator mode
00418   //
00419 
00420   /* We need to save the transpose state here, and then reset it after 
00421    * application. The reason for this is that if we later apply the 
00422    * operator outside Thyra (in Aztec, for instance), it will remember
00423    * the transpose flag set here. */
00424   bool oldState = op_->UseTranspose();
00425   op_->SetUseTranspose(
00426     real_trans(trans_trans(opTrans_,M_trans)) == NOTRANS ? false : true );
00427 
00428   //
00429   // Perform the apply operation
00430   //
00431   {
00432     THYRA_FUNC_TIME_MONITOR_DIFF(
00433       "Thyra::EpetraLinearOp::euclideanApply: Apply", Apply);
00434     if( beta == 0.0 ) {
00435       // Y = M * X
00436       if( applyAs_ == EPETRA_OP_APPLY_APPLY ) {
00437         THYRA_FUNC_TIME_MONITOR_DIFF(
00438           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Apply",
00439           ApplyApply);
00440         op_->Apply( *X, *Y );
00441       }
00442       else if( applyAs_ == EPETRA_OP_APPLY_APPLY_INVERSE ) {
00443         THYRA_FUNC_TIME_MONITOR_DIFF(
00444           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): ApplyInverse",
00445           ApplyApplyInverse);
00446         op_->ApplyInverse( *X, *Y );
00447       }
00448       else {
00449 #ifdef TEUCHOS_DEBUG
00450         TEUCHOS_TEST_FOR_EXCEPT(true);
00451 #endif
00452       }
00453       // Y = alpha * Y
00454       if( alpha != 1.0 ) {
00455         THYRA_FUNC_TIME_MONITOR_DIFF(
00456           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Scale Y",
00457           Scale);
00458         Y->Scale(alpha);
00459       }
00460     }
00461     else {  // beta != 0.0
00462       // Y_inout = beta * Y_inout
00463       if(beta != 0.0) {
00464         THYRA_FUNC_TIME_MONITOR_DIFF(
00465           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Scale Y",
00466           Scale);
00467         scale( beta, Y_inout );
00468       }
00469       else {
00470         THYRA_FUNC_TIME_MONITOR_DIFF(
00471           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Y=0",
00472           Apply2);
00473         assign( Y_inout, 0.0 );
00474       }
00475       // T = M * X
00476       Epetra_MultiVector T(op_->OperatorRangeMap(), numCols, false);
00477       // NOTE: Above, op_->OperatorRange() will be right for either
00478       // non-transpose or transpose because we have already set the
00479       // UseTranspose flag correctly.
00480       if( applyAs_ == EPETRA_OP_APPLY_APPLY ) {
00481         THYRA_FUNC_TIME_MONITOR_DIFF(
00482           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Apply",
00483           Apply2);
00484         op_->Apply( *X, T );
00485       }
00486       else if( applyAs_ == EPETRA_OP_APPLY_APPLY_INVERSE ) {
00487         THYRA_FUNC_TIME_MONITOR_DIFF(
00488           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): ApplyInverse",
00489           ApplyInverse);
00490         op_->ApplyInverse( *X, T );
00491       }
00492       else {
00493 #ifdef TEUCHOS_DEBUG
00494         TEUCHOS_TEST_FOR_EXCEPT(true);
00495 #endif
00496       }
00497       // Y_inout += alpha * T
00498       {
00499         THYRA_FUNC_TIME_MONITOR_DIFF(
00500           "Thyra::EpetraLinearOp::euclideanApply: Apply(beta!=0): Update Y",
00501           Update);
00502         update(
00503           alpha,
00504           *create_MultiVector(
00505             Teuchos::rcp(&Teuchos::getConst(T),false),
00506             Y_inout->range(),
00507             XY_domain
00508             ),
00509           Y_inout
00510           );
00511       }
00512     }
00513   }
00514 
00515   // Reset the transpose state
00516   op_->SetUseTranspose(oldState);
00517 
00518   // 2009/04/14: ToDo: This will not reset the transpose flag correctly if an
00519   // exception is thrown!
00520 
00521 }
00522 
00523 
00524 // Protected member functions overridden from ScaledLinearOpBase
00525 
00526 
00527 bool EpetraLinearOp::supportsScaleLeftImpl() const
00528 {
00529   return nonnull(rowMatrix_);
00530 }
00531 
00532 
00533 bool EpetraLinearOp::supportsScaleRightImpl() const
00534 {
00535   return nonnull(rowMatrix_);
00536 }
00537 
00538 
00539 void EpetraLinearOp::scaleLeftImpl(const VectorBase<double> &row_scaling_in)
00540 {
00541   using Teuchos::rcpFromRef;
00542   const RCP<const Epetra_Vector> row_scaling =
00543     get_Epetra_Vector(getRangeMap(), rcpFromRef(row_scaling_in));
00544   rowMatrix_->LeftScale(*row_scaling);
00545 }
00546 
00547 
00548 void EpetraLinearOp::scaleRightImpl(const VectorBase<double> &col_scaling_in)
00549 {
00550   using Teuchos::rcpFromRef;
00551   const RCP<const Epetra_Vector> col_scaling =
00552     get_Epetra_Vector(getDomainMap(), rcpFromRef(col_scaling_in));
00553   rowMatrix_->RightScale(*col_scaling);
00554 }
00555 
00556 
00557 // Protected member functions overridden from RowStatLinearOpBase
00558 
00559 
00560 bool EpetraLinearOp::rowStatIsSupportedImpl(
00561   const RowStatLinearOpBaseUtils::ERowStat rowStat) const
00562 {
00563   if (is_null(rowMatrix_)) {
00564     return false;
00565   }
00566   switch (rowStat) {
00567     case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
00568       return true;
00569     default:
00570       TEUCHOS_TEST_FOR_EXCEPT(true);
00571   }
00572   return false; // Will never be called!
00573 }
00574 
00575 
00576 void EpetraLinearOp::getRowStatImpl(
00577   const RowStatLinearOpBaseUtils::ERowStat rowStat,
00578   const Ptr<VectorBase<double> > &rowStatVec_in
00579   ) const
00580 {
00581   using Teuchos::rcpFromPtr;
00582   const RCP<Epetra_Vector> rowStatVec =
00583     get_Epetra_Vector(getRangeMap(), rcpFromPtr(rowStatVec_in));
00584   switch (rowStat) {
00585     case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
00586       rowMatrix_->InvRowSums(*rowStatVec);
00587       break;
00588     default:
00589       TEUCHOS_TEST_FOR_EXCEPT(true);
00590   }
00591 }
00592 
00593 
00594 // Allocators for domain and range spaces
00595 
00596 
00597 RCP< const SpmdVectorSpaceBase<double> > 
00598 EpetraLinearOp::allocateDomain(
00599   const RCP<Epetra_Operator> &op,
00600   EOpTransp op_trans
00601   ) const
00602 {
00603   return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00604     create_VectorSpace(Teuchos::rcp(&op->OperatorDomainMap(),false))
00605     );
00606   // ToDo: What about the transpose argument???, test this!!!
00607 }
00608 
00609 
00610 RCP<const SpmdVectorSpaceBase<double> > 
00611 EpetraLinearOp::allocateRange(
00612   const RCP<Epetra_Operator> &op,
00613   EOpTransp op_trans
00614   ) const
00615 {
00616   return Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00617     create_VectorSpace(Teuchos::rcp(&op->OperatorRangeMap(),false))
00618     );
00619   // ToDo: What about the transpose argument???, test this!!!
00620 }
00621 
00622 
00623 // private
00624 
00625 
00626 const Epetra_Map& EpetraLinearOp::getRangeMap() const
00627 {
00628   return ( applyAs_ == EPETRA_OP_APPLY_APPLY
00629     ? op_->OperatorRangeMap() : op_->OperatorDomainMap() );
00630   // ToDo: What about the transpose argument???, test this!!!
00631 }
00632 
00633 
00634 const Epetra_Map& EpetraLinearOp::getDomainMap() const
00635 {
00636   return ( applyAs_ == EPETRA_OP_APPLY_APPLY
00637     ? op_->OperatorDomainMap() : op_->OperatorRangeMap() );
00638   // ToDo: What about the transpose argument???, test this!!!
00639 }
00640 
00641 
00642 } // end namespace Thyra
00643 
00644 
00645 // Nonmembers
00646 
00647 
00648 Teuchos::RCP<Thyra::EpetraLinearOp>
00649 Thyra::nonconstEpetraLinearOp()
00650 {
00651   return Teuchos::rcp(new EpetraLinearOp());
00652 }
00653 
00654 
00655 Teuchos::RCP<Thyra::EpetraLinearOp>
00656 Thyra::partialNonconstEpetraLinearOp(
00657   const RCP<const VectorSpaceBase<double> > &range,
00658   const RCP<const VectorSpaceBase<double> > &domain,
00659   const RCP<Epetra_Operator> &op,
00660   EOpTransp opTrans,
00661   EApplyEpetraOpAs applyAs,
00662   EAdjointEpetraOp adjointSupport
00663   )
00664 {
00665   RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
00666   thyraEpetraOp->partiallyInitialize(
00667     range, domain,op,opTrans, applyAs, adjointSupport
00668     );
00669   return thyraEpetraOp;
00670 }
00671 
00672 
00673 Teuchos::RCP<Thyra::EpetraLinearOp>
00674 Thyra::nonconstEpetraLinearOp(
00675   const RCP<Epetra_Operator> &op,
00676   EOpTransp opTrans,
00677   EApplyEpetraOpAs applyAs,
00678   EAdjointEpetraOp adjointSupport,
00679   const RCP< const VectorSpaceBase<double> > &range,
00680   const RCP< const VectorSpaceBase<double> > &domain
00681   )
00682 {
00683   RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
00684   thyraEpetraOp->initialize(
00685     op,opTrans, applyAs, adjointSupport, range, domain
00686     );
00687   return thyraEpetraOp;
00688 }
00689 
00690 
00691 Teuchos::RCP<const Thyra::EpetraLinearOp>
00692 Thyra::epetraLinearOp(
00693   const RCP<const Epetra_Operator> &op,
00694   EOpTransp opTrans,
00695   EApplyEpetraOpAs applyAs,
00696   EAdjointEpetraOp adjointSupport,
00697   const RCP<const VectorSpaceBase<double> > &range,
00698   const RCP<const VectorSpaceBase<double> > &domain
00699   )
00700 {
00701   RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
00702   thyraEpetraOp->initialize(
00703     Teuchos::rcp_const_cast<Epetra_Operator>(op), // Safe cast due to return type!
00704     opTrans, applyAs, adjointSupport, range, domain
00705     );
00706   return thyraEpetraOp;
00707 }
00708 
00709 
00710 Teuchos::RCP<Thyra::EpetraLinearOp>
00711 Thyra::nonconstEpetraLinearOp(
00712   const RCP<Epetra_Operator> &op,
00713   const std::string &label,
00714   EOpTransp opTrans,
00715   EApplyEpetraOpAs applyAs,
00716   EAdjointEpetraOp adjointSupport,
00717   const RCP<const VectorSpaceBase<double> > &range,
00718   const RCP<const VectorSpaceBase<double> > &domain
00719   )
00720 {
00721   RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
00722   thyraEpetraOp->initialize(
00723     op,opTrans, applyAs, adjointSupport, range, domain
00724     );
00725   thyraEpetraOp->setObjectLabel(label);
00726   return thyraEpetraOp;
00727 }
00728 
00729 
00730 Teuchos::RCP<const Thyra::EpetraLinearOp>
00731 Thyra::epetraLinearOp(
00732   const RCP<const Epetra_Operator> &op,
00733   const std::string &label,
00734   EOpTransp opTrans,
00735   EApplyEpetraOpAs applyAs,
00736   EAdjointEpetraOp adjointSupport,
00737   const RCP< const SpmdVectorSpaceBase<double> > &range,
00738   const RCP< const SpmdVectorSpaceBase<double> > &domain
00739   )
00740 {
00741   RCP<EpetraLinearOp> thyraEpetraOp = Teuchos::rcp(new EpetraLinearOp());
00742   thyraEpetraOp->initialize(
00743     Teuchos::rcp_const_cast<Epetra_Operator>(op), // Safe cast due to return type!
00744     opTrans, applyAs, adjointSupport, range, domain
00745     );
00746   thyraEpetraOp->setObjectLabel(label);
00747   return thyraEpetraOp;
00748 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines