Thyra_LinearNonlinearSolver.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_LINEAR_NONLINEAR_SOLVER_BASE_HPP
00030 #define THYRA_LINEAR_NONLINEAR_SOLVER_BASE_HPP
00031 
00032 
00033 #include "Thyra_NonlinearSolverBase.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00036 #include "Teuchos_StandardParameterEntryValidators.hpp"
00037 #include "Teuchos_as.hpp"
00038 
00039 
00040 namespace Thyra {
00041 
00042 
00049 template <class Scalar>
00050 class LinearNonlinearSolver : public NonlinearSolverBase<Scalar> {
00051 public:
00052 
00055 
00057   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00059   RCP<Teuchos::ParameterList> getNonconstParameterList();
00061   RCP<Teuchos::ParameterList> unsetParameterList();
00063   RCP<const Teuchos::ParameterList> getParameterList() const;
00065   RCP<const Teuchos::ParameterList> getValidParameters() const;
00066 
00068 
00071 
00073   void setModel(
00074     const RCP<const ModelEvaluator<Scalar> > &model
00075     );
00077   RCP<const ModelEvaluator<Scalar> > getModel() const;
00079   SolveStatus<Scalar> solve(
00080     VectorBase<Scalar> *x,
00081     const SolveCriteria<Scalar> *solveCriteria,
00082     VectorBase<Scalar> *delta
00083     );
00085   RCP<LinearOpWithSolveBase<Scalar> > get_nonconst_W(const bool forceUpToDate);
00087   RCP<const LinearOpWithSolveBase<Scalar> > get_W() const;
00088 
00090 
00091 private:
00092 
00093   RCP<Teuchos::ParameterList> paramList_;
00094   RCP<const ModelEvaluator<Scalar> > model_;
00095   RCP<LinearOpWithSolveBase<Scalar> > J_;
00096 
00097 };
00098 
00099 
00104 template <class Scalar>
00105 RCP<LinearNonlinearSolver<Scalar> > linearNonlinearSolver()
00106 {
00107   return Teuchos::rcp(new LinearNonlinearSolver<Scalar>());
00108 }
00109 
00110 
00111 // ////////////////////////
00112 // Defintions
00113 
00114 
00115 // Overridden from Teuchos::ParameterListAcceptor
00116 
00117 
00118 template<class Scalar>
00119 void LinearNonlinearSolver<Scalar>::setParameterList(
00120   RCP<Teuchos::ParameterList> const& paramList
00121   )
00122 {
00123   using Teuchos::get;
00124   TEST_FOR_EXCEPT(is_null(paramList));
00125   paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
00126   paramList_ = paramList;
00127   // ToDo: Accept some parameters if this makes sense!
00128   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00129 #ifdef TEUCHOS_DEBUG
00130   paramList_->validateParameters(*getValidParameters(),0);
00131 #endif // TEUCHOS_DEBUG
00132 }
00133 
00134 
00135 template<class Scalar>
00136 RCP<Teuchos::ParameterList>
00137 LinearNonlinearSolver<Scalar>::getNonconstParameterList()
00138 {
00139   return paramList_;
00140 }
00141 
00142 
00143 template<class Scalar>
00144 RCP<Teuchos::ParameterList>
00145 LinearNonlinearSolver<Scalar>::unsetParameterList()
00146 {
00147   RCP<Teuchos::ParameterList> _paramList = paramList_;
00148   paramList_ = Teuchos::null;
00149   return _paramList;
00150 }
00151 
00152 
00153 template<class Scalar>
00154 RCP<const Teuchos::ParameterList>
00155 LinearNonlinearSolver<Scalar>::getParameterList() const
00156 {
00157   return paramList_;
00158 }
00159 
00160 
00161 template<class Scalar>
00162 RCP<const Teuchos::ParameterList>
00163 LinearNonlinearSolver<Scalar>::getValidParameters() const
00164 {
00165   using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
00166   static RCP<const Teuchos::ParameterList> validPL;
00167   if (is_null(validPL)) {
00168     RCP<Teuchos::ParameterList>
00169       pl = Teuchos::parameterList();
00170     // ToDo: Set up some parameters when needed!
00171     Teuchos::setupVerboseObjectSublist(&*pl);
00172     validPL = pl;
00173   }
00174   return validPL;
00175 }
00176 
00177 
00178 // Overridden from NonlinearSolverBase
00179 
00180 
00181 template <class Scalar>
00182 void LinearNonlinearSolver<Scalar>::setModel(
00183   const RCP<const ModelEvaluator<Scalar> > &model
00184   )
00185 {
00186   TEST_FOR_EXCEPT(model.get()==NULL);
00187   model_ = model;
00188   J_ = Teuchos::null;
00189 }
00190 
00191 
00192 template <class Scalar>
00193 RCP<const ModelEvaluator<Scalar> >
00194 LinearNonlinearSolver<Scalar>::getModel() const
00195 {
00196   return model_;
00197 }
00198 
00199 
00200 template <class Scalar>
00201 SolveStatus<Scalar> LinearNonlinearSolver<Scalar>::solve(
00202   VectorBase<Scalar> *x,
00203   const SolveCriteria<Scalar> *solveCriteria,
00204   VectorBase<Scalar> *delta
00205   )
00206 {
00207 
00208   using std::endl;
00209   using Teuchos::incrVerbLevel;
00210   using Teuchos::describe;
00211   using Teuchos::as;
00212   using Teuchos::rcp;
00213   using Teuchos::OSTab;
00214   using Teuchos::getFancyOStream;
00215   typedef Teuchos::ScalarTraits<Scalar> ST;
00216   typedef Thyra::ModelEvaluatorBase MEB;
00217   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00218   typedef Thyra::LinearOpWithSolveBase<Scalar> LOWSB;
00219   typedef Teuchos::VerboseObjectTempState<LOWSB> VOTSLOWSB;
00220 
00221 #ifdef TEUCHOS_DEBUG
00222   TEST_FOR_EXCEPT(0==x);
00223   THYRA_ASSERT_VEC_SPACES(
00224     "TimeStepNonlinearSolver<Scalar>::solve(...)",
00225     *x->space(),*model_->get_x_space() );
00226   TEST_FOR_EXCEPT(
00227     0!=solveCriteria && "ToDo: Support passed in solve criteria!" );
00228 #endif
00229   
00230   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00231   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00232   const bool showTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW));
00233   const bool dumpAll = (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME)); 
00234   TEUCHOS_OSTAB;
00235   VOTSME stateModel_outputTempState(model_,out,incrVerbLevel(verbLevel,-1));
00236   if(out.get() && showTrace)
00237     *out
00238       << "\nEntering LinearNonlinearSolver::solve(...) ...\n"
00239       << "\nmodel = " << describe(*model_,verbLevel);
00240 
00241   if(out.get() && dumpAll) {
00242     *out << "\nInitial guess:\n";
00243     *out << "\nx = " << *x;
00244   }
00245 
00246   // Compute the Jacobian and the residual at the input point!
00247   if(!J_.get()) J_ = model_->create_W();
00248   RCP<VectorBase<Scalar> >
00249     f = createMember(model_->get_f_space());
00250   if(out.get() && showTrace)
00251     *out << "\nEvaluating the model f and W ...\n";
00252   eval_f_W( *model_, *x,  &*f, &*J_ );
00253 
00254   // Solve the system: J*dx = -f
00255   RCP<VectorBase<Scalar> >
00256     dx = createMember(model_->get_x_space());
00257   if(out.get() && showTrace)
00258     *out << "\nSolving the system J*dx = -f ...\n";
00259   VOTSLOWSB J_outputTempState(J_,out,incrVerbLevel(verbLevel,-1));
00260   assign( &*dx, ST::zero() );
00261   Thyra::SolveStatus<Scalar>
00262     linearSolveStatus = Thyra::solve( *J_, NOTRANS, *f, &*dx );
00263   if(out.get() && showTrace)
00264     *out << "\nLinear solve status:\n" << linearSolveStatus;
00265   Vt_S( &*dx, Scalar(-ST::one()) );
00266   if(out.get() && dumpAll)
00267     *out << "\ndx = " << Teuchos::describe(*dx,verbLevel);
00268   if (delta != NULL) {
00269     Thyra::assign( delta, *dx );
00270     if(out.get() && dumpAll)
00271       *out << "\ndelta = " << Teuchos::describe(*delta,verbLevel);
00272   }
00273 
00274   // Update the solution: x += dx
00275   Vp_V( x, *dx );
00276   if(out.get() && dumpAll)
00277     *out << "\nUpdated solution x = " << Teuchos::describe(*x,verbLevel);
00278   
00279   if(out.get() && showTrace)
00280     *out << "\nLeaving LinearNonlinearSolver::solve(...) ...\n";
00281   
00282   // Return default status
00283   return SolveStatus<Scalar>();
00284 
00285 }
00286 
00287 
00288 template <class Scalar>
00289 RCP<LinearOpWithSolveBase<Scalar> >
00290 LinearNonlinearSolver<Scalar>::get_nonconst_W(const bool forceUpToDate)
00291 {
00292   if (forceUpToDate) {
00293     TEST_FOR_EXCEPT(forceUpToDate);
00294   }
00295   return J_;
00296 }
00297 
00298 
00299 template <class Scalar>
00300 RCP<const LinearOpWithSolveBase<Scalar> >
00301 LinearNonlinearSolver<Scalar>::get_W() const
00302 {
00303   return J_;
00304 }
00305 
00306 
00307 } // namespace Thyra
00308 
00309 
00310 #endif // THYRA_LINEAR_NONLINEAR_SOLVER_BASE_HPP

Generated on Wed May 12 21:42:54 2010 for Thyra Nonlinear Solver Support by  doxygen 1.4.7