00001
00002
00003
00004
00005 #ifndef PLAYA_NONLINEAROPERATORBASE_HPP
00006 #define PLAYA_NONLINEAROPERATORBASE_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaHandleable.hpp"
00010 #include "PlayaOut.hpp"
00011 #include "PlayaObjectWithVerbosity.hpp"
00012 #include "PlayaVectorDecl.hpp"
00013 #include "PlayaLinearOperatorDecl.hpp"
00014 #include "PlayaLinearCombinationDecl.hpp"
00015
00016 namespace Playa
00017 {
00018 using namespace Teuchos;
00019
00020
00021
00022
00023 template <class Scalar>
00024 class NonlinearOperatorBase
00025 : public Handleable<NonlinearOperatorBase<Scalar> >,
00026 public ObjectWithVerbosity
00027 {
00028 public:
00029
00030
00031 NonlinearOperatorBase()
00032 : domain_(), range_(),
00033 jacobianIsValid_(false),
00034 residualIsValid_(false),
00035 currentEvalPt_(),
00036 currentFunctionValue_(),
00037 currentJ_()
00038 {;}
00039
00040
00041 NonlinearOperatorBase(const VectorSpace<Scalar>& domain,
00042 const VectorSpace<Scalar>& range)
00043 : domain_(domain.ptr()), range_(range.ptr()),
00044 jacobianIsValid_(false),
00045 residualIsValid_(false),
00046 currentEvalPt_(),
00047 currentFunctionValue_(),
00048 currentJ_()
00049 {;}
00050
00051
00052 const RCP<const VectorSpaceBase<Scalar> >& domain() const
00053 {return domain_;}
00054
00055
00056 const RCP<const VectorSpaceBase<Scalar> >& range() const
00057 {return range_;}
00058
00059
00060 void setEvalPt(const Vector<Scalar>& x) const
00061 {
00062 if (this->verb() >= 1)
00063 {
00064 Out::os() << "NonlinearOperatorBase Setting new eval pt";
00065 if (this->verb() > 3)
00066 {
00067 Out::os() << " to " << std::endl ;
00068 x.print(Out::os());
00069 }
00070 Out::os() << std::endl;
00071 }
00072 jacobianIsValid_ = false;
00073 residualIsValid_ = false;
00074
00075 currentEvalPt_.acceptCopyOf(x);
00076 }
00077
00078
00079
00080 const Vector<double>& currentEvalPt() const {return currentEvalPt_;}
00081
00082
00083 LinearOperator<double> getJacobian() const
00084 {
00085 if (this->verb() > 1)
00086 {
00087 Out::os() << "NonlinearOperatorBase getting Jacobian" << std::endl;
00088 }
00089 if (!jacobianIsValid_)
00090 {
00091 if (this->verb() > 3)
00092 {
00093 Out::os() << "...computing new J and F" << std::endl;
00094 }
00095 currentJ_
00096 = computeJacobianAndFunction(currentFunctionValue_);
00097 jacobianIsValid_ = true;
00098 residualIsValid_ = true;
00099 }
00100 else
00101 {
00102 if (this->verb() > 1)
00103 {
00104 Out::os() << "...reusing valid J" << std::endl;
00105 }
00106 }
00107 if (this->verb() > 3)
00108 {
00109 Out::os() << "J is " << std::endl;
00110 currentJ_.print(Out::os());
00111 Out::os() << std::endl;
00112 }
00113 return currentJ_;
00114 }
00115
00116
00117
00118
00119 Vector<double> getFunctionValue() const
00120 {
00121 if (this->verb() > 1)
00122 {
00123 Out::os() << "NonlinearOperatorBase getting function value" << std::endl;
00124 }
00125 if (!residualIsValid_)
00126 {
00127 if (this->verb() > 1)
00128 {
00129 Out::os() << "...computing new F" << std::endl;
00130 }
00131 currentFunctionValue_ = computeFunctionValue();
00132 residualIsValid_ = true;
00133 }
00134 else
00135 {
00136 if (this->verb() > 1)
00137 {
00138 Out::os() << "...reusing valid F" << std::endl;
00139 }
00140 }
00141
00142 if (this->verb() > 3)
00143 {
00144 Out::os() << "F is " << std::endl;
00145 currentFunctionValue_.print(Out::os());
00146 Out::os() << std::endl;
00147 }
00148 return currentFunctionValue_;
00149 }
00150
00151
00152
00153 virtual Vector<double> getInitialGuess() const = 0 ;
00154
00155
00156 protected:
00157
00158
00159 virtual LinearOperator<Scalar> computeJacobianAndFunction(Vector<double>& functionValue) const = 0 ;
00160
00161
00162 virtual Vector<Scalar> computeFunctionValue() const
00163 {
00164 computeJacobianAndFunction(currentFunctionValue_);
00165 return currentFunctionValue_;
00166 }
00167
00168
00169
00170
00171 void setDomainAndRange(const VectorSpace<Scalar>& domain,
00172 const VectorSpace<Scalar>& range)
00173 {
00174 domain_ = domain.ptr();
00175 range_ = range.ptr();
00176 }
00177
00178
00179 private:
00180
00181 RCP<const VectorSpaceBase<Scalar> > domain_;
00182
00183
00184 RCP<const VectorSpaceBase<Scalar> > range_;
00185
00186
00187 mutable bool jacobianIsValid_;
00188
00189
00190 mutable bool residualIsValid_;
00191
00192
00193 mutable Vector<double> currentEvalPt_;
00194
00195
00196 mutable Vector<double> currentFunctionValue_;
00197
00198
00199 mutable LinearOperator<double> currentJ_;
00200 };
00201
00202
00203
00204
00205 }
00206
00207
00208 #endif