Teko Version of the Day
Teko_NeumannSeriesPreconditionerFactory.hpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #ifndef __Teko_NeumannSeriesPreconditionerFactory_hpp__
00048 #define __Teko_NeumannSeriesPreconditionerFactory_hpp__
00049 
00050 #include "Teko_NeumannSeriesPreconditionerFactoryDecl.hpp"
00051 
00052 #include "Thyra_DefaultPreconditioner.hpp"
00053 #include "Thyra_DefaultPreconditioner.hpp"
00054 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00055 #include "Thyra_DefaultAddedLinearOp.hpp"
00056 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00057 #include "Thyra_DefaultIdentityLinearOp.hpp"
00058 
00059 #include "Teuchos_Array.hpp"
00060 #include "Teuchos_StandardParameterEntryValidators.hpp"
00061 
00062 namespace Teko {
00063 
00064 using Teuchos::RCP;
00065 using Teuchos::rcp;
00066 
00067 static RCP<Teuchos::StringToIntegralParameterEntryValidator<Teko::DiagonalType> > scalingTypeVdtor;
00068 
00069 template <typename ScalarT>
00070 NeumannSeriesPreconditionerFactory<ScalarT>::NeumannSeriesPreconditionerFactory()
00071    : numberOfTerms_(1), scalingType_(Teko::NotDiag)
00072 {
00073 }
00074 
00076 template <typename ScalarT>
00077 bool NeumannSeriesPreconditionerFactory<ScalarT>::isCompatible(const Thyra::LinearOpSourceBase<ScalarT> &fwdOpSrc) const
00078 {
00079    return true;
00080 }
00081 
00083 template <typename ScalarT>
00084 RCP<Thyra::PreconditionerBase<ScalarT> > NeumannSeriesPreconditionerFactory<ScalarT>::createPrec() const
00085 {
00086    return rcp(new Thyra::DefaultPreconditioner<ScalarT>()); 
00087 }
00088 
00097 template <typename ScalarT>
00098 void NeumannSeriesPreconditionerFactory<ScalarT>::initializePrec(const RCP<const Thyra::LinearOpSourceBase<ScalarT> > & fwdOpSrc,
00099                     Thyra::PreconditionerBase<ScalarT> * prec,
00100                     const Thyra::ESupportSolveUse supportSolveUse) const
00101 {
00102    using Thyra::scale;
00103    using Thyra::add;
00104    using Thyra::multiply;
00105 
00106    RCP<const Thyra::LinearOpBase<ScalarT> > M; // left-preconditioner
00107    RCP<const Thyra::LinearOpBase<ScalarT> > A = fwdOpSrc->getOp();
00108    if(scalingType_!=Teko::NotDiag) {
00109       M = Teko::getInvDiagonalOp(A,scalingType_);
00110       A = Thyra::multiply(M,A); 
00111    }
00112 
00113    RCP<const Thyra::LinearOpBase<ScalarT> > id = Thyra::identity<ScalarT>(A->range()); // I
00114    RCP<const Thyra::LinearOpBase<ScalarT> > idMA = add(id,scale(-1.0,A)); // I - A
00115  
00116 
00117    RCP<const Thyra::LinearOpBase<ScalarT> > precOp;
00118    if(numberOfTerms_==1) {
00119       // no terms requested, just return identity matrix
00120       precOp = id;
00121    }
00122    else {
00123       int iters = numberOfTerms_-1;
00124       // use Horner's method to computed higher order polynomials
00125       precOp = add(scale(2.0,id),scale(-1.0,A)); // I + (I - A)
00126       for(int i=0;i<iters;i++)
00127          precOp = add(id,multiply(idMA,precOp)); // I + (I - A) * p
00128    }
00129 
00130    // multiply by the preconditioner if it exists
00131    if(M!=Teuchos::null)
00132       precOp = Thyra::multiply(precOp,M);
00133    
00134    // must first cast that to be initialized
00135    Thyra::DefaultPreconditioner<ScalarT> & dPrec = Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
00136 
00137    // this const-cast is unfortunate...needs to be fixed (larger than it seems!) ECC FIXME!
00138    dPrec.initializeUnspecified(Teuchos::rcp_const_cast<Thyra::LinearOpBase<ScalarT> >(precOp));
00139 }
00140 
00142 template <typename ScalarT>
00143 void NeumannSeriesPreconditionerFactory<ScalarT>::uninitializePrec(Thyra::PreconditionerBase<ScalarT> * prec, 
00144                       RCP<const Thyra::LinearOpSourceBase<ScalarT> > * fwdOpSrc,
00145                       Thyra::ESupportSolveUse *supportSolveUse) const
00146 {
00147    Thyra::DefaultPreconditioner<ScalarT> & dPrec = Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
00148  
00149    // wipe out any old preconditioner
00150    dPrec.uninitialize();
00151 }
00152 
00153 // for ParameterListAcceptor
00154 
00156 template <typename ScalarT>
00157 void NeumannSeriesPreconditionerFactory<ScalarT>::setParameterList(const RCP<Teuchos::ParameterList> & paramList)
00158 {
00159    TEST_FOR_EXCEPT(paramList==Teuchos::null);
00160  
00161    // check the parameters are correct
00162    paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
00163 
00164    // store the parameter list
00165    paramList_ = paramList;
00166 
00167    numberOfTerms_ = paramList_->get<int>("Number of Terms");
00168 
00169    // get the scaling type
00170    scalingType_ = Teko::NotDiag;
00171    const Teuchos::ParameterEntry * entry = paramList_->getEntryPtr("Scaling Type");
00172    if(entry!=NULL)
00173       scalingType_ = scalingTypeVdtor->getIntegralValue(*entry);
00174 }
00175 
00177 template <typename ScalarT>
00178 RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getValidParameters() const
00179 {
00180    static RCP<Teuchos::ParameterList> validPL;
00181 
00182    // only fill valid parameter list once
00183    if(validPL==Teuchos::null) {
00184       RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00185    
00186       // build the validator for scaling type
00187       scalingTypeVdtor = Teuchos::stringToIntegralParameterEntryValidator<DiagonalType>(
00188             Teuchos::tuple<std::string>("Diagonal","Lumped","AbsRowSum","None"),
00189             Teuchos::tuple<Teko::DiagonalType>(Teko::Diagonal,Teko::Lumped,Teko::AbsRowSum,Teko::NotDiag),
00190             "Scaling Type");
00191 
00192       pl->set<int>("Number of Terms",1,
00193                    "The number of terms to use in the Neumann series expansion.");
00194       pl->set("Scaling Type","None","The number of terms to use in the Neumann series expansion.",
00195               scalingTypeVdtor);
00196 
00197       validPL = pl;
00198    }
00199 
00200    return validPL;
00201 }
00202 
00204 template <typename ScalarT>
00205 RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::unsetParameterList()
00206 {
00207   Teuchos::RCP<Teuchos::ParameterList> oldList = paramList_;
00208   paramList_ = Teuchos::null;
00209   return oldList;
00210 }
00211 
00213 template <typename ScalarT>
00214 RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getParameterList() const
00215 {
00216   return paramList_;
00217 }
00218 
00220 template <typename ScalarT>
00221 RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getNonconstParameterList()
00222 {
00223    return paramList_;
00224 }
00225 
00226 template <typename ScalarT>
00227 std::string NeumannSeriesPreconditionerFactory<ScalarT>::description() const
00228 {
00229   std::ostringstream oss;
00230   oss << "Teko::NeumannSeriesPreconditionerFactory";
00231   return oss.str();
00232 }
00233 
00234 } // end namespace Teko
00235 
00236 #endif
 All Classes Files Functions Variables