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_Epetra_Group.H"
00043
00044 #include "Teuchos_ParameterList.hpp"
00045 #include "NOX_Utils.H"
00046 #include "NOX_Epetra_Interface_Required.H"
00047 #include "Epetra_Vector.h"
00048 #include "Epetra_Operator.h"
00049 #include "AztecOO_ConditionNumber.h"
00050
00051 using namespace NOX;
00052 using namespace NOX::Epetra;
00053 using namespace Teuchos;
00054
00055 Group::Group(Teuchos::ParameterList& printParams,
00056 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00057 const NOX::Epetra::Vector& x):
00058 utils(printParams),
00059 xVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(DeepCopy))),
00060 xVector(*xVectorPtr),
00061 RHSVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00062 RHSVector(*RHSVectorPtr),
00063 gradVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00064 gradVector(*gradVectorPtr),
00065 NewtonVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00066 NewtonVector(*NewtonVectorPtr),
00067 normNewtonSolveResidual(0),
00068 conditionNumber(0.0),
00069 sharedLinearSystem(*sharedLinearSystemPtr),
00070 userInterfacePtr(i)
00071 {
00072
00073 resetIsValid();
00074 }
00075
00076 Group::Group(Teuchos::ParameterList& printParams,
00077 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00078 const NOX::Epetra::Vector& x,
00079 const Teuchos::RCP<NOX::Epetra::LinearSystem>& linSys):
00080 utils(printParams),
00081 xVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(DeepCopy))),
00082 xVector(*xVectorPtr),
00083 RHSVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00084 RHSVector(*RHSVectorPtr),
00085 gradVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00086 gradVector(*gradVectorPtr),
00087 NewtonVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(x.clone(ShapeCopy))),
00088 NewtonVector(*NewtonVectorPtr),
00089 normNewtonSolveResidual(0),
00090 conditionNumber(0.0),
00091 sharedLinearSystemPtr(Teuchos::rcp(new SharedObject<NOX::Epetra::LinearSystem, NOX::Epetra::Group>(linSys))),
00092 sharedLinearSystem(*sharedLinearSystemPtr),
00093 userInterfacePtr(i)
00094 {
00095
00096 resetIsValid();
00097 }
00098
00099 Group::Group(const Group& source, CopyType type) :
00100 utils(source.utils),
00101 xVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(source.xVector.clone(type))),
00102 xVector(*xVectorPtr),
00103 RHSVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(source.RHSVector.clone(type))),
00104 RHSVector(*RHSVectorPtr),
00105 gradVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(source.gradVector.clone(type))),
00106 gradVector(*gradVectorPtr),
00107 NewtonVectorPtr(rcp_dynamic_cast<NOX::Epetra::Vector>(source.NewtonVector.clone(type))),
00108 NewtonVector(*NewtonVectorPtr),
00109 sharedLinearSystemPtr(source.sharedLinearSystemPtr),
00110 sharedLinearSystem(*sharedLinearSystemPtr),
00111 userInterfacePtr(source.userInterfacePtr)
00112 {
00113
00114 switch (type) {
00115
00116 case DeepCopy:
00117
00118 isValidRHS = source.isValidRHS;
00119 isValidJacobian = source.isValidJacobian;
00120 isValidGrad = source.isValidGrad;
00121 isValidNewton = source.isValidNewton;
00122 isValidNormNewtonSolveResidual = source.isValidNormNewtonSolveResidual;
00123 isValidConditionNumber = source.isValidConditionNumber;
00124 normNewtonSolveResidual = source.normNewtonSolveResidual;
00125 conditionNumber = source.conditionNumber;
00126 isValidPreconditioner = source.isValidPreconditioner;
00127 isValidSolverJacOp = source.isValidSolverJacOp;
00128
00129
00130 if (isValidJacobian)
00131 sharedLinearSystem.getObject(this);
00132
00133 break;
00134
00135 case ShapeCopy:
00136 resetIsValid();
00137 break;
00138
00139 default:
00140 cerr << "ERROR: Invalid ConstructorType for group copy constructor." << endl;
00141 throw "NOX Error";
00142 }
00143
00144 }
00145
00146 Group::~Group()
00147 {
00148 }
00149
00150 void Group::resetIsValid()
00151 {
00152 isValidRHS = false;
00153 isValidJacobian = false;
00154 isValidGrad = false;
00155 isValidNewton = false;
00156 isValidNormNewtonSolveResidual = false;
00157 isValidPreconditioner = false;
00158 isValidSolverJacOp = false;
00159 isValidConditionNumber = false;
00160 return;
00161 }
00162
00163 Teuchos::RCP<NOX::Abstract::Group> Group::clone(CopyType type) const
00164 {
00165 Teuchos::RCP<NOX::Abstract::Group> newgrp =
00166 Teuchos::rcp(new NOX::Epetra::Group(*this, type));
00167 return newgrp;
00168 }
00169
00170 Abstract::Group& Group::operator=(const NOX::Abstract::Group& source)
00171 {
00172 return operator=(dynamic_cast<const Group&> (source));
00173 }
00174
00175 Abstract::Group& Group::operator=(const Group& source)
00176 {
00177
00178 xVector = source.xVector;
00179
00180
00181
00182
00183
00184 isValidRHS = source.isValidRHS;
00185 isValidGrad = source.isValidGrad;
00186 isValidNewton = source.isValidNewton;
00187 isValidJacobian = source.isValidJacobian;
00188 isValidNormNewtonSolveResidual = source.isValidNormNewtonSolveResidual;
00189 isValidPreconditioner = source.isValidPreconditioner;
00190 isValidSolverJacOp = source.isValidSolverJacOp;
00191 isValidConditionNumber = source.isValidConditionNumber;
00192
00193
00194 if (isValidRHS) {
00195 RHSVector = source.RHSVector;
00196 }
00197
00198 if (isValidGrad)
00199 gradVector = source.gradVector;
00200
00201 if (isValidNewton)
00202 NewtonVector = source.NewtonVector;
00203
00204 if (isValidNormNewtonSolveResidual)
00205 normNewtonSolveResidual = source.normNewtonSolveResidual;
00206
00207
00208 if (isValidJacobian)
00209 sharedLinearSystem.getObject(this);
00210
00211
00212 if (isValidConditionNumber)
00213 conditionNumber = source.conditionNumber;
00214
00215 return *this;
00216 }
00217
00218 void Group::setX(const Abstract::Vector& y)
00219 {
00220 setX(dynamic_cast<const NOX::Epetra::Vector&> (y));
00221 return;
00222 }
00223
00224 void Group::setX(const NOX::Epetra::Vector& y)
00225 {
00226 if (isPreconditioner()) {
00227 sharedLinearSystem.getObject(this)->destroyPreconditioner();
00228 }
00229 resetIsValid();
00230 xVector = y;
00231 return;
00232 }
00233
00234 void Group::computeX(const NOX::Abstract::Group& grp,
00235 const NOX::Abstract::Vector& d,
00236 double step)
00237 {
00238
00239 const NOX::Epetra::Group& epetragrp = dynamic_cast<const Group&> (grp);
00240 const NOX::Epetra::Vector& epetrad =
00241 dynamic_cast<const NOX::Epetra::Vector&> (d);
00242 computeX(epetragrp, epetrad, step);
00243 return;
00244 }
00245
00246 void Group::computeX(const Group& grp,
00247 const NOX::Epetra::Vector& d,
00248 double step)
00249 {
00250 if (isPreconditioner())
00251 sharedLinearSystem.getObject(this)->destroyPreconditioner();
00252 resetIsValid();
00253 xVector.update(1.0, grp.xVector, step, d);
00254 return;
00255 }
00256
00257 Abstract::Group::ReturnType Group::computeF()
00258 {
00259 if (isF())
00260 return Abstract::Group::Ok;
00261
00262 bool status = false;
00263
00264 status = userInterfacePtr->computeF(xVector.getEpetraVector(), RHSVector.getEpetraVector(), NOX::Epetra::Interface::Required::Residual);
00265
00266 if (status == false) {
00267 cout << "ERROR: Epetra::Group::computeF() - fill failed!!!"
00268 << endl;
00269 throw "NOX Error: Fill Failed";
00270 }
00271
00272 isValidRHS = true;
00273
00274 return Abstract::Group::Ok;
00275 }
00276
00277 Abstract::Group::ReturnType Group::computeJacobian()
00278 {
00279
00280 if (isJacobian())
00281 return Abstract::Group::Ok;
00282
00283
00284 bool status = false;
00285
00286 status = sharedLinearSystem.getObject(this)->
00287 computeJacobian(xVector);
00288
00289 if (status == false) {
00290 cout << "ERROR: NOX::Epetra::Group::computeJacobian() - fill failed!!!"
00291 << endl;
00292 throw "NOX Error: Fill Failed";
00293 }
00294
00295
00296 isValidJacobian = true;
00297
00298 return Abstract::Group::Ok;
00299 }
00300
00301 Abstract::Group::ReturnType Group::computeGradient()
00302 {
00303 if (isGradient())
00304 return Abstract::Group::Ok;
00305
00306 if (!isF()) {
00307 cerr << "ERROR: NOX::Epetra::Group::computeGradient() - RHS is out of date wrt X!" << endl;
00308 throw "NOX Error";
00309 }
00310
00311 if (!isJacobian()) {
00312 cerr << "ERROR: NOX::Epetra::Group::computeGradient() - Jacobian is out of date wrt X!" << endl;
00313 throw "NOX Error";
00314 }
00315
00316
00317 sharedLinearSystem.getObject(this)->applyJacobianTranspose(RHSVector,
00318 gradVector);
00319
00320
00321 isValidGrad = true;
00322
00323
00324 return Abstract::Group::Ok;
00325 }
00326
00327 Abstract::Group::ReturnType Group::computeNewton(Teuchos::ParameterList& p)
00328 {
00329 if (isNewton())
00330 return Abstract::Group::Ok;
00331
00332 if (!isF()) {
00333 cerr << "ERROR: NOX::Epetra::Group::computeNewton() - invalid RHS" << endl;
00334 throw "NOX Error";
00335 }
00336
00337 if (!isJacobian()) {
00338 cerr << "ERROR: NOX::Epetra::Group::computeNewton() - invalid Jacobian" << endl;
00339 throw "NOX Error";
00340 }
00341
00342 Abstract::Group::ReturnType status;
00343
00344
00345 NewtonVector.init(0.0);
00346
00347
00348 status = applyJacobianInverse(p, RHSVector, NewtonVector);
00349
00350
00351 NewtonVector.scale(-1.0);
00352
00353
00354
00355 isValidNewton = true;
00356
00357
00358 computeNormNewtonSolveResidual();
00359
00360
00361 return status;
00362 }
00363
00364 Abstract::Group::ReturnType Group::applyJacobian(const Abstract::Vector& input, Abstract::Vector& result) const
00365 {
00366 const NOX::Epetra::Vector& epetrainput =
00367 dynamic_cast<const NOX::Epetra::Vector&> (input);
00368 NOX::Epetra::Vector& epetraresult =
00369 dynamic_cast<NOX::Epetra::Vector&> (result);
00370 return applyJacobian(epetrainput, epetraresult);
00371 }
00372
00373 Abstract::Group::ReturnType Group::applyJacobian(const NOX::Epetra::Vector& input, NOX::Epetra::Vector& result) const
00374 {
00375
00376 if (!isJacobian())
00377 return Abstract::Group::BadDependency;
00378
00379
00380 bool status = sharedLinearSystem.getObject()->applyJacobian(input, result);
00381
00382 return status == true ? Abstract::Group::Ok : Abstract::Group::Failed;
00383 }
00384
00385 Abstract::Group::ReturnType Group::applyJacobianInverse (Teuchos::ParameterList &p, const Abstract::Vector &input, Abstract::Vector &result) const
00386 {
00387 const NOX::Epetra::Vector& epetraInput = dynamic_cast<const NOX::Epetra::Vector&>(input);
00388 NOX::Epetra::Vector& epetraResult = dynamic_cast<NOX::Epetra::Vector&>(result);
00389 return applyJacobianInverse(p, epetraInput, epetraResult);
00390 }
00391
00392 Abstract::Group::ReturnType Group::applyJacobianInverse (Teuchos::ParameterList &p, const NOX::Epetra::Vector &input, NOX::Epetra::Vector &result) const
00393 {
00394 if (!isJacobian())
00395 return Abstract::Group::BadDependency;
00396
00397 if (!isValidSolverJacOp) {
00398 sharedLinearSystem.getObject(this)->setJacobianOperatorForSolve(sharedLinearSystem.getObject(this)->getJacobianOperator());
00399 isValidSolverJacOp = true;
00400 }
00401
00402
00403 NOX::Epetra::LinearSystem::PreconditionerReusePolicyType precPolicy =
00404 sharedLinearSystem.getObject(this)->getPreconditionerPolicy();
00405
00406 if (!isPreconditioner()) {
00407 if (precPolicy == NOX::Epetra::LinearSystem::PRPT_REBUILD) {
00408 sharedLinearSystem.getObject(this)->destroyPreconditioner();
00409 sharedLinearSystem.getObject(this)->
00410 createPreconditioner(xVector, p, false);
00411 isValidPreconditioner = true;
00412 }
00413 else if (precPolicy == NOX::Epetra::LinearSystem::PRPT_RECOMPUTE) {
00414 sharedLinearSystem.getObject(this)->recomputePreconditioner(xVector, p);
00415 isValidPreconditioner = true;
00416 }
00417 else if (precPolicy == NOX::Epetra::LinearSystem::PRPT_REUSE) {
00418
00419 }
00420 }
00421
00422 bool status = sharedLinearSystem.getObject(this)->applyJacobianInverse(p, input, result);
00423
00424 return status == true ? Abstract::Group::Ok : Abstract::Group::NotConverged;
00425 }
00426
00427 Abstract::Group::ReturnType Group::applyJacobianTranspose(const Abstract::Vector& input, Abstract::Vector& result) const
00428 {
00429 const NOX::Epetra::Vector& epetrainput = dynamic_cast<const NOX::Epetra::Vector&> (input);
00430 NOX::Epetra::Vector& epetraresult = dynamic_cast<NOX::Epetra::Vector&> (result);
00431 return applyJacobianTranspose(epetrainput, epetraresult);
00432 }
00433
00434 Abstract::Group::ReturnType Group::applyJacobianTranspose(const NOX::Epetra::Vector& input, NOX::Epetra::Vector& result) const
00435 {
00436
00437 if (!isJacobian())
00438 return Abstract::Group::BadDependency;
00439
00440 bool status = sharedLinearSystem.getObject()->applyJacobianTranspose(input, result);
00441
00442 return status == true ? Abstract::Group::Ok : Abstract::Group::Failed;
00443 }
00444
00445 Abstract::Group::ReturnType Group::applyRightPreconditioning(
00446 bool useTranspose,
00447 Teuchos::ParameterList& params,
00448 const Abstract::Vector& input,
00449 Abstract::Vector& result) const
00450 {
00451 const NOX::Epetra::Vector& epetraInput = dynamic_cast<const NOX::Epetra::Vector&>(input);
00452 NOX::Epetra::Vector& epetraResult = dynamic_cast<NOX::Epetra::Vector&>(result);
00453
00454 return applyRightPreconditioning(useTranspose, params, epetraInput, epetraResult);
00455 }
00456
00457 Abstract::Group::ReturnType Group::applyRightPreconditioning(
00458 bool useTranspose,
00459 Teuchos::ParameterList& linearSolverParams,
00460 const NOX::Epetra::Vector& input,
00461 NOX::Epetra::Vector& result) const
00462 {
00463
00464 bool success = false;
00465
00466 if (!isPreconditioner()) {
00467 sharedLinearSystem.getObject(this)->destroyPreconditioner();
00468 sharedLinearSystem.getObject(this)->
00469 createPreconditioner(xVector, linearSolverParams, false);
00470 isValidPreconditioner = true;
00471 }
00472
00473 success = sharedLinearSystem.getObject()->
00474 applyRightPreconditioning(useTranspose, linearSolverParams, input, result);
00475
00476 if (success == true)
00477 return Abstract::Group::Ok;
00478 else
00479 return Abstract::Group::Failed;
00480 }
00481
00482 bool Group::isF() const
00483 {
00484 return isValidRHS;
00485 }
00486
00487 bool Group::isJacobian() const
00488 {
00489 return ((sharedLinearSystem.isOwner(this)) && (isValidJacobian));
00490 }
00491
00492 bool Group::isGradient() const
00493 {
00494 return isValidGrad;
00495 }
00496
00497 bool Group::isNewton() const
00498 {
00499 return isValidNewton;
00500 }
00501
00502 bool Group::isNormNewtonSolveResidual() const
00503 {
00504 return isValidNormNewtonSolveResidual;
00505 }
00506
00507 bool Group::isPreconditioner() const
00508 {
00509 return ((sharedLinearSystem.isOwner(this)) && (isValidPreconditioner) &&
00510 (sharedLinearSystem.getObject(this)->isPreconditionerConstructed()));
00511 }
00512
00513 bool Group::isConditionNumber() const
00514 {
00515 return isValidConditionNumber;
00516 }
00517
00518 const Abstract::Vector& Group::getX() const
00519 {
00520 return xVector;
00521 }
00522
00523 const Abstract::Vector& Group::getF() const
00524 {
00525 if (!isF()) {
00526 cerr << "ERROR: NOX::Epetra::Group::getF() - invalid RHS" << endl;
00527 throw "NOX Error";
00528 }
00529
00530 return RHSVector;
00531 }
00532
00533 double Group::getNormF() const
00534 {
00535 if (!isF()) {
00536 cerr << "ERROR: NOX::Epetra::Group::getNormF() - invalid RHS" << endl;
00537 throw "NOX Error";
00538 }
00539
00540 return RHSVectorPtr->norm();
00541 }
00542
00543 const Abstract::Vector& Group::getGradient() const
00544 {
00545 if (!isGradient()) {
00546 cerr << "ERROR: NOX::Epetra::Group::getGradient() - invalid gradient" << endl;
00547 throw "NOX Error";
00548 }
00549
00550 return gradVector;
00551 }
00552
00553 const Abstract::Vector& Group::getNewton() const
00554 {
00555 if (!isNewton()) {
00556 cerr << "ERROR: NOX::Epetra::Group::getNewton() - invalid Newton vector" << endl;
00557 throw "NOX Error";
00558 }
00559
00560 return NewtonVector;
00561 }
00562
00563 Abstract::Group::ReturnType NOX::Epetra::Group::getNormLastLinearSolveResidual(double& residual) const
00564 {
00565
00566 if (isValidNormNewtonSolveResidual) {
00567 residual = normNewtonSolveResidual;
00568 return NOX::Abstract::Group::Ok;
00569 }
00570
00571
00572
00573 if (utils.isPrintType(Utils::Warning)) {
00574 cout << "ERROR: NOX::Epetra::Group::getNormLastLinearSolveResidual() - "
00575 << "Group has not performed a Newton solve corresponding to this "
00576 << "solution vector!" << endl;
00577 }
00578 return NOX::Abstract::Group::BadDependency;
00579 }
00580
00581 Teuchos::RCP<NOX::Epetra::Interface::Required> Group::
00582 getRequiredInterface()
00583 {
00584 return userInterfacePtr;
00585 }
00586
00587 Teuchos::RCP<const NOX::Epetra::LinearSystem> Group::
00588 getLinearSystem() const
00589 {
00590 return sharedLinearSystem.getObject();
00591 }
00592
00593 Teuchos::RCP<NOX::Epetra::LinearSystem> Group::getLinearSystem()
00594 {
00595 return sharedLinearSystem.getObject(this);
00596 }
00597
00598 bool Group::computeNormNewtonSolveResidual ()
00599 {
00600
00601 if (isValidNormNewtonSolveResidual)
00602 return true;
00603
00604
00605
00606 if (!isValidRHS) {
00607 cerr << "ERROR: NOX::Epetra::Group::computeNormNewtonSolveResidual() - invalid RHS"
00608 << endl;
00609 throw "NOX Error";
00610 }
00611 if (!isValidNewton) {
00612 cerr << "ERROR: NOX::Epetra::Group::computeNormNewtonSolveResidual() - invalid "
00613 << "Newton direction" << endl;
00614 throw "NOX Error";
00615 }
00616
00617
00618 if (Teuchos::is_null(tmpVectorPtr)) {
00619 tmpVectorPtr =
00620 Teuchos::rcp(new Epetra_Vector(RHSVector.getEpetraVector()));
00621 }
00622 NOX::Epetra::Vector tmpNoxVector(*tmpVectorPtr, ShapeCopy);
00623
00624 sharedLinearSystem.getObject()->applyJacobian(NewtonVector, tmpNoxVector);
00625 tmpNoxVector.update(1.0, RHSVector, 1.0);
00626 normNewtonSolveResidual = tmpNoxVector.norm();
00627
00628 isValidNormNewtonSolveResidual = true;
00629
00630 return true;
00631 }
00632
00633 Abstract::Group::ReturnType NOX::Epetra::Group::
00634 computeJacobianConditionNumber(int maxIters, double tolerance,
00635 int krylovSubspaceSize, bool printOutput)
00636 {
00637 if (!isConditionNumber()) {
00638 if (!isJacobian()) {
00639 cerr << "ERROR: NOX::Epetra::Group::computeJacobianConditionNumber()"
00640 << " - Jacobian is invalid wrt the solution." << endl;
00641 throw "NOX Error";
00642 }
00643
00644 if (Teuchos::is_null(azConditionNumberPtr))
00645 azConditionNumberPtr = Teuchos::rcp(new AztecOOConditionNumber);
00646
00647 azConditionNumberPtr->
00648 initialize(*(sharedLinearSystem.getObject()->getJacobianOperator()),
00649 AztecOOConditionNumber::GMRES_, krylovSubspaceSize,
00650 printOutput);
00651
00652 azConditionNumberPtr->computeConditionNumber(maxIters, tolerance);
00653
00654 conditionNumber = azConditionNumberPtr->getConditionNumber();
00655
00656 isValidConditionNumber = true;
00657 }
00658 return NOX::Abstract::Group::Ok;
00659 }
00660
00661 double NOX::Epetra::Group::getJacobianConditionNumber() const
00662 {
00663 if (!isConditionNumber()) {
00664 cerr << "ERROR: NOX::Epetra::Group::getJacobianConditionNumber()"
00665 << " - condition number has not yet been computed!" << endl;
00666 throw "NOX Error";
00667 }
00668
00669 return conditionNumber;
00670 }