Teko Version of the Day
Teko_SmootherPreconditionerFactory.cpp
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 // Teko includes
00048 #include "Teko_SmootherPreconditionerFactory.hpp"
00049 
00050 #include "Teko_PreconditionerInverseFactory.hpp"
00051 #include "Thyra_MultiVectorStdOps.hpp"
00052 
00053 namespace Teko {
00054 
00055 // A small utility class to help distinguish smoothers
00056 class SmootherRequestMesg : public RequestMesg {
00057 public:
00058    SmootherRequestMesg(unsigned int block)
00059       : RequestMesg("__smoother_request_message__")
00060       , block_(block) {}
00061 
00062    unsigned int getBlock() const 
00063    { return block_; }
00064 
00065 private:
00066    unsigned int block_;
00067 };
00068 
00072 
00073 SmootherLinearOp::SmootherLinearOp(const LinearOp & A, const LinearOp & invM,unsigned int applications,bool useDestAsInitialGuess)
00074    : A_(A), invM_(invM), applications_(applications), initialGuessType_(Unspecified), requestMesg_(Teuchos::null)
00075 {
00076    // set initial guess type
00077    initialGuessType_ = useDestAsInitialGuess ? DestAsInitialGuess : NoInitialGuess;
00078 }
00079 
00080 SmootherLinearOp::SmootherLinearOp(const LinearOp & A, const LinearOp & invM,unsigned int applications,unsigned int block)
00081    : A_(A), invM_(invM), applications_(applications), initialGuessType_(RequestInitialGuess), requestMesg_(Teuchos::null)
00082 {
00083    requestMesg_ = Teuchos::rcp(new SmootherRequestMesg(block));
00084 }
00085 
00086 void SmootherLinearOp::implicitApply(const MultiVector & b, MultiVector & x,
00087                                      const double alpha, const double beta) const
00088 {
00089    using Teuchos::RCP;
00090 
00091    MultiVector residual = deepcopy(b); // residual = b
00092    MultiVector scrap = deepcopy(b);    // scrap = b
00093    MultiVector error;                  // location for initial guess
00094 
00095    // construct initial guess: required to assign starting point for destination
00096    // vector appropriately
00097    switch(initialGuessType_) {
00098    case RequestInitialGuess:
00099       // get initial guess from request handler
00100       error = deepcopy(getRequestHandler()->request<MultiVector>(*requestMesg_));
00101       Thyra::assign<double>(x.ptr(),*error); // x = error (initial guess)
00102       break;
00103    case DestAsInitialGuess:
00104       error = deepcopy(x);    // error = x
00105       break;
00106    case NoInitialGuess:
00107       Thyra::assign<double>(x.ptr(),0.0); // x = 0
00108       error = deepcopy(x);                // error = x
00109       break;
00110    case Unspecified:
00111    default:
00112       TEUCHOS_ASSERT(false);
00113    }
00114 
00115    for(unsigned int current=0;current<applications_;++current) {
00116       // compute current residual
00117       Teko::applyOp(A_,error,scrap);
00118       Teko::update(-1.0,scrap,1.0,residual); // residual = residual-A*error
00119 
00120       // compute appoximate correction using residual
00121       Thyra::assign(error.ptr(),0.0); // set error to zero
00122       Teko::applyOp(invM_,residual,error);
00123 
00124       // update solution with error
00125       Teko::update(1.0,error,1.0,x); // x = x+error
00126    }
00127 }
00128 
00130 void SmootherLinearOp::setRequestHandler(const Teuchos::RCP<RequestHandler> & rh)
00131 { 
00132    Teko_DEBUG_SCOPE("SmootherLinearOp::setRequestHandler",10);
00133    requestHandler_ = rh;
00134 }
00135 
00137 Teuchos::RCP<RequestHandler> SmootherLinearOp::getRequestHandler() const
00138 { 
00139    Teko_DEBUG_SCOPE("SmootherLinearOp::getRequestHandler",10);
00140    return requestHandler_;
00141 }
00142 
00143 LinearOp buildSmootherLinearOp(const LinearOp & A,const LinearOp & invM,unsigned int applications,bool useDestAsInitialGuess)
00144 {
00145    return Teuchos::rcp(new SmootherLinearOp(A,invM,applications,useDestAsInitialGuess));
00146 }
00147 
00148 LinearOp buildSmootherLinearOp(const LinearOp & A,const LinearOp & invM,unsigned int applications,unsigned int initialGuessBlock)
00149 {
00150    return Teuchos::rcp(new SmootherLinearOp(A,invM,applications,initialGuessBlock));
00151 }
00152 
00153 
00157 
00159 SmootherPreconditionerFactory::SmootherPreconditionerFactory()
00160    : sweepCount_(0), initialGuessType_(Unspecified), initialGuessBlock_(0), precFactory_(Teuchos::null) 
00161 { }
00162 
00166 LinearOp SmootherPreconditionerFactory::buildPreconditionerOperator(LinearOp & lo,PreconditionerState & state) const
00167 {
00168    TEUCHOS_TEST_FOR_EXCEPTION(precFactory_==Teuchos::null,std::runtime_error,
00169                       "ERROR: Teko::SmootherPreconditionerFactory::buildPreconditionerOperator requires that a "
00170                    << "preconditioner factory has been set. Currently it is null!");
00171 
00172    // preconditions
00173    TEUCHOS_ASSERT(sweepCount_>0);
00174    TEUCHOS_ASSERT(initialGuessType_!=Unspecified);
00175    TEUCHOS_ASSERT(precFactory_!=Teuchos::null);
00176 
00177    // build user specified preconditioner
00178    ModifiableLinearOp & invM = state.getModifiableOp("prec");
00179    if(invM==Teuchos::null)
00180       invM = Teko::buildInverse(*precFactory_,lo);
00181    else
00182       Teko::rebuildInverse(*precFactory_,lo,invM);
00183 
00184    // conditional on initial guess type, build the smoother
00185    switch(initialGuessType_) {
00186    case RequestInitialGuess:
00187       return buildSmootherLinearOp(lo,invM,sweepCount_,initialGuessBlock_);
00188    case DestAsInitialGuess:
00189       return buildSmootherLinearOp(lo,invM,sweepCount_,true); // use an initial guess
00190    case NoInitialGuess:
00191       return buildSmootherLinearOp(lo,invM,sweepCount_,false); // no initial guess
00192    case Unspecified:
00193    default:
00194       TEUCHOS_ASSERT(false);
00195    }
00196 
00197    // should never get here
00198    TEUCHOS_ASSERT(false);
00199    return Teuchos::null;
00200 }
00201 
00205 void SmootherPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & settings)
00206 {
00207    // declare all strings used by this initialization routine
00209 
00210    const std::string str_sweepCount = "Sweep Count";
00211    const std::string str_initialGuessBlock = "Initial Guess Block";
00212    const std::string str_destAsInitialGuess = "Destination As Initial Guess";
00213    const std::string str_precType = "Preconditioner Type";
00214 
00215    // default parameters
00217 
00218    initialGuessType_ = Unspecified;
00219    initialGuessBlock_ = 0;
00220    sweepCount_ = 0;
00221    precFactory_ = Teuchos::null;
00222 
00223    // get sweep counts
00225 
00226    if(settings.isParameter(str_sweepCount))
00227       sweepCount_ = settings.get<int>(str_sweepCount);
00228 
00229    // get initial guess
00231 
00232    if(settings.isParameter(str_initialGuessBlock)) {
00233       initialGuessBlock_ = settings.get<int>(str_initialGuessBlock);
00234       initialGuessType_ = RequestInitialGuess;
00235    }
00236 
00237    if(settings.isParameter(str_destAsInitialGuess)) {
00238       bool useDest = settings.get<bool>(str_destAsInitialGuess);
00239       if(useDest) {
00240          TEUCHOS_TEST_FOR_EXCEPTION(initialGuessType_!=Unspecified, std::runtime_error,
00241                             "Cannot set both \"" << str_initialGuessBlock  <<  
00242                             "\" and \""          << str_destAsInitialGuess << "\"");
00243 
00244          initialGuessType_ = DestAsInitialGuess;
00245       }
00246    }
00247 
00248    // safe to assume if the other values are not turned on there is no initial guess
00249    if(initialGuessType_==Unspecified) 
00250       initialGuessType_ = NoInitialGuess;
00251 
00252    // get preconditioner factory
00254  
00255    TEUCHOS_TEST_FOR_EXCEPTION(not settings.isParameter(str_precType),std::runtime_error,
00256                       "Parameter \"" << str_precType << "\" is required by a Teko::SmootherPreconditionerFactory");
00257       
00258    // grab library and preconditioner name
00259    Teuchos::RCP<const InverseLibrary> il = getInverseLibrary();
00260    std::string precName = settings.get<std::string>(str_precType);
00261 
00262    // build preconditioner factory
00263    precFactory_ = il->getInverseFactory(precName);
00264    TEUCHOS_TEST_FOR_EXCEPTION(precFactory_==Teuchos::null,std::runtime_error,
00265                       "ERROR: \"" << str_precType << "\" = " << precName 
00266                    << " could not be found");
00267 
00268    // post conditions required to be satisfied
00270 
00271    TEUCHOS_ASSERT(sweepCount_>0);
00272    TEUCHOS_ASSERT(initialGuessType_!=Unspecified);
00273    TEUCHOS_ASSERT(precFactory_!=Teuchos::null);
00274 }
00275 
00276 } // end namespace Teko
 All Classes Files Functions Variables