00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "NOX_Common.H"
00043
00044 #include "NOX_Direction_Broyden.H"
00045 #include "NOX_Abstract_Vector.H"
00046 #include "NOX_Abstract_Group.H"
00047 #include "NOX_Solver_Generic.H"
00048 #include "NOX_Solver_LineSearchBased.H"
00049 #include "NOX_Utils.H"
00050 #include "NOX_GlobalData.H"
00051
00052
00053
00054 NOX::Direction::Broyden::BroydenMemoryUnit::BroydenMemoryUnit()
00055 {
00056 lambda = 0;
00057 snormsqr = 0;
00058 }
00059
00060 NOX::Direction::Broyden::BroydenMemoryUnit::~BroydenMemoryUnit()
00061 {
00062
00063 }
00064
00065 void NOX::Direction::Broyden::BroydenMemoryUnit::
00066 reset(const NOX::Abstract::Vector& d)
00067 {
00068
00069 if (Teuchos::is_null(sptr))
00070 sptr = d.clone(DeepCopy);
00071 else
00072 *sptr = d;
00073
00074
00075 snormsqr = sptr->norm();
00076 snormsqr = snormsqr * snormsqr;
00077
00078 lambda = 0.0;
00079 }
00080
00081 void NOX::Direction::Broyden::BroydenMemoryUnit::setStep(double step)
00082 {
00083 lambda = step;
00084 if (step != 1.0)
00085 {
00086 sptr->scale(step);
00087 snormsqr = step * step * snormsqr;
00088 }
00089 }
00090
00091 Teuchos::RCP<const NOX::Abstract::Vector> NOX::Direction::Broyden::
00092 BroydenMemoryUnit::sPtr() const
00093 {
00094 return sptr;
00095 }
00096
00097 double NOX::Direction::Broyden::BroydenMemoryUnit::step() const
00098 {
00099 return lambda;
00100 }
00101
00102 double NOX::Direction::Broyden::BroydenMemoryUnit::sNormSqr() const
00103 {
00104 return snormsqr;
00105 }
00106
00107
00108
00109 NOX::Direction::Broyden::BroydenMemory::BroydenMemory()
00110 {
00111 }
00112
00113 NOX::Direction::Broyden::BroydenMemory::~BroydenMemory()
00114 {
00115 }
00116
00117 void NOX::Direction::Broyden::BroydenMemory::reset(int m)
00118 {
00119 mMax = m;
00120
00121 if (memory.size() < (unsigned int) mMax)
00122 memory.resize(mMax);
00123
00124 if (index.capacity() < (unsigned int) mMax)
00125 index.reserve(mMax);
00126
00127 index.resize(0);
00128 }
00129
00130 void NOX::Direction::Broyden::BroydenMemory::reset()
00131 {
00132 index.resize(0);
00133 }
00134
00135 void NOX::Direction::Broyden::BroydenMemory::
00136 push(const NOX::Abstract::Vector& d)
00137 {
00138
00139 if (mMax == 0)
00140 return;
00141
00142
00143 int m = index.size();
00144
00145
00146
00147 if (m < mMax)
00148 {
00149 index.push_back(m);
00150 }
00151 else
00152 {
00153 int k = index[0];
00154 for (int i = 0; i < m - 1; i ++)
00155 index[i] = index[i+1];
00156 index[m-1] = k;
00157 }
00158
00159
00160 memory[index.back()].reset(d);
00161 }
00162
00163 bool NOX::Direction::Broyden::BroydenMemory::empty() const
00164 {
00165 return index.empty();
00166 }
00167
00168 int NOX::Direction::Broyden::BroydenMemory::size() const
00169 {
00170 return index.size();
00171 }
00172
00173 NOX::Direction::Broyden::BroydenMemoryUnit&
00174 NOX::Direction::Broyden::BroydenMemory::operator[](int i)
00175 {
00176 return memory[index[i]];
00177 }
00178
00179
00180
00181 NOX::Direction::Broyden::
00182 Broyden(const Teuchos::RCP<NOX::GlobalData>& gd,
00183 Teuchos::ParameterList& p) :
00184 lsParamsPtr(NULL),
00185 inexactNewtonUtils(gd, p)
00186 {
00187 reset(gd, p);
00188 }
00189
00190 NOX::Direction::Broyden::~Broyden()
00191 {
00192
00193 }
00194
00195 bool NOX::Direction::Broyden::
00196 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00197 Teuchos::ParameterList& params)
00198 {
00199 globalDataPtr = gd;
00200 utils = gd->getUtils();
00201
00202 Teuchos::ParameterList& p = params.sublist("Broyden");
00203
00204
00205 lsParamsPtr = &p.sublist("Linear Solver");
00206
00207
00208
00209
00210
00211
00212 inexactNewtonUtils.reset(gd, params);
00213
00214
00215 cntMax = p.get("Restart Frequency", 10);
00216
00217
00218 maxConvRate = p.get("Max Convergence Rate", 1.0);
00219
00220
00221 memorySizeMax = p.get("Memory", cntMax);
00222
00223
00224 memory.reset(memorySizeMax);
00225
00226 return true;
00227 }
00228
00229 bool NOX::Direction::Broyden::compute(NOX::Abstract::Vector& dir,
00230 NOX::Abstract::Group& soln,
00231 const NOX::Solver::Generic& solver)
00232 {
00233 throwError("compute", "This direction can only be used with a line search based solver.");
00234 return false;
00235 }
00236
00237 bool NOX::Direction::Broyden::compute(NOX::Abstract::Vector& dir,
00238 NOX::Abstract::Group& soln,
00239 const NOX::Solver::LineSearchBased& solver)
00240 {
00241
00242 NOX::Abstract::Group::ReturnType status;
00243
00244
00245 status = soln.computeF();
00246 if (status != NOX::Abstract::Group::Ok)
00247 throwError("compute", "Unable to compute F");
00248
00249
00250 if (doRestart(soln, solver))
00251 {
00252
00253 memory.reset();
00254
00255
00256 if (Teuchos::is_null(oldJacobianGrpPtr))
00257 oldJacobianGrpPtr = soln.clone(NOX::DeepCopy);
00258 else
00259
00260
00261
00262 *oldJacobianGrpPtr = soln;
00263
00264
00265 if (utils->isPrintType(NOX::Utils::Details))
00266 utils->out() << " Recomputing Jacobian" << endl;
00267
00268 status = oldJacobianGrpPtr->computeJacobian();
00269 if (status != NOX::Abstract::Group::Ok)
00270 throwError("compute", "Unable to compute Jacobian");
00271
00272
00273 cnt = 0;
00274 }
00275
00276
00277 if (!memory.empty())
00278 {
00279 double step = solver.getStepSize();
00280 memory[memory.size() - 1].setStep(step);
00281 }
00282
00283
00284
00285
00286 inexactNewtonUtils.computeForcingTerm(soln,
00287 solver.getPreviousSolutionGroup(),
00288 solver.getNumIterations(),
00289 solver);
00290
00291
00292 cnt ++;
00293 status = oldJacobianGrpPtr->applyJacobianInverse(*lsParamsPtr,
00294 soln.getF(),
00295 dir);
00296 if (status != NOX::Abstract::Group::Ok)
00297 throwError("compute", "Unable to apply Jacobian inverse");
00298 dir.scale(-1.0);
00299
00300
00301 if (!memory.empty())
00302 {
00303
00304 int m = memory.size();
00305
00306
00307 double step;
00308 Teuchos::RCP<const NOX::Abstract::Vector> sPtr;
00309
00310
00311
00312 double stepNext = memory[0].step();
00313 Teuchos::RCP<const NOX::Abstract::Vector> sPtrNext =
00314 memory[0].sPtr();
00315
00316
00317 double a, b, c, denom;
00318
00319 for (int i = 0; i < m-1; i ++)
00320 {
00321 step = stepNext;
00322 sPtr = sPtrNext;
00323 stepNext = memory[i+1].step();
00324 sPtrNext = memory[i+1].sPtr();
00325
00326 a = step / stepNext;
00327 b = step - 1;
00328 c = sPtr->innerProduct(dir) / memory[i].sNormSqr();
00329
00330 dir.update(a * c, *sPtrNext, b * c, *sPtr, 1.0);
00331 }
00332
00333 step = stepNext;
00334 sPtr = sPtrNext;
00335
00336 a = sPtr->innerProduct(dir);
00337 b = memory[m-1].sNormSqr();
00338 c = (step - 1) * a;
00339 denom = b - step * a;
00340
00341 dir.update(c / denom, *sPtr, b / denom);
00342 }
00343
00345 memory.push(dir);
00346
00347 return true;
00348 }
00349
00350 bool NOX::Direction::Broyden::doRestart(NOX::Abstract::Group& soln,
00351 const NOX::Solver::LineSearchBased& solver)
00352 {
00353
00354 if (solver.getNumIterations() == 0)
00355 return true;
00356
00357
00358 if (cnt >= cntMax)
00359 return true;
00360
00361
00362 if (solver.getStepSize() == 0.0)
00363 return true;
00364
00365
00366 convRate = soln.getNormF() / solver.getPreviousSolutionGroup().getNormF();
00367 if (convRate > maxConvRate)
00368 return true;
00369
00370 return false;
00371 }
00372
00373
00374 void NOX::Direction::Broyden::throwError(const string& functionName, const string& errorMsg)
00375 {
00376 if (utils->isPrintType(NOX::Utils::Error))
00377 utils->err() << "NOX::Direction::Broyden::" << functionName << " - " << errorMsg << endl;
00378 throw "NOX Error";
00379 }