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_Belos_Group.H"
00043 #include "NOX_Belos_MultiVector.H"
00044 #include "NOX_Belos_JacobianOperator.H"
00045 #include "NOX_Belos_PreconditionOperator.H"
00046 #include "NOX_Parameter_List.H"
00047
00048 #include "BelosConfigDefs.hpp"
00049 #include "BelosLinearProblemManager.hpp"
00050 #include "BelosOutputManager.hpp"
00051 #include "BelosStatusTestMaxIters.hpp"
00052 #include "BelosStatusTestMaxRestarts.hpp"
00053 #include "BelosStatusTestResNorm.hpp"
00054 #include "BelosStatusTestCombo.hpp"
00055 #include "BelosBlockGmres.hpp"
00056 #include "BelosBlockCG.hpp"
00057
00058 NOX::Belos::Group::Group(const Teuchos::RCP<NOX::Abstract::Group>& g,
00059 NOX::Parameter::List& printParams)
00060 : grpPtr(g),
00061 newtonVecPtr(g->getX().clone(NOX::ShapeCopy)),
00062 isValidNewton(false),
00063 utils(printParams),
00064 myPID(printParams.getParameter("MyPID", 0))
00065 {
00066 }
00067
00068 NOX::Belos::Group::Group(const NOX::Belos::Group& source,
00069 NOX::CopyType type)
00070 : grpPtr(source.grpPtr->clone(type)),
00071 newtonVecPtr(source.newtonVecPtr->clone(type)),
00072 isValidNewton(false),
00073 utils(source.utils),
00074 myPID(source.myPID)
00075 {
00076 if (type == NOX::DeepCopy)
00077 isValidNewton = source.isValidNewton;
00078 }
00079
00080
00081 NOX::Belos::Group::~Group()
00082 {
00083
00084 }
00085
00086 NOX::Belos::Group&
00087 NOX::Belos::Group::operator=(const NOX::Belos::Group& source)
00088 {
00089
00090
00091 if (this != &source) {
00092
00093
00094 *grpPtr = *source.grpPtr;
00095 *newtonVecPtr = *source.newtonVecPtr;
00096 isValidNewton = source.isValidNewton;
00097 utils = source.utils;
00098 myPID = source.myPID;
00099
00100 }
00101
00102 return *this;
00103 }
00104
00105 NOX::Abstract::Group&
00106 NOX::Belos::Group::operator=(const NOX::Abstract::Group& source)
00107 {
00108 return *this =
00109 dynamic_cast<const NOX::Belos::Group&>(source);
00110 }
00111
00112 Teuchos::RCP<NOX::Abstract::Group>
00113 NOX::Belos::Group::clone(NOX::CopyType type) const
00114 {
00115 Teuchos::RCP<NOX::Belos::Group> newGrp =
00116 Teuchos::rcp(new NOX::Belos::Group(*this, type));
00117 return newGrp;
00118 }
00119
00120 void
00121 NOX::Belos::Group::setX(const NOX::Abstract::Vector& y)
00122 {
00123 grpPtr->setX(y);
00124 resetIsValid();
00125 }
00126
00127 void
00128 NOX::Belos::Group::computeX(const NOX::Abstract::Group& g,
00129 const NOX::Abstract::Vector& d,
00130 double step)
00131 {
00132 const NOX::Belos::Group& belos_g = dynamic_cast<const NOX::Belos::Group&>(g);
00133 grpPtr->computeX(*(belos_g.grpPtr),d,step);
00134 resetIsValid();
00135 }
00136
00137 NOX::Abstract::Group::ReturnType
00138 NOX::Belos::Group::computeF()
00139 {
00140 return grpPtr->computeF();
00141 }
00142
00143 NOX::Abstract::Group::ReturnType
00144 NOX::Belos::Group::computeJacobian()
00145 {
00146 return grpPtr->computeJacobian();
00147 }
00148
00149 NOX::Abstract::Group::ReturnType
00150 NOX::Belos::Group::computeGradient()
00151 {
00152 return grpPtr->computeGradient();
00153 }
00154
00155 NOX::Abstract::Group::ReturnType
00156 NOX::Belos::Group::computeNewton(NOX::Parameter::List& params)
00157 {
00158 if (isValidNewton)
00159 return NOX::Abstract::Group::Ok;
00160
00161 if (!isF()) {
00162 cerr << "ERROR: NOX::Belos::Group::computeNewton() - invalid RHS" << endl;
00163 throw "NOX Error";
00164 }
00165
00166 if (!isJacobian()) {
00167 cerr << "ERROR: NOX::Belos::Group::computeNewton() - invalid Jacobian"
00168 << endl;
00169 throw "NOX Error";
00170 }
00171
00172 Abstract::Group::ReturnType status;
00173
00174
00175 newtonVecPtr->init(0.0);
00176
00177 status = applyJacobianInverse(params, getF(), *newtonVecPtr);
00178
00179 newtonVecPtr->scale(-1.0);
00180
00181
00182
00183 isValidNewton = true;
00184
00185 return status;
00186 }
00187
00188 NOX::Abstract::Group::ReturnType
00189 NOX::Belos::Group::applyJacobian(const NOX::Abstract::Vector& input,
00190 NOX::Abstract::Vector& result) const
00191 {
00192 return grpPtr->applyJacobian(input, result);
00193 }
00194
00195 NOX::Abstract::Group::ReturnType
00196 NOX::Belos::Group::applyJacobianTranspose(
00197 const NOX::Abstract::Vector& input,
00198 NOX::Abstract::Vector& result) const
00199 {
00200 return grpPtr->applyJacobianTranspose(input, result);
00201 }
00202
00203 NOX::Abstract::Group::ReturnType
00204 NOX::Belos::Group::applyJacobianInverse(NOX::Parameter::List& params,
00205 const NOX::Abstract::Vector& input,
00206 NOX::Abstract::Vector& result) const
00207 {
00208
00209 Teuchos::RCP<NOX::Abstract::MultiVector> inputs =
00210 input.createMultiVector(NULL, 0, NOX::DeepCopy);
00211 Teuchos::RCP<NOX::Abstract::MultiVector> results =
00212 result.createMultiVector(NULL, 0, NOX::DeepCopy);
00213
00214
00215 NOX::Abstract::Group::ReturnType res =
00216 applyJacobianInverseMultiVector(params, *inputs, *results);
00217
00218
00219 result = (*results)[0];
00220
00221 return res;
00222 }
00223
00224 NOX::Abstract::Group::ReturnType
00225 NOX::Belos::Group::applyRightPreconditioning(
00226 bool useTranspose,
00227 NOX::Parameter::List& params,
00228 const NOX::Abstract::Vector& input,
00229 NOX::Abstract::Vector& result) const
00230 {
00231 return grpPtr->applyRightPreconditioning(useTranspose, params, input,
00232 result);
00233 }
00234
00235 NOX::Abstract::Group::ReturnType
00236 NOX::Belos::Group::applyJacobianMultiVector(
00237 const NOX::Abstract::MultiVector& input,
00238 NOX::Abstract::MultiVector& result) const
00239 {
00240 return grpPtr->applyJacobianMultiVector(input, result);
00241 }
00242
00243 NOX::Abstract::Group::ReturnType
00244 NOX::Belos::Group::applyJacobianTransposeMultiVector(
00245 const NOX::Abstract::MultiVector& input,
00246 NOX::Abstract::MultiVector& result) const
00247 {
00248 return grpPtr->applyJacobianTransposeMultiVector(input, result);
00249 }
00250
00251 NOX::Abstract::Group::ReturnType
00252 NOX::Belos::Group::applyJacobianInverseMultiVector(
00253 NOX::Parameter::List& params,
00254 const NOX::Abstract::MultiVector& input,
00255 NOX::Abstract::MultiVector& result) const
00256 {
00257
00258 NOX::Abstract::MultiVector& nonConstInput =
00259 const_cast<NOX::Abstract::MultiVector&>(input);
00260
00261
00262 NOX::Belos::JacobianOperator belosJacOp(*grpPtr);
00263
00264
00265 NOX::Belos::PreconditionOperator belosPrecOp(*grpPtr, params);
00266
00267
00268 NOX::Belos::MultiVector belosInput(nonConstInput);
00269 NOX::Belos::MultiVector belosResult(result);
00270
00271
00272 ::Belos::LinearProblemManager<double> belosProblem(&belosJacOp,
00273 &belosResult,
00274 &belosInput);
00275 belosProblem.SetRightPrec(&belosPrecOp);
00276
00277
00278 int maxits = params.getParameter("Max Iterations", 400);
00279 double tol = params.getParameter("Tolerance", 1.0e-6);
00280 int length = params.getParameter("Size of Krylov Subspace", 300);
00281 int numrestarts = params.getParameter("Number of Restarts", 20);
00282 int maxblocksize = params.getParameter("Maximum block size", 10);
00283 string method = params.getParameter("Belos Solver", "GMRES");
00284 int verbLevel = params.getParameter("Verbosity Level", 1);
00285
00286
00287 ::Belos::StatusTestMaxIters<double> test1(maxits);
00288 ::Belos::StatusTestMaxRestarts<double> test2(numrestarts);
00289 ::Belos::StatusTestCombo<double> belosTest(
00290 ::Belos::StatusTestCombo<double>::OR,
00291 test1, test2 );
00292 ::Belos::StatusTestResNorm<double> test3( tol );
00293 test3.DefineScaleForm(::Belos::StatusTestResNorm<double>::NormOfPrecInitRes,
00294 ::Belos::TwoNorm );
00295 belosTest.AddStatusTest( test3 );
00296
00297
00298 int blocksize = input.numVectors();
00299 if (blocksize > maxblocksize)
00300 blocksize = maxblocksize;
00301 belosProblem.SetBlockSize(blocksize);
00302
00303
00304 ::Belos::OutputManager<double> belosOutputManager(myPID, verbLevel);
00305
00306 if (method == "GMRES") {
00307 ::Belos::BlockGmres<double> belosGMRES(belosProblem, belosTest,
00308 belosOutputManager, length);
00309 belosGMRES.Solve();
00310 }
00311 else if (method == "CG") {
00312 ::Belos::BlockCG<double> belosCG(belosProblem, belosTest,
00313 belosOutputManager);
00314 belosCG.Solve();
00315 }
00316 else {
00317 cout << "ERROR: NOX::Belos::Group::applyJacobianInverseMultiVector" << endl
00318 << "\"Belos Solver\" parameter \"" << method
00319 << "\" is invalid!" << endl;
00320 throw "NOX Error";
00321 }
00322
00323 ::Belos::StatusType status = belosTest.GetStatus();
00324 if (status == ::Belos::Unconverged)
00325 return NOX::Abstract::Group::NotConverged;
00326 else if (status == ::Belos::Converged)
00327 return NOX::Abstract::Group::Ok;
00328 else
00329 return NOX::Abstract::Group::Failed;
00330 }
00331
00332 NOX::Abstract::Group::ReturnType
00333 NOX::Belos::Group::applyRightPreconditioningMultiVector(
00334 bool useTranspose,
00335 NOX::Parameter::List& params,
00336 const NOX::Abstract::MultiVector& input,
00337 NOX::Abstract::MultiVector& result) const
00338 {
00339 return grpPtr->applyRightPreconditioningMultiVector(useTranspose, params,
00340 input, result);
00341 }
00342
00343 bool
00344 NOX::Belos::Group::isF() const
00345 {
00346 return grpPtr->isF();
00347 }
00348
00349 bool
00350 NOX::Belos::Group::isJacobian() const
00351 {
00352 return grpPtr->isJacobian();
00353 }
00354
00355 bool
00356 NOX::Belos::Group::isGradient() const
00357 {
00358 return grpPtr->isGradient();
00359 }
00360
00361 bool
00362 NOX::Belos::Group::isNewton() const
00363 {
00364 return isValidNewton;
00365 }
00366
00367 const NOX::Abstract::Vector&
00368 NOX::Belos::Group::getX() const
00369 {
00370 return grpPtr->getX();
00371 }
00372
00373 const NOX::Abstract::Vector&
00374 NOX::Belos::Group::getF() const
00375 {
00376 return grpPtr->getF();
00377 }
00378
00379 double
00380 NOX::Belos::Group::getNormF() const
00381 {
00382 return grpPtr->getNormF();
00383 }
00384
00385 const NOX::Abstract::Vector&
00386 NOX::Belos::Group::getGradient() const
00387 {
00388 return grpPtr->getGradient();
00389 }
00390
00391 const NOX::Abstract::Vector&
00392 NOX::Belos::Group::getNewton() const
00393 {
00394 return *newtonVecPtr;
00395 }
00396
00397 double
00398 NOX::Belos::Group::getNormNewtonSolveResidual() const
00399 {
00400 NOX::Abstract::Group::ReturnType status;
00401 Teuchos::RCP<NOX::Abstract::Vector> residual =
00402 getF().clone(NOX::DeepCopy);
00403
00404 status = applyJacobian(*newtonVecPtr, *residual);
00405 if (status != NOX::Abstract::Group::Ok) {
00406 cerr << "Error: NOX::Belos::Group::getNormNewtonSolveResidual() -- applyJacobian failed!" << endl;
00407 throw "NOX Error";
00408 }
00409
00410 residual->update(1.0, getF(), 1.0);
00411 double resid_norm = residual->norm();
00412
00413 return resid_norm;
00414 }
00415
00416 void
00417 NOX::Belos::Group::resetIsValid()
00418 {
00419 isValidNewton = false;
00420 }