00001
00002
00003
00004
00005
00006
00007
00008 #include "NOX_Common.H"
00009 #include "NOX_Playa_Group.hpp"
00010 #include "PlayaMPIComm.hpp"
00011 #include "PlayaOut.hpp"
00012 #include "PlayaTabs.hpp"
00013
00014 using Playa::Out;
00015 using Playa::Tabs;
00016 using Teuchos::RCP;
00017
00018
00019 NOX::NOXPlaya::Group::Group(const Playa::Vector<double>& initcond,
00020 const Playa::NonlinearOperator<double>& nonlinOp,
00021 const Playa::LinearSolver<double>& linsolver)
00022 :
00023 precision(3),
00024 xVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, DeepCopy))),
00025 fVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00026 newtonVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00027 gradientVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00028 solver(linsolver),
00029 jacobian(),
00030 nonlinearOp(nonlinOp),
00031 normF(0.0)
00032 {
00033 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00034 resetIsValid();
00035 }
00036
00037 NOX::NOXPlaya::Group::Group(const Playa::NonlinearOperator<double>& nonlinOp,
00038 const Playa::LinearSolver<double>& linsolver)
00039 :
00040 precision(3),
00041 xVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, DeepCopy))),
00042 fVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00043 newtonVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00044 gradientVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00045 solver(linsolver),
00046 jacobian(),
00047 nonlinearOp(nonlinOp),
00048 normF(0.0)
00049 {
00050 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00051 resetIsValid();
00052 }
00053
00054 NOX::NOXPlaya::Group::Group(const Playa::Vector<double>& initcond,
00055 const Playa::NonlinearOperator<double>& nonlinOp,
00056 const Playa::LinearSolver<double>& linsolver,
00057 int numdigits)
00058 :
00059 precision(numdigits),
00060 xVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, DeepCopy))),
00061 fVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00062 newtonVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00063 gradientVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00064 solver(linsolver),
00065 jacobian(),
00066 nonlinearOp(nonlinOp),
00067 normF(0.0)
00068 {
00069 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00070 resetIsValid();
00071 }
00072
00073 NOX::NOXPlaya::Group::Group(const Playa::NonlinearOperator<double>& nonlinOp,
00074 const Playa::LinearSolver<double>& linsolver,
00075 int numdigits)
00076 :
00077 precision(numdigits),
00078 xVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, DeepCopy))),
00079 fVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00080 newtonVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00081 gradientVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00082 solver(linsolver),
00083 jacobian(),
00084 nonlinearOp(nonlinOp),
00085 normF(0.0)
00086 {
00087 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00088 resetIsValid();
00089 }
00090
00091
00092 NOX::NOXPlaya::Group::Group(const NOX::NOXPlaya::Group& source, NOX::CopyType type) :
00093 precision(source.precision),
00094 xVector(rcp(new NOX::NOXPlaya::Vector(*(source.xVector), precision, type))),
00095 fVector(rcp(new NOX::NOXPlaya::Vector(*(source.fVector), precision, type))),
00096 newtonVector(rcp(new NOX::NOXPlaya::Vector(*(source.newtonVector), precision, type))),
00097 gradientVector(rcp(new NOX::NOXPlaya::Vector(*(source.gradientVector), precision, type))),
00098 solver(source.solver),
00099 jacobian(source.jacobian),
00100 nonlinearOp(source.nonlinearOp),
00101 isValidF(false),
00102 isValidJacobian(false),
00103 isValidGradient(false),
00104 isValidNewton(false),
00105 normF(0.0)
00106 {
00107 switch (type)
00108 {
00109
00110 case NOX::DeepCopy:
00111
00112 isValidF = source.isValidF;
00113 isValidGradient = source.isValidGradient;
00114 isValidNewton = source.isValidNewton;
00115 isValidJacobian = source.isValidJacobian;
00116 normF = source.normF;
00117 break;
00118
00119 case NOX::ShapeCopy:
00120 resetIsValid();
00121 normF = 0.0;
00122 break;
00123
00124 default:
00125 Out::os() << "NOX:Playa::Group - invalid CopyType for copy constructor." << std::endl;
00126 throw "NOX Playa Error";
00127 }
00128
00129 }
00130
00131 NOX::NOXPlaya::Group::~Group()
00132 {
00133 }
00134
00135 void NOX::NOXPlaya::Group::resetIsValid()
00136 {
00137 isValidF = false;
00138 isValidJacobian = false;
00139 isValidGradient = false;
00140 isValidNewton = false;
00141 }
00142
00143
00144 RCP<NOX::Abstract::Group> NOX::NOXPlaya::Group::clone(NOX::CopyType type) const
00145 {
00146 return rcp(new NOX::NOXPlaya::Group(*this, type));
00147 }
00148
00149 NOX::Abstract::Group& NOX::NOXPlaya::Group::operator=(const NOX::Abstract::Group& source)
00150 {
00151 return operator=(dynamic_cast<const NOX::NOXPlaya::Group&> (source));
00152 }
00153
00154 NOX::Abstract::Group& NOX::NOXPlaya::Group::operator=(const NOX::NOXPlaya::Group& source)
00155 {
00156 if (this != &source)
00157 {
00158
00159
00160 *xVector = *(source.xVector);
00161 nonlinearOp = source.nonlinearOp;
00162 solver = source.solver;
00163 jacobian = source.jacobian;
00164 precision = source.precision;
00165
00166
00167 isValidF = source.isValidF;
00168 isValidGradient = source.isValidGradient;
00169 isValidNewton = source.isValidNewton;
00170 isValidJacobian = source.isValidJacobian;
00171
00172
00173 if (isValidF)
00174 {
00175 *fVector = *(source.fVector);
00176 normF = source.normF;
00177 }
00178
00179 if (isValidGradient)
00180 *gradientVector = *(source.gradientVector);
00181
00182 if (isValidNewton)
00183 *newtonVector = *(source.newtonVector);
00184
00185 if (isValidJacobian)
00186 jacobian = source.jacobian;
00187 }
00188 return *this;
00189 }
00190
00191 void NOX::NOXPlaya::Group::setX(const NOX::Abstract::Vector& y)
00192 {
00193 setX(dynamic_cast<const NOX::NOXPlaya::Vector&> (y));
00194 }
00195
00196 void NOX::NOXPlaya::Group::setX(const NOX::NOXPlaya::Vector& y)
00197 {
00198 resetIsValid();
00199 nonlinearOp.setEvalPt(y.getPlayaVector());
00200 *xVector = y;
00201 }
00202
00203 void NOX::NOXPlaya::Group::computeX(const NOX::Abstract::Group& grp,
00204 const NOX::Abstract::Vector& d,
00205 double step)
00206 {
00207
00208 const Group& tsfgrp = dynamic_cast<const NOX::NOXPlaya::Group&> (grp);
00209 const Vector& tsfd = dynamic_cast<const NOX::NOXPlaya::Vector&> (d);
00210 computeX(tsfgrp, tsfd, step);
00211 }
00212
00213 void NOX::NOXPlaya::Group::computeX(const Group& grp, const Vector& d, double step)
00214 {
00215 Tabs tab;
00216 resetIsValid();
00217 xVector->update(1.0, *(grp.xVector), step, d);
00218 }
00219
00220 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeF()
00221 {
00222
00223 if (nonlinearOp.verb() > 2)
00224 {
00225 Out::os() << "calling computeF()" << std::endl;
00226 }
00227
00228 if (isValidF)
00229 {
00230 if (nonlinearOp.verb() > 2)
00231 {
00232 Out::os() << "reusing F" << std::endl;
00233 }
00234 return NOX::Abstract::Group::Ok;
00235 }
00236 else
00237 {
00238 if (nonlinearOp.verb() > 2)
00239 {
00240 Out::os() << "computing new F" << std::endl;
00241 }
00242 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00243 *fVector = nonlinearOp.getFunctionValue();
00244 isValidF = true;
00245 normF = fVector->norm();
00246 }
00247
00248 return (isValidF) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00249 }
00250
00251 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeJacobian()
00252 {
00253
00254 if (nonlinearOp.verb() > 2)
00255 {
00256 Out::os() << "calling computeJ()" << std::endl;
00257 }
00258
00259
00260 if (isValidJacobian)
00261 {
00262 return NOX::Abstract::Group::Ok;
00263 }
00264 else
00265 {
00266 nonlinearOp.setEvalPt(xVector->getPlayaVector());
00267 jacobian = nonlinearOp.getJacobian();
00268
00269 isValidJacobian = true;
00270 }
00271 return (isValidJacobian) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00272 }
00273
00274 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeGradient()
00275 {
00276 if (nonlinearOp.verb() > 2)
00277 {
00278 Out::os() << "calling computeGrad()" << std::endl;
00279 }
00280 if (isValidGradient)
00281 {
00282 return NOX::Abstract::Group::Ok;
00283 }
00284 else
00285 {
00286 if (!isF())
00287 {
00288 Out::os() << "ERROR: NOX::NOXPlaya::Group::computeGrad() - F is out of date wrt X!" << std::endl;
00289 return NOX::Abstract::Group::BadDependency;
00290 }
00291
00292 if (!isJacobian())
00293 {
00294 Out::os() << "ERROR: NOX::NOXPlaya::Group::computeGrad() - Jacobian is out of date wrt X!" << std::endl;
00295 return NOX::Abstract::Group::BadDependency;
00296 }
00297
00298
00299
00300 NOX::Abstract::Group::ReturnType status
00301 = applyJacobianTranspose(*fVector,*gradientVector);
00302 isValidGradient = (status == NOX::Abstract::Group::Ok);
00303
00304
00305 return status;
00306 }
00307 }
00308
00309 NOX::Abstract::Group::ReturnType
00310 NOX::NOXPlaya::Group::computeNewton(Teuchos::ParameterList& p)
00311 {
00312 if (isNewton())
00313 {
00314 return NOX::Abstract::Group::Ok;
00315 }
00316 else
00317 {
00318 if (!isF())
00319 {
00320 Out::os() << "ERROR: NOX::Example::Group::computeNewton() - invalid F"
00321 << std::endl;
00322 throw "NOX Error";
00323 }
00324
00325 if (!isJacobian())
00326 {
00327 Out::os() << "ERROR: NOX::Example::Group::computeNewton() - invalid Jacobian" << std::endl;
00328 throw "NOX Error";
00329 }
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 NOX::Abstract::Group::ReturnType status
00340 = applyJacobianInverse(p, *fVector, *newtonVector);
00341 isValidNewton = (status == NOX::Abstract::Group::Ok);
00342
00343
00344
00345 newtonVector->scale(-1.0);
00346
00347 if (nonlinearOp.verb() > 0)
00348 {
00349 Out::os() << "newton step" << std::endl;
00350 newtonVector->getPlayaVector().print(Out::os());
00351 }
00352
00353
00354 return status;
00355 }
00356 }
00357
00358
00359 NOX::Abstract::Group::ReturnType
00360 NOX::NOXPlaya::Group::applyJacobian(const Abstract::Vector& input,
00361 NOX::Abstract::Vector& result) const
00362 {
00363 const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00364 NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result);
00365 return applyJacobian(tsfinput, tsfresult);
00366 }
00367
00368 NOX::Abstract::Group::ReturnType
00369 NOX::NOXPlaya::Group::applyJacobian(const NOX::NOXPlaya::Vector& input,
00370 NOX::NOXPlaya::Vector& result) const
00371 {
00372
00373 if (!isJacobian())
00374 {
00375 return NOX::Abstract::Group::BadDependency;
00376 }
00377 else
00378 {
00379
00380 jacobian.apply(input.getPlayaVector(),result.getPlayaVector());
00381 return NOX::Abstract::Group::Ok;
00382 }
00383 }
00384
00385 NOX::Abstract::Group::ReturnType
00386 NOX::NOXPlaya::Group::applyJacobianTranspose(const NOX::Abstract::Vector& input,
00387 NOX::Abstract::Vector& result) const
00388 {
00389 const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00390 NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result);
00391 return applyJacobianTranspose(tsfinput, tsfresult);
00392 }
00393
00394 NOX::Abstract::Group::ReturnType
00395 NOX::NOXPlaya::Group::applyJacobianTranspose(const NOX::NOXPlaya::Vector& input, NOX::NOXPlaya::Vector& result) const
00396 {
00397
00398 if (!isJacobian())
00399 {
00400 return NOX::Abstract::Group::BadDependency;
00401 }
00402 else
00403 {
00404
00405 jacobian.applyTranspose(input.getPlayaVector(), result.getPlayaVector());
00406
00407 return NOX::Abstract::Group::Ok;
00408 }
00409 }
00410
00411 NOX::Abstract::Group::ReturnType
00412 NOX::NOXPlaya::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00413 const Abstract::Vector& input,
00414 NOX::Abstract::Vector& result) const
00415 {
00416 const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00417 NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result);
00418 return applyJacobianInverse(p, tsfinput, tsfresult);
00419 }
00420
00421 NOX::Abstract::Group::ReturnType
00422 NOX::NOXPlaya::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00423 const NOX::NOXPlaya::Vector& input,
00424 NOX::NOXPlaya::Vector& result) const
00425 {
00426
00427 if (!isJacobian())
00428 {
00429 Out::os() << "ERROR: NOX::NOXPlaya::Group::applyJacobianInverse() - invalid Jacobian" << std::endl;
00430 throw "NOX Error";
00431
00432 }
00433
00434
00435
00436
00437
00438
00439
00440 if (nonlinearOp.verb() > 4)
00441 {
00442 Out::os() << "---------------- applying J^-1 ------------------" << std::endl;
00443 Out::os() << "J=" << std::endl;
00444 jacobian.print(Out::os());
00445 Out::os() << "F=" << std::endl;
00446 input.getPlayaVector().print(Out::os());
00447 }
00448
00449 Playa::SolverState<double> status
00450 = solver.solve(jacobian, input.getPlayaVector(),
00451 result.getPlayaVector());
00452
00453 if (status.finalState() != Playa::SolveConverged)
00454 {
00455 return NOX::Abstract::Group::Failed;
00456 }
00457 else
00458 {
00459 if (nonlinearOp.verb() > 2)
00460 {
00461 Out::os() << "soln=" << std::endl;
00462 result.getPlayaVector().print(Out::os());
00463 }
00464 return NOX::Abstract::Group::Ok;
00465 }
00466
00467 }
00468
00469 bool NOX::NOXPlaya::Group::isF() const
00470 {
00471 return isValidF;
00472 }
00473
00474 bool NOX::NOXPlaya::Group::isJacobian() const
00475 {
00476 return isValidJacobian;
00477 }
00478
00479 bool NOX::NOXPlaya::Group::isGradient() const
00480 {
00481 return isValidGradient;
00482 }
00483
00484 bool NOX::NOXPlaya::Group::isNewton() const
00485 {
00486 return isValidNewton;
00487 }
00488
00489 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getX() const
00490 {
00491 return *xVector;
00492 }
00493
00494 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getF() const
00495 {
00496 if (nonlinearOp.verb() > 2)
00497 {
00498 Out::os() << "calling getF()" << std::endl;
00499 }
00500 TEUCHOS_TEST_FOR_EXCEPTION(!isF(), runtime_error,
00501 "calling getF() with invalid function value");
00502 return *fVector;
00503 }
00504
00505 double NOX::NOXPlaya::Group::getNormF() const
00506 {
00507 if (nonlinearOp.verb() > 2)
00508 {
00509 Out::os() << "normF = " << normF << std::endl;
00510 }
00511 TEUCHOS_TEST_FOR_EXCEPTION(!isF(), runtime_error,
00512 "calling normF() with invalid function value");
00513 return normF;
00514 }
00515
00516 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getGradient() const
00517 {
00518 return *gradientVector;
00519 }
00520
00521 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getNewton() const
00522 {
00523 return *newtonVector;
00524 }
00525
00526 void NOX::NOXPlaya::Group::print() const
00527 {
00528 cout << "x = " << *xVector << "\n";
00529
00530 if (isValidF) {
00531 cout << "F(x) = " << *fVector << "\n";
00532 cout << "|| F(x) || = " << normF << "\n";
00533 }
00534 else
00535 cout << "F(x) has not been computed" << "\n";
00536
00537 cout << std::endl;
00538 }
00539
00540
00541