Thyra_DefaultScaledAdjointLinearOp.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 // 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_SCALED_ADJOINT_LINEAR_OP_HPP
00030 #define THYRA_SCALED_ADJOINT_LINEAR_OP_HPP
00031 
00032 #include "Thyra_DefaultScaledAdjointLinearOpDecl.hpp"
00033 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00034 
00035 namespace Thyra {
00036 
00037 //Constructors/initializers/accessors
00038 
00039 template<class Scalar>
00040 void DefaultScaledAdjointLinearOp<Scalar>::initialize(
00041   const Scalar                                               &scalar
00042   ,const ETransp                                             &transp
00043   ,const Teuchos::RefCountPtr<LinearOpBase<Scalar> >         &Op
00044   )
00045 {
00046   initializeImpl(scalar,transp,Op,false);
00047 }
00048 
00049 template<class Scalar>
00050 void DefaultScaledAdjointLinearOp<Scalar>::initialize(
00051   const Scalar                                               &scalar
00052   ,const ETransp                                             &transp
00053   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &Op
00054   )
00055 {
00056   initializeImpl(scalar,transp,Op,false);
00057 }
00058 
00059 template<class Scalar>
00060 Teuchos::RefCountPtr<LinearOpBase<Scalar> >
00061 DefaultScaledAdjointLinearOp<Scalar>::getNonconstOp()
00062 {
00063   return getOpImpl().getNonconstObj();
00064 }
00065 
00066 template<class Scalar>
00067 Teuchos::RefCountPtr<const LinearOpBase<Scalar> >
00068 DefaultScaledAdjointLinearOp<Scalar>::getOp() const
00069 {
00070   return getOpImpl().getConstObj();
00071 }
00072 
00073 template<class Scalar>
00074 void DefaultScaledAdjointLinearOp<Scalar>::uninitialize()
00075 {
00076   typedef Teuchos::ScalarTraits<Scalar>  ST;
00077   origOp_.uninitialize();
00078   overallScalar_    = ST::zero();
00079   overallTransp_    = NOTRANS;
00080   allScalarETransp_ = Teuchos::null;
00081 
00082 }
00083 
00084 // Overridden from Teuchos::Describable
00085                                                 
00086 template<class Scalar>
00087 std::string DefaultScaledAdjointLinearOp<Scalar>::description() const
00088 {
00089   assertInitialized();
00090   typedef Teuchos::ScalarTraits<Scalar>  ST;
00091   std::ostringstream oss;
00092   oss << "DefaultScaledAdjointLinearOp<" << ST::name() << ">{"
00093       << overallScalar() << ","<<toString(overallTransp())<<","
00094       << origOp_.getConstObj()->description() << "}";
00095   return oss.str();
00096 }
00097 
00098 template<class Scalar>
00099 void DefaultScaledAdjointLinearOp<Scalar>::describe(
00100   Teuchos::FancyOStream                &out_arg
00101   ,const Teuchos::EVerbosityLevel      verbLevel
00102   ) const
00103 {
00104   typedef Teuchos::ScalarTraits<Scalar>  ST;
00105   using Teuchos::RefCountPtr;
00106   using Teuchos::FancyOStream;
00107   using Teuchos::OSTab;
00108   assertInitialized();
00109   RefCountPtr<FancyOStream> out = rcp(&out_arg,false);
00110   OSTab tab(out);
00111   switch(verbLevel) {
00112     case Teuchos::VERB_DEFAULT:
00113     case Teuchos::VERB_LOW:
00114       *out << this->description() << std::endl;
00115       break;
00116     case Teuchos::VERB_MEDIUM:
00117     case Teuchos::VERB_HIGH:
00118     case Teuchos::VERB_EXTREME:
00119     {
00120       *out
00121         << "DefaultScaledAdjointLinearOp<" << ST::name() << ">{"
00122         << "rangeDim=" << this->range()->dim()
00123         << ",domainDim=" << this->domain()->dim() << "}\n";
00124       OSTab tab(out);
00125       *out
00126         <<  "overallScalar="<< overallScalar() << std::endl
00127         <<  "overallTransp="<<toString(overallTransp()) << std::endl
00128         <<  "Constituent transformations:\n";
00129       for( int i = 0; i <= my_index_; ++i ) {
00130         const ScalarETransp<Scalar> &scalar_transp = (*allScalarETransp_)[my_index_-i];
00131         OSTab tab(out,i+1);
00132         if(scalar_transp.scalar != ST::one() && scalar_transp.transp != NOTRANS)
00133           *out << "scalar="<<scalar_transp.scalar<<",transp="<<toString(scalar_transp.transp)<<std::endl;
00134         else if(scalar_transp.scalar != ST::one())
00135           *out << "scalar="<<scalar_transp.scalar<<std::endl;
00136         else if( scalar_transp.transp != NOTRANS )
00137           *out << "transp="<<toString(scalar_transp.transp)<<std::endl;
00138         else
00139           *out << "no-transformation\n";
00140       }
00141       tab.incrTab(my_index_+2);
00142       *out << "origOp =\n" << Teuchos::describe(*origOp_.getConstObj(),verbLevel);
00143       break;
00144     }
00145     default:
00146       TEST_FOR_EXCEPT(true); // Should never be called!
00147   }
00148 }
00149 
00150 // Overridden from OpBase
00151 
00152 template<class Scalar>
00153 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00154 DefaultScaledAdjointLinearOp<Scalar>::range() const
00155 {
00156   assertInitialized();
00157   return ( real_trans(this->overallTransp())==NOTRANS ? this->getOrigOp()->range() : this->getOrigOp()->domain() );
00158 }
00159 
00160 template<class Scalar>
00161 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00162 DefaultScaledAdjointLinearOp<Scalar>::domain() const
00163 {
00164   assertInitialized();
00165   return ( real_trans(this->overallTransp())==NOTRANS ? this->getOrigOp()->domain() : this->getOrigOp()->range() );
00166 }
00167 
00168 template<class Scalar>
00169 Teuchos::RefCountPtr<const LinearOpBase<Scalar> >
00170 DefaultScaledAdjointLinearOp<Scalar>::clone() const
00171 {
00172   return Teuchos::null; // Not supported yet but could be
00173 }
00174 
00175 template<class Scalar>
00176 bool DefaultScaledAdjointLinearOp<Scalar>::opSupported(ETransp M_trans) const
00177 {
00178   assertInitialized();
00179   return Thyra::opSupported(*this->getOrigOp(),trans_trans(this->overallTransp(),M_trans));
00180 }
00181 
00182 // Overridden from LinearOpBase
00183 
00184 template<class Scalar>
00185 void DefaultScaledAdjointLinearOp<Scalar>::apply(
00186   const ETransp                     M_trans
00187   ,const MultiVectorBase<Scalar>    &X
00188   ,MultiVectorBase<Scalar>          *Y
00189   ,const Scalar                     alpha
00190   ,const Scalar                     beta
00191   ) const
00192 {
00193   assertInitialized();
00194   Thyra::apply(
00195     *this->getOrigOp(),trans_trans(M_trans,this->overallTransp())
00196     ,X,Y,Scalar(this->overallScalar()*alpha),beta
00197     );
00198 }
00199 
00200 // Overridden from ScaledAdointLinearOpBase
00201 
00202 template<class Scalar>
00203 Scalar DefaultScaledAdjointLinearOp<Scalar>::overallScalar() const
00204 {
00205   return overallScalar_;
00206 }
00207 
00208 template<class Scalar>
00209 ETransp DefaultScaledAdjointLinearOp<Scalar>::overallTransp() const
00210 {
00211   return overallTransp_;
00212 }
00213 
00214 template<class Scalar>
00215 Teuchos::RefCountPtr<LinearOpBase<Scalar> >
00216 DefaultScaledAdjointLinearOp<Scalar>::getNonconstOrigOp()
00217 {
00218   return origOp_.getNonconstObj();
00219 }
00220 
00221 template<class Scalar>
00222 Teuchos::RefCountPtr<const LinearOpBase<Scalar> >
00223 DefaultScaledAdjointLinearOp<Scalar>::getOrigOp() const
00224 {
00225   return origOp_.getConstObj();
00226 }
00227 
00228 // private
00229 
00230 template<class Scalar>
00231 void DefaultScaledAdjointLinearOp<Scalar>::initializeImpl(
00232   const Scalar                                               &scalar
00233   ,const ETransp                                             &transp
00234   ,const Teuchos::RefCountPtr<const LinearOpBase<Scalar> >   &Op
00235   ,const bool                                                isConst
00236   )
00237 {
00238   typedef Teuchos::ScalarTraits<Scalar>  ST;
00239 #ifdef TEUCHOS_DEBUG
00240   TEST_FOR_EXCEPTION(
00241     Op.get()==NULL, std::invalid_argument
00242     ,"DefaultScaledAdjointLinearOp<"<<ST::name()<<">::initialize(scalar,transp,Op): "
00243     "Error!, Op.get()==NULL is not allowed!"
00244     );
00245 #endif // TEUCHOS_DEBUG
00246   Teuchos::RefCountPtr<const DefaultScaledAdjointLinearOp<Scalar> >
00247     saOp = Teuchos::rcp_dynamic_cast<const DefaultScaledAdjointLinearOp<Scalar> >(Op);
00248   if(saOp.get()) {
00249     origOp_            = saOp->origOp_;
00250     overallScalar_     = saOp->overallScalar_*scalar;
00251     overallTransp_     = trans_trans(saOp->overallTransp_,transp) ;
00252     my_index_          = saOp->my_index_ + 1;
00253     allScalarETransp_  = saOp->allScalarETransp_;
00254   }
00255   else {
00256     if(isConst)
00257       origOp_.initialize(Op);
00258     else
00259       origOp_.initialize(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(Op));
00260     overallScalar_     = scalar;
00261     overallTransp_     = transp;
00262     my_index_          = 0;
00263     allScalarETransp_  = Teuchos::rcp(new allScalarETransp_t());
00264   }
00265   allScalarETransp_->push_back(ScalarETransp<Scalar>(scalar,transp));
00266 }
00267 
00268 
00269 template<class Scalar>
00270 typename DefaultScaledAdjointLinearOp<Scalar>::CNLOC
00271 DefaultScaledAdjointLinearOp<Scalar>::getOpImpl() const
00272 {
00273   assertInitialized();
00274   if( my_index_ > 0 ) {
00275     const ScalarETransp<Scalar> &scalar_transp = allScalarETransp_->at(my_index_);
00276     Teuchos::RefCountPtr<DefaultScaledAdjointLinearOp<Scalar> >
00277       Op = Teuchos::rcp(new DefaultScaledAdjointLinearOp<Scalar>());
00278     Op->origOp_            = origOp_;
00279     Op->overallScalar_     = overallScalar_/scalar_transp.scalar;
00280     Op->overallTransp_     = trans_trans(overallTransp_,scalar_transp.transp);
00281     Op->my_index_          = my_index_ - 1;
00282     Op->allScalarETransp_  = allScalarETransp_;
00283     return CNLOC(
00284       Teuchos::rcp_implicit_cast<LinearOpBase<Scalar> >(Op)
00285       );
00286   }
00287   return origOp_;
00288 }
00289 
00290 
00291 } // namespace Thyra
00292 
00293 #endif  // THYRA_SCALED_ADJOINT_LINEAR_OP_HPP

Generated on Thu Sep 18 12:32:31 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1