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