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
00043 #include "NOX_Epetra_LinearSystem_AztecOO.H"
00044
00045
00046 #include "NOX_Epetra_Interface_Required.H"
00047 #include "NOX_Epetra_Interface_Jacobian.H"
00048 #include "NOX_Epetra_Interface_Preconditioner.H"
00049 #include "NOX_Epetra_MatrixFree.H"
00050 #include "NOX_Epetra_FiniteDifference.H"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "NOX_Epetra_Scaling.H"
00053 #include "NOX_Utils.H"
00054
00055
00056 #include "Epetra_Map.h"
00057 #include "Epetra_Vector.h"
00058 #include "Epetra_Operator.h"
00059 #include "Epetra_RowMatrix.h"
00060 #include "Epetra_VbrMatrix.h"
00061 #include "Epetra_CrsMatrix.h"
00062 #include "Epetra_LinearProblem.h"
00063 #include "AztecOO.h"
00064 #include "AztecOO_Operator.h"
00065 #include "AztecOO_StatusTest.h"
00066 #include "AztecOO_StatusTestCombo.h"
00067 #include "AztecOO_StatusTestMaxIters.h"
00068 #include "AztecOO_StatusTestResNorm.h"
00069 #include "Ifpack.h"
00070 #include "Ifpack_IlukGraph.h"
00071 #include "Ifpack_CrsRiluk.h"
00072 #include "Teuchos_ParameterList.hpp"
00073
00074
00075 #ifdef HAVE_NOX_DEBUG
00076 #ifdef HAVE_NOX_EPETRAEXT
00077 #include "EpetraExt_BlockMapOut.h"
00078 #include "EpetraExt_MultiVectorOut.h"
00079 #include "EpetraExt_RowMatrixOut.h"
00080 #endif
00081 #endif
00082
00083 #ifdef HAVE_NOX_ML_EPETRA
00084 #include "Teuchos_ParameterList.hpp"
00085 #endif
00086
00087 #include <typeinfo>
00088
00089
00090 NOX::Epetra::LinearSystemAztecOO::
00091 LinearSystemAztecOO(
00092 Teuchos::ParameterList& printParams,
00093 Teuchos::ParameterList& linearSolverParams,
00094 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
00095 const NOX::Epetra::Vector& cloneVector,
00096 const Teuchos::RCP<NOX::Epetra::Scaling> s):
00097 utils(printParams),
00098 jacType(EpetraOperator),
00099 precType(EpetraOperator),
00100 precMatrixSource(UseJacobian),
00101 scaling(s),
00102 conditionNumberEstimate(0.0),
00103 isPrecConstructed(false),
00104 precQueryCounter(0),
00105 maxAgeOfPrec(1),
00106 timer(cloneVector.getEpetraVector().Comm()),
00107 timeCreatePreconditioner(0.0),
00108 timeApplyJacbianInverse(0.0)
00109 {
00110
00111 aztecSolverPtr = Teuchos::rcp(new AztecOO());
00112 tmpVectorPtr = Teuchos::rcp(new NOX::Epetra::Vector(cloneVector));
00113
00114
00115 createJacobianOperator(printParams, linearSolverParams, iReq, cloneVector);
00116 createPrecOperator(printParams, linearSolverParams, iReq, cloneVector);
00117
00118 reset(linearSolverParams);
00119 }
00120
00121
00122 NOX::Epetra::LinearSystemAztecOO::
00123 LinearSystemAztecOO(
00124 Teuchos::ParameterList& printParams,
00125 Teuchos::ParameterList& linearSolverParams,
00126 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
00127 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
00128 const Teuchos::RCP<Epetra_Operator>& jacobian,
00129 const NOX::Epetra::Vector& cloneVector,
00130 const Teuchos::RCP<NOX::Epetra::Scaling> s):
00131 utils(printParams),
00132 jacInterfacePtr(iJac),
00133 jacType(EpetraOperator),
00134 jacPtr(jacobian),
00135 precType(EpetraOperator),
00136 precMatrixSource(UseJacobian),
00137 scaling(s),
00138 conditionNumberEstimate(0.0),
00139 isPrecConstructed(false),
00140 precQueryCounter(0),
00141 maxAgeOfPrec(1),
00142 timer(cloneVector.getEpetraVector().Comm()),
00143 timeCreatePreconditioner(0.0),
00144 timeApplyJacbianInverse(0.0)
00145 {
00146
00147 aztecSolverPtr = Teuchos::rcp(new AztecOO());
00148 tmpVectorPtr = Teuchos::rcp(new NOX::Epetra::Vector(cloneVector));
00149
00150
00151 jacType = getOperatorType(*jacPtr);
00152
00153 createPrecOperator(printParams, linearSolverParams, iReq, cloneVector);
00154
00155 reset(linearSolverParams);
00156 }
00157
00158
00159 NOX::Epetra::LinearSystemAztecOO::
00160 LinearSystemAztecOO(
00161 Teuchos::ParameterList& printParams,
00162 Teuchos::ParameterList& linearSolverParams,
00163 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
00164 const Teuchos::RCP<NOX::Epetra::Interface::Preconditioner>& iPrec,
00165 const Teuchos::RCP<Epetra_Operator>& preconditioner,
00166 const NOX::Epetra::Vector& cloneVector,
00167 const Teuchos::RCP<NOX::Epetra::Scaling> s):
00168 utils(printParams),
00169 jacType(EpetraOperator),
00170 precInterfacePtr(iPrec),
00171 precType(EpetraOperator),
00172 precPtr(preconditioner),
00173 precMatrixSource(SeparateMatrix),
00174 scaling(s),
00175 conditionNumberEstimate(0.0),
00176 isPrecConstructed(false),
00177 precQueryCounter(0),
00178 maxAgeOfPrec(1),
00179 timer(cloneVector.getEpetraVector().Comm()),
00180 timeCreatePreconditioner(0.0),
00181 timeApplyJacbianInverse(0.0)
00182 {
00183
00184 aztecSolverPtr = Teuchos::rcp(new AztecOO());
00185 tmpVectorPtr = Teuchos::rcp(new NOX::Epetra::Vector(cloneVector));
00186
00187
00188 createJacobianOperator(printParams, linearSolverParams, iReq, cloneVector);
00189
00190 precType = getOperatorType(*precPtr);
00191
00192 reset(linearSolverParams);
00193 }
00194
00195
00196 NOX::Epetra::LinearSystemAztecOO::
00197 LinearSystemAztecOO(
00198 Teuchos::ParameterList& printParams,
00199 Teuchos::ParameterList& linearSolverParams,
00200 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
00201 const Teuchos::RCP<Epetra_Operator>& jacobian,
00202 const Teuchos::RCP<NOX::Epetra::Interface::Preconditioner>& iPrec,
00203 const Teuchos::RCP<Epetra_Operator>& preconditioner,
00204 const NOX::Epetra::Vector& cloneVector,
00205 const Teuchos::RCP<NOX::Epetra::Scaling> s):
00206 utils(printParams),
00207 jacInterfacePtr(iJac),
00208 jacType(EpetraOperator),
00209 jacPtr(jacobian),
00210 precInterfacePtr(iPrec),
00211 precType(EpetraOperator),
00212 precPtr(preconditioner),
00213 precMatrixSource(SeparateMatrix),
00214 scaling(s),
00215 conditionNumberEstimate(0.0),
00216 isPrecConstructed(false),
00217 precQueryCounter(0),
00218 maxAgeOfPrec(1),
00219 timer(cloneVector.getEpetraVector().Comm()),
00220 timeCreatePreconditioner(0.0),
00221 timeApplyJacbianInverse(0.0)
00222 {
00223
00224 aztecSolverPtr = Teuchos::rcp(new AztecOO());
00225 tmpVectorPtr = Teuchos::rcp(new NOX::Epetra::Vector(cloneVector));
00226
00227
00228 jacType = getOperatorType(*jacPtr);
00229 precType = getOperatorType(*precPtr);
00230
00231 reset(linearSolverParams);
00232 }
00233
00234
00235 NOX::Epetra::LinearSystemAztecOO::~LinearSystemAztecOO()
00236 {
00237 destroyPreconditioner();
00238 }
00239
00240
00241 void NOX::Epetra::LinearSystemAztecOO::
00242 reset(Teuchos::ParameterList& linearSolverParams)
00243 {
00244
00245
00246 destroyPreconditioner();
00247
00248
00249 string prec = linearSolverParams.get("Preconditioner", "None");
00250 if (prec == "AztecOO")
00251 precAlgorithm = AztecOO_;
00252 else if (prec == "Ifpack")
00253 precAlgorithm = Ifpack_;
00254 else if (prec == "New Ifpack")
00255 precAlgorithm = NewIfpack_;
00256 #ifdef HAVE_NOX_ML_EPETRA
00257 else if (prec == "ML")
00258 precAlgorithm = ML_;
00259 #endif
00260 else if (prec == "User Defined")
00261 precAlgorithm = UserDefined_;
00262 else if (prec == "None")
00263 precAlgorithm = None_;
00264 else {
00265 string errorMessage = "Option for \"Preconditioner\" is invalid!";
00266 throwError("reset()", errorMessage);
00267 }
00268
00269
00270
00271 checkPreconditionerValidity();
00272
00273 zeroInitialGuess =
00274 linearSolverParams.get("Zero Initial Guess", false);
00275
00276 manualScaling =
00277 linearSolverParams.get("Compute Scaling Manually", true);
00278
00279
00280
00281 outputSolveDetails =
00282 linearSolverParams.get("Output Solver Details", true);
00283
00284 throwErrorOnPrecFailure =
00285 linearSolverParams.get("Throw Error on Prec Failure", true);
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310 setAztecOptions(linearSolverParams, *aztecSolverPtr);
00311
00312
00313 std::string preReusePolicyName =
00314 linearSolverParams.get("Preconditioner Reuse Policy", "Rebuild");
00315 if (preReusePolicyName == "Rebuild")
00316 precReusePolicy = PRPT_REBUILD;
00317 else if (preReusePolicyName == "Recompute")
00318 precReusePolicy = PRPT_RECOMPUTE;
00319 else if (preReusePolicyName == "Reuse")
00320 precReusePolicy = PRPT_REUSE;
00321 else {
00322 string errorMessage = "Option for \"Preconditioner Reuse Policy\" is invalid! \nPossible options are \"Reuse\", \"Rebuild\", and \"Recompute\".";
00323 throwError("reset()", errorMessage);
00324 }
00325 maxAgeOfPrec = linearSolverParams.get("Max Age Of Prec", 1);
00326 precQueryCounter = 0;
00327
00328 #ifdef HAVE_NOX_DEBUG
00329 #ifdef HAVE_NOX_EPETRAEXT
00330 linearSolveCount = 0;
00331 #endif
00332 #endif
00333
00334 }
00335
00336
00337 void NOX::Epetra::LinearSystemAztecOO::
00338 setAztecOptions(Teuchos::ParameterList& p, AztecOO& aztec) const
00339 {
00340
00341 aztec.SetOutputStream(utils.pout(NOX::Utils::LinearSolverDetails));
00342 aztec.SetErrorStream(utils.perr());
00343
00344
00345 string linearSolver = p.get("Aztec Solver", "GMRES");
00346 if (linearSolver == "CG")
00347 aztec.SetAztecOption(AZ_solver, AZ_cg);
00348 else if (linearSolver == "GMRES")
00349 aztec.SetAztecOption(AZ_solver, AZ_gmres);
00350 else if (linearSolver == "CGS")
00351 aztec.SetAztecOption(AZ_solver, AZ_cgs);
00352 else if (linearSolver == "TFQMR")
00353 aztec.SetAztecOption(AZ_solver, AZ_tfqmr);
00354 else if (linearSolver == "BiCGStab")
00355 aztec.SetAztecOption(AZ_solver, AZ_bicgstab);
00356 else if (linearSolver == "LU")
00357 aztec.SetAztecOption(AZ_solver, AZ_lu);
00358 else {
00359 utils.out() << "ERROR: NOX::Epetra::Group::setAztecOptions" << endl
00360 << "\"Aztec Solver\" parameter \"" << linearSolver
00361 << "\" is invalid!" << endl;
00362 throw "NOX Error";
00363 }
00364
00365
00366 if (precAlgorithm == AztecOO_) {
00367
00368 string aztecPreconditioner = p.get("Aztec Preconditioner", "ilu");
00369
00370 if (aztecPreconditioner == "ilu") {
00371 aztec.SetAztecOption(AZ_precond, AZ_dom_decomp);
00372 aztec.SetAztecOption(AZ_overlap, p.get("Overlap", 0));
00373 aztec.SetAztecOption(AZ_subdomain_solve, AZ_ilu);
00374 aztec.SetAztecOption(AZ_graph_fill, p.get("Graph Fill", 0));
00375 }
00376 else if (aztecPreconditioner == "ilut") {
00377 aztec.SetAztecOption(AZ_precond, AZ_dom_decomp);
00378 aztec.SetAztecOption(AZ_overlap, p.get("Overlap", 0));
00379 aztec.SetAztecOption(AZ_subdomain_solve, AZ_ilut);
00380 aztec.SetAztecParam(AZ_drop, p.get("Drop Tolerance", 0.0));
00381 aztec.SetAztecParam(AZ_ilut_fill, p.get("Fill Factor", 1.0));
00382 }
00383 else if (aztecPreconditioner == "Jacobi") {
00384 aztec.SetAztecOption(AZ_precond, AZ_Jacobi);
00385 aztec.SetAztecOption(AZ_poly_ord, p.get("Steps", 3));
00386 }
00387 else if (aztecPreconditioner == "Symmetric Gauss-Seidel") {
00388 aztec.SetAztecOption(AZ_precond, AZ_sym_GS);
00389 aztec.SetAztecOption(AZ_poly_ord, p.get("Steps", 3));
00390 }
00391 else if (aztecPreconditioner == "Polynomial") {
00392 aztec.SetAztecOption(AZ_precond, AZ_Neumann);
00393 aztec.SetAztecOption(AZ_poly_ord, p.get("Polynomial Order", 3));
00394 }
00395 else if (aztecPreconditioner == "Least-squares Polynomial") {
00396 aztec.SetAztecOption(AZ_precond, AZ_ls);
00397 aztec.SetAztecOption(AZ_poly_ord, p.get("Polynomial Order", 3));
00398 }
00399 else {
00400 string errorMessage = "\"Aztec Preconditioner\" parameter is invalid!";
00401 throwError("setAztecOptions", errorMessage);
00402 }
00403
00404 }
00405 else
00406 aztec.SetAztecOption(AZ_precond, AZ_none);
00407
00408
00409
00410 string rcmReordering = p.get("RCM Reordering", "Disabled");
00411 if (rcmReordering == "Enabled")
00412 aztec.SetAztecOption(AZ_reorder, 1);
00413 else if (rcmReordering == "Disabled")
00414 aztec.SetAztecOption(AZ_reorder, 0);
00415 else {
00416 string errorMessage = "\"RCM Reordering\" parameter is invalid!";
00417 throwError("setAztecOptions", errorMessage);
00418 }
00419
00420
00421 string orthog = p.get("Orthogonalization", "Classical");
00422 if (orthog == "Classical")
00423 aztec.SetAztecOption(AZ_orthog, AZ_classic);
00424 else if (orthog == "Modified")
00425 aztec.SetAztecOption(AZ_orthog, AZ_modified);
00426 else {
00427 string errorMessage = "\"Orthogonalization\" parameter is invalid!";
00428 throwError("setAztecOptions()", errorMessage);
00429 }
00430
00431
00432 aztec.SetAztecOption(AZ_kspace,
00433 p.get("Size of Krylov Subspace", 300));
00434
00435
00436 string convCriteria = p.get("Convergence Test", "r0");
00437 if (convCriteria == "r0")
00438 aztec.SetAztecOption(AZ_conv, AZ_r0);
00439 else if (convCriteria == "rhs")
00440 aztec.SetAztecOption(AZ_conv, AZ_rhs);
00441 else if (convCriteria == "Anorm")
00442 aztec.SetAztecOption(AZ_conv, AZ_Anorm);
00443 else if (convCriteria == "no scaling")
00444 aztec.SetAztecOption(AZ_conv, AZ_noscaled);
00445 else if (convCriteria == "sol")
00446 aztec.SetAztecOption(AZ_conv, AZ_sol);
00447 else {
00448 string errorMessage = "\"Convergence Test\" parameter is invalid!";
00449 throwError("setAztecOptions()", errorMessage);
00450 }
00451
00452
00453 if (p.isParameter("Ill-Conditioning Threshold")) {
00454 aztec.SetAztecParam(AZ_ill_cond_thresh,
00455 p.get("Ill-Conditioning Threshold", 1.0e+11));
00456 }
00457
00458
00459 if (utils.isPrintType(Utils::LinearSolverDetails))
00460 aztec.SetAztecOption(AZ_output, p.get("Output Frequency",
00461 AZ_last));
00462 else
00463 aztec.SetAztecOption(AZ_output, p.get("Output Frequency", 0));
00464
00465
00466 if (utils.isPrintType(Utils::Debug)) {
00467
00468 }
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 return;
00482 }
00483
00484
00485 bool NOX::Epetra::LinearSystemAztecOO::createJacobianOperator(
00486 Teuchos::ParameterList& printParams,
00487 Teuchos::ParameterList& lsParams,
00488 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
00489 const NOX::Epetra::Vector& cloneVector)
00490 {
00491 string choice = lsParams.get("Jacobian Operator", "Matrix-Free");
00492
00493 if (choice == "Matrix-Free") {
00494 jacPtr =
00495 Teuchos::rcp(new MatrixFree(printParams, iReq, cloneVector));
00496 jacInterfacePtr =
00497 Teuchos::rcp_dynamic_cast<NOX::Epetra::Interface::Jacobian>(jacPtr);
00498 jacType = EpetraOperator;
00499 }
00500 else if (choice == "Finite Difference") {
00501 jacPtr =
00502 Teuchos::rcp(new FiniteDifference(printParams, iReq, cloneVector));
00503 jacInterfacePtr =
00504 Teuchos::rcp_dynamic_cast<NOX::Epetra::Interface::Jacobian>(jacPtr);
00505 jacType = EpetraRowMatrix;
00506 }
00507 else
00508 throwError("createJacobianOperator",
00509 "The specified value for parameter \" Jacobian Operator\" is not valid");
00510
00511 return true;
00512 }
00513
00514
00515 bool NOX::Epetra::LinearSystemAztecOO::createPrecOperator(
00516 Teuchos::ParameterList& printParams,
00517 Teuchos::ParameterList& lsParams,
00518 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
00519 const NOX::Epetra::Vector& cloneVector)
00520 {
00521 string choice = lsParams.get("Preconditioner Operator",
00522 "Use Jacobian");
00523
00524 if (choice == "Use Jacobian") {
00525 precType = jacType;
00526 precMatrixSource = UseJacobian;
00527 }
00528 else if (choice == "Finite Difference") {
00529 precPtr =
00530 Teuchos::rcp(new FiniteDifference(printParams, iReq, cloneVector));
00531 precInterfacePtr = Teuchos::rcp_dynamic_cast
00532 <NOX::Epetra::Interface::Preconditioner>(precPtr);
00533 precType = EpetraRowMatrix;
00534 precMatrixSource = SeparateMatrix;
00535 }
00536 else
00537 throwError("createPreconditionerOperator",
00538 "The value for the parameter \" Preconditioner Operator\" is not valid");
00539
00540 return true;
00541 }
00542
00543
00544 bool NOX::Epetra::LinearSystemAztecOO::
00545 applyJacobian(const NOX::Epetra::Vector& input,
00546 NOX::Epetra::Vector& result) const
00547 {
00548 jacPtr->SetUseTranspose(false);
00549 int status = jacPtr->Apply(input.getEpetraVector(),
00550 result.getEpetraVector());
00551 return (status == 0);
00552 }
00553
00554
00555 bool NOX::Epetra::LinearSystemAztecOO::
00556 applyJacobianTranspose(const NOX::Epetra::Vector& input,
00557 NOX::Epetra::Vector& result) const
00558 {
00559
00560 jacPtr->SetUseTranspose(true);
00561 int status = jacPtr->Apply(input.getEpetraVector(),
00562 result.getEpetraVector());
00563 jacPtr->SetUseTranspose(false);
00564
00565 return (status == 0);
00566 }
00567
00568
00569 bool NOX::Epetra::LinearSystemAztecOO::
00570 applyJacobianInverse(Teuchos::ParameterList &p,
00571 const NOX::Epetra::Vector& input,
00572 NOX::Epetra::Vector& result)
00573 {
00574
00575 if ( p.get("Use Preconditioner as Solver", false) )
00576 return applyRightPreconditioning(false, p, input, result);
00577
00578
00579 double startTime = timer.WallTime();
00580
00581
00582
00583
00584 NOX::Epetra::Vector& nonConstInput = const_cast<NOX::Epetra::Vector&>(input);
00585
00586
00587 if (zeroInitialGuess)
00588 result.init(0.0);
00589
00590
00591 Epetra_LinearProblem Problem(jacPtr.get(),
00592 &(result.getEpetraVector()),
00593 &(nonConstInput.getEpetraVector()));
00594
00595
00596
00597
00598
00599 this->setAztecOOJacobian();
00600 this->setAztecOOPreconditioner();
00601 aztecSolverPtr->SetLHS(&(result.getEpetraVector()));
00602 aztecSolverPtr->SetRHS(&(nonConstInput.getEpetraVector()));
00603
00604
00605 if ( !Teuchos::is_null(scaling) ) {
00606
00607 if ( !manualScaling )
00608 scaling->computeScaling(Problem);
00609
00610 scaling->scaleLinearSystem(Problem);
00611
00612 if (utils.isPrintType(Utils::Details)) {
00613 utils.out() << *scaling << endl;
00614 }
00615 }
00616
00617
00618
00619 #ifdef HAVE_NOX_DEBUG
00620 #ifdef HAVE_NOX_EPETRAEXT
00621
00622 ++linearSolveCount;
00623 std::ostringstream iterationNumber;
00624 iterationNumber << linearSolveCount;
00625
00626 std::string prefixName = p.get("Write Linear System File Prefix",
00627 "NOX_LinSys");
00628 std::string postfixName = iterationNumber.str();
00629 postfixName += ".mm";
00630
00631 if (p.get("Write Linear System", false)) {
00632
00633 std::string mapFileName = prefixName + "_Map_" + postfixName;
00634 std::string jacFileName = prefixName + "_Jacobian_" + postfixName;
00635 std::string rhsFileName = prefixName + "_RHS_" + postfixName;
00636
00637 Epetra_RowMatrix* printMatrix = NULL;
00638 printMatrix = dynamic_cast<Epetra_RowMatrix*>(jacPtr.get());
00639
00640 if (printMatrix == NULL) {
00641 cout << "Error: NOX::Epetra::LinearSystemAztecOO::applyJacobianInverse() - "
00642 << "Could not cast the Jacobian operator to an Epetra_RowMatrix!"
00643 << "Please set the \"Write Linear System\" parameter to false."
00644 << endl;
00645 throw "NOX Error";
00646 }
00647
00648 EpetraExt::BlockMapToMatrixMarketFile(mapFileName.c_str(),
00649 printMatrix->RowMatrixRowMap());
00650 EpetraExt::RowMatrixToMatrixMarketFile(jacFileName.c_str(), *printMatrix,
00651 "test matrix", "Jacobian XXX");
00652 EpetraExt::MultiVectorToMatrixMarketFile(rhsFileName.c_str(),
00653 nonConstInput.getEpetraVector());
00654
00655 }
00656 #endif
00657 #endif
00658
00659
00660 if (!isPrecConstructed && (precAlgorithm != None_)) {
00661 throwError("applyJacobianInverse",
00662 "Preconditioner is not constructed! Call createPreconditioner() first.");
00663 }
00664
00665
00666 int maxit = p.get("Max Iterations", 400);
00667 double tol = p.get("Tolerance", 1.0e-6);
00668
00669 int aztecStatus = -1;
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 if (precAlgorithm == AztecOO_ &&
00681 precReusePolicy == PRPT_REUSE)
00682 aztecSolverPtr->SetAztecOption(AZ_pre_calc, AZ_calc);
00683
00684 aztecStatus = aztecSolverPtr->Iterate(maxit, tol);
00685
00686
00687 if ( !Teuchos::is_null(scaling) )
00688 scaling->unscaleLinearSystem(Problem);
00689
00690
00691 if (outputSolveDetails) {
00692 Teuchos::ParameterList& outputList = p.sublist("Output");
00693 int prevLinIters =
00694 outputList.get("Total Number of Linear Iterations", 0);
00695 int curLinIters = 0;
00696 double achievedTol = -1.0;
00697 curLinIters = aztecSolverPtr->NumIters();
00698 achievedTol = aztecSolverPtr->ScaledResidual();
00699
00700 outputList.set("Number of Linear Iterations", curLinIters);
00701 outputList.set("Total Number of Linear Iterations",
00702 (prevLinIters + curLinIters));
00703 outputList.set("Achieved Tolerance", achievedTol);
00704 }
00705
00706
00707 #ifdef HAVE_NOX_DEBUG
00708 #ifdef HAVE_NOX_EPETRAEXT
00709 if (p.get("Write Linear System", false)) {
00710 std::string lhsFileName = prefixName + "_LHS_" + postfixName;
00711 EpetraExt::MultiVectorToMatrixMarketFile(lhsFileName.c_str(),
00712 result.getEpetraVector());
00713 }
00714 #endif
00715 #endif
00716
00717 double endTime = timer.WallTime();
00718 timeApplyJacbianInverse += (endTime - startTime);
00719
00720 if (aztecStatus != 0)
00721 return false;
00722
00723 return true;
00724 }
00725
00726
00727 bool NOX::Epetra::LinearSystemAztecOO::
00728 applyRightPreconditioning(bool useTranspose,
00729 Teuchos::ParameterList& params,
00730 const NOX::Epetra::Vector& input,
00731 NOX::Epetra::Vector& result) const
00732 {
00733 int errorCode = 1;
00734
00735
00736 if (!isPrecConstructed) {
00737 throwError("applyRightPreconditioning",
00738 "Preconditioner is not constructed! Call createPreconditioner() first.");
00739 }
00740
00741 if (precAlgorithm == None_) {
00742 if (&result != &input)
00743 result = input;
00744 return true;
00745 }
00746 else if (precAlgorithm == AztecOO_) {
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756 tmpVectorPtr->init(0.0);
00757
00758
00759 aztecSolverPtr->SetAztecOption(AZ_output,AZ_none);
00760
00761
00762 int numIters = params.get("AztecOO Preconditioner Iterations", 1);
00763
00764 AztecOO_Operator prec(aztecSolverPtr.get(), numIters);
00765
00766 errorCode = prec.ApplyInverse(input.getEpetraVector(),
00767 result.getEpetraVector());
00768 }
00769 else if (precAlgorithm == Ifpack_){
00770
00771 if (useTranspose)
00772 ifpackPreconditionerPtr->SetUseTranspose(useTranspose);
00773
00774 errorCode = ifpackPreconditionerPtr->ApplyInverse(input.getEpetraVector(),
00775 result.getEpetraVector());
00776
00777 if (useTranspose)
00778 ifpackPreconditionerPtr->SetUseTranspose(false);
00779
00780 }
00781 else if (precAlgorithm == NewIfpack_){
00782
00783 if (useTranspose)
00784 newIfpackPreconditionerPtr->SetUseTranspose(useTranspose);
00785
00786 errorCode = newIfpackPreconditionerPtr->ApplyInverse(input.getEpetraVector(),
00787 result.getEpetraVector());
00788
00789 if (useTranspose)
00790 newIfpackPreconditionerPtr->SetUseTranspose(false);
00791
00792 }
00793 else if (precAlgorithm == UserDefined_) {
00794
00795 if (useTranspose)
00796 precPtr->SetUseTranspose(true);
00797
00798 errorCode = precPtr->ApplyInverse(input.getEpetraVector(),
00799 result.getEpetraVector());
00800 if (useTranspose)
00801 precPtr->SetUseTranspose(false);
00802
00803 }
00804 else
00805 throwError("applyRightPreconditioning",
00806 "Parameter \"preconditioner\" is not vaild for this method");
00807
00808 if (errorCode != 0) {
00809 std::string msg = "Error - NOX::Epetra::LinearSystemAztecOO::applyRightPreconditioning() - A non-zero error code has been returned from the preconditioner.";
00810 if (throwErrorOnPrecFailure) {
00811 TEST_FOR_EXCEPTION(true, std::logic_error, msg);
00812 }
00813 else {
00814 if (utils.isPrintType(NOX::Utils::Warning))
00815 utils.out() << msg << endl;
00816 }
00817 return false;
00818 }
00819
00820 return true;
00821 }
00822
00823
00824 bool NOX::Epetra::LinearSystemAztecOO::checkPreconditionerValidity()
00825 {
00826 if (precAlgorithm == None_)
00827 return true;
00828
00829
00830 else if ((precAlgorithm == AztecOO_) || (precAlgorithm == Ifpack_) ||
00831 (precAlgorithm == NewIfpack_) || (precAlgorithm == ML_)) {
00832
00833
00834 if ((precType == EpetraRowMatrix) ||
00835 (precType == EpetraVbrMatrix) ||
00836 (precType == EpetraCrsMatrix)) {
00837 return true;
00838 }
00839 else {
00840 throwError("checkPreconditionerValidity", "Preconditioner must be an Epetra_RowMatrix derived object for AztecOO and Ifpack preconditioners");
00841 }
00842
00843 }
00844 else if (precAlgorithm == UserDefined_){
00845
00846
00847 if (Teuchos::is_null(precPtr)) {
00848 throwError("checkPreconditionerValidity", "Preconditioiner is NULL!");
00849 }
00850 return true;
00851 }
00852
00853 return true;
00854 }
00855
00856
00857 bool NOX::Epetra::LinearSystemAztecOO::
00858 createPreconditioner(const NOX::Epetra::Vector& x, Teuchos::ParameterList& p,
00859 bool recomputeGraph) const
00860 {
00861 double startTime = timer.WallTime();
00862
00863 if (utils.isPrintType(Utils::LinearSolverDetails))
00864 utils.out() << "\n Creating a new preconditioner" << endl;;
00865
00866 if (precAlgorithm == None_) {
00867 return true;
00868 }
00869
00870
00871 Epetra_LinearProblem Problem(jacPtr.get(),
00872 &(tmpVectorPtr->getEpetraVector()),
00873 &(tmpVectorPtr->getEpetraVector()));
00874 if ( !Teuchos::is_null(scaling) ) {
00875 if (!manualScaling)
00876 scaling->computeScaling(Problem);
00877
00878 scaling->scaleLinearSystem(Problem);
00879 }
00880
00881
00882 if (precAlgorithm == AztecOO_) {
00883
00884
00885
00886
00887
00888 if (precMatrixSource == UseJacobian) {
00889
00890
00891 this->setAztecOOJacobian();
00892 aztecSolverPtr->SetPrecMatrix(dynamic_cast<Epetra_RowMatrix*>(jacPtr.get()));
00893 aztecSolverPtr->ConstructPreconditioner(conditionNumberEstimate);
00894 }
00895 else if (precMatrixSource == SeparateMatrix) {
00896 Epetra_RowMatrix& precMatrix = dynamic_cast<Epetra_RowMatrix&>(*precPtr);
00897 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
00898 *precPtr, &p);
00899 this->setAztecOOJacobian();
00900 aztecSolverPtr->SetPrecMatrix(&precMatrix);
00901 aztecSolverPtr->ConstructPreconditioner(conditionNumberEstimate);
00902 }
00903
00904 }
00905 else if (precAlgorithm == Ifpack_) {
00906
00907 if (precMatrixSource == UseJacobian) {
00908 createIfpackPreconditioner(p);
00909 solvePrecOpPtr = ifpackPreconditionerPtr;
00910
00911 }
00912 else if (precMatrixSource == SeparateMatrix) {
00913
00914 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
00915 *precPtr, &p);
00916 createIfpackPreconditioner(p);
00917 solvePrecOpPtr = ifpackPreconditionerPtr;
00918 }
00919
00920 }
00921 else if (precAlgorithm == NewIfpack_) {
00922
00923 if (precMatrixSource == UseJacobian) {
00924 createNewIfpackPreconditioner(p);
00925 solvePrecOpPtr = newIfpackPreconditionerPtr;
00926 }
00927 else if (precMatrixSource == SeparateMatrix) {
00928
00929 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
00930 *precPtr, &p);
00931 createNewIfpackPreconditioner(p);
00932 solvePrecOpPtr = newIfpackPreconditionerPtr;
00933 }
00934
00935 }
00936 #ifdef HAVE_NOX_ML_EPETRA
00937 else if (precAlgorithm == ML_) {
00938
00939 if (precMatrixSource == UseJacobian) {
00940 createMLPreconditioner(p);
00941 solvePrecOpPtr = MLPreconditionerPtr;
00942 }
00943 else if (precMatrixSource == SeparateMatrix) {
00944
00945 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
00946 *precPtr, &p);
00947 createMLPreconditioner(p);
00948 solvePrecOpPtr = MLPreconditionerPtr;
00949 }
00950
00951 }
00952 #endif
00953 else if (precAlgorithm == UserDefined_) {
00954
00955 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
00956 *precPtr, &p);
00957 solvePrecOpPtr = precPtr;
00958
00959 }
00960
00961 isPrecConstructed = true;
00962
00963
00964 if ( !Teuchos::is_null(scaling) )
00965 scaling->unscaleLinearSystem(Problem);
00966
00967 double endTime = timer.WallTime();
00968 timeCreatePreconditioner += (endTime - startTime);
00969
00970 if (utils.isPrintType(Utils::LinearSolverDetails))
00971 utils.out() << "\n Time required to create precondtioner : "
00972 << (endTime - startTime) << " (sec.)" << endl;;
00973
00974 return true;
00975 }
00976
00977
00978 bool NOX::Epetra::LinearSystemAztecOO::
00979 createIfpackPreconditioner(Teuchos::ParameterList& p) const
00980 {
00981
00982
00983 if ( !Teuchos::is_null(ifpackGraphPtr) )
00984 throwError("createIfpackPreconditioner", "Ifpack Graph NOT NULL");
00985 if ( !Teuchos::is_null(ifpackPreconditionerPtr))
00986 throwError("createIfpackPreconditioner", "Ifpack Prec NOT NULL");
00987
00988 if (utils.isPrintType(Utils::Debug))
00989 utils.out() << "NOX::Epetra::LinearSolverAztecOO : createIfpackPrecon - \n"
00990 << " using Fill Factor --> " << p.get("Fill Factor", 1)
00991 << endl;
00992
00993
00994 if (precType == EpetraVbrMatrix) {
00995
00996 Epetra_VbrMatrix* vbr = 0;
00997
00998 if (precMatrixSource == UseJacobian)
00999 vbr = dynamic_cast<Epetra_VbrMatrix*>(jacPtr.get());
01000 else if (precMatrixSource == SeparateMatrix)
01001 vbr = dynamic_cast<Epetra_VbrMatrix*>(precPtr.get());
01002
01003 if (vbr == 0)
01004 throwError("createIfpackPreconditioner",
01005 "Dynamic cast to VBR Matrix failed!");
01006
01007 ifpackGraphPtr =
01008 Teuchos::rcp(new Ifpack_IlukGraph(vbr->Graph(),
01009 p.get("Fill Factor", 1),
01010 p.get("Overlap", 0)));
01011 ifpackGraphPtr->ConstructFilledGraph();
01012 ifpackPreconditionerPtr =
01013 Teuchos::rcp(new Ifpack_CrsRiluk(*ifpackGraphPtr));
01014 ifpackPreconditionerPtr->SetAbsoluteThreshold(p.get("Absolute Threshold", 0.0));
01015 ifpackPreconditionerPtr->SetRelativeThreshold(p.get("Relative Threshold", 1.0));
01016 int err = ifpackPreconditionerPtr->InitValues(*vbr);
01017 if (err != 0)
01018 precError(err, "createIfpackPreconditioner()", "Ifpack", "InitValues");
01019
01020 err = ifpackPreconditionerPtr->Factor();
01021 if (err != 0)
01022 precError(err, "createIfpackPreconditioner()", "Ifpack", "Factor");
01023
01024 }
01025
01026
01027 else if (precType == EpetraCrsMatrix) {
01028
01029 Epetra_CrsMatrix* crs = 0;
01030
01031 if (precMatrixSource == UseJacobian)
01032 crs = dynamic_cast<Epetra_CrsMatrix*>(jacPtr.get());
01033 else if (precMatrixSource == SeparateMatrix)
01034 crs = dynamic_cast<Epetra_CrsMatrix*>(precPtr.get());
01035
01036 if (crs == 0)
01037 throwError("createIfpackPreconditioner",
01038 "Dynamic cast to CRS Matrix failed!");
01039
01040 ifpackGraphPtr =
01041 Teuchos::rcp(new Ifpack_IlukGraph(crs->Graph(),
01042 p.get("Fill Factor", 1),
01043 p.get("Overlap", 0)));
01044 ifpackGraphPtr->ConstructFilledGraph();
01045 ifpackPreconditionerPtr =
01046 Teuchos::rcp(new Ifpack_CrsRiluk(*ifpackGraphPtr));
01047 ifpackPreconditionerPtr->SetAbsoluteThreshold(p.get("Absolute Threshold", 0.0));
01048 ifpackPreconditionerPtr->SetRelativeThreshold(p.get("Relative Threshold", 1.0));
01049 int err = ifpackPreconditionerPtr->InitValues(*crs);
01050 if (err != 0)
01051 precError(err, "createIfpackPreconditioner()", "Ifpack", "InitValues");
01052
01053 err = ifpackPreconditionerPtr->Factor();
01054 if (err != 0)
01055 precError(err, "createIfpackPreconditioner()", "Ifpack", "Factor");
01056
01057 }
01058
01059
01060 else if (precType == EpetraRowMatrix) {
01061
01062
01063
01064 Epetra_CrsMatrix* crs = 0;
01065 NOX::Epetra::FiniteDifference* FDoperator = 0;
01066
01067 if (precMatrixSource == UseJacobian)
01068 FDoperator = dynamic_cast<NOX::Epetra::FiniteDifference*>(jacPtr.get());
01069 else if (precMatrixSource == SeparateMatrix)
01070 FDoperator = dynamic_cast<NOX::Epetra::FiniteDifference*>(precPtr.get());
01071
01072 if (FDoperator != 0)
01073 crs = &(FDoperator->getUnderlyingMatrix());
01074
01075 if (crs == 0)
01076 throwError("createIfpackPreconditioner",
01077 "FiniteDifference: Underlying matrix NOT CRS Matrix!");
01078
01079 ifpackGraphPtr =
01080 Teuchos::rcp(new Ifpack_IlukGraph(crs->Graph(),
01081 p.get("Fill Factor", 1),
01082 p.get("Overlap", 0)));
01083 ifpackGraphPtr->ConstructFilledGraph();
01084 ifpackPreconditionerPtr =
01085 Teuchos::rcp(new Ifpack_CrsRiluk(*ifpackGraphPtr));
01086 ifpackPreconditionerPtr->SetAbsoluteThreshold(p.get("Absolute Threshold", 0.0));
01087 ifpackPreconditionerPtr->SetRelativeThreshold(p.get("Relative Threshold", 1.0));
01088 int err = ifpackPreconditionerPtr->InitValues(*crs);
01089 if (err != 0)
01090 precError(err, "createIfpackPreconditioner()", "Ifpack", "InitValues");
01091
01092 err = ifpackPreconditionerPtr->Factor();
01093 if (err != 0)
01094 precError(err, "createIfpackPreconditioner()", "Ifpack", "InitValues");
01095
01096 }
01097
01098 else {
01099
01100
01101
01102
01103 throwError("createIfpackPreconditioner",
01104 "Ifpack requires a CRS or VBR matrix!");
01105 return false;
01106 }
01107
01108
01109 if (utils.isPrintType(Utils::Details)) {
01110 double kappa;
01111 ifpackPreconditionerPtr->Condest(false, kappa);
01112 utils.out() <<
01113 "\n Condition number estimate of preconditioner is " <<
01114 utils.sciformat(kappa) << std::endl;
01115 }
01116
01117 return true;
01118 }
01119
01120
01121 bool NOX::Epetra::LinearSystemAztecOO::
01122 createNewIfpackPreconditioner(Teuchos::ParameterList& p) const
01123 {
01124
01125
01126 if ( !Teuchos::is_null(newIfpackPreconditionerPtr) )
01127 throwError("createNewIfpackPreconditioner", "Ifpack Prec NOT NULL");
01128
01129
01130 Teuchos::ParameterList& teuchosParams =
01131 p.sublist("Ifpack");
01132
01133 if (utils.isPrintType(Utils::Debug)) {
01134 utils.out() << "NOX::Epetra::LinearSolverAztecOO : createNewIfpackPrecon - \n"
01135 << " using Teuchos parameter : " << endl;
01136 teuchosParams.print(utils.out());
01137 }
01138
01139 Ifpack Factory;
01140
01141
01142 if (precType == EpetraVbrMatrix || precType == EpetraCrsMatrix ||
01143 precType == EpetraRowMatrix) {
01144
01145 Epetra_RowMatrix* row = 0;
01146
01147 if (precMatrixSource == UseJacobian)
01148 row = dynamic_cast<Epetra_RowMatrix*>(jacPtr.get());
01149 else if (precMatrixSource == SeparateMatrix)
01150 row = dynamic_cast<Epetra_RowMatrix*>(precPtr.get());
01151
01152 if (row == 0)
01153 throwError("createIfpackPreconditioner",
01154 "Dynamic cast to Row Matrix failed!");
01155
01156 newIfpackPreconditionerPtr = Teuchos::rcp(Factory.Create(
01157 p.get("Ifpack Preconditioner", "ILU"),
01158 row,
01159 p.get("Overlap", 0) ));
01160 newIfpackPreconditionerPtr->SetParameters(teuchosParams);
01161 int err = newIfpackPreconditionerPtr->Initialize();
01162 if (err != 0)
01163 precError(err, "createNewIfpackPreconditioner()", "Ifpack", "Initialize");
01164
01165
01166 err = newIfpackPreconditionerPtr->Compute();
01167 if (err != 0)
01168 precError(err, "createNewIfpackPreconditioner()", "Ifpack", "Compute");
01169
01170 return true;
01171 }
01172
01173
01174
01175
01176 throwError("createNewIfpackPreconditioner",
01177 "Ifpack requires a CRS or VBR matrix!");
01178
01179 return false;
01180 }
01181
01182 #ifdef HAVE_NOX_ML_EPETRA
01183
01184 bool NOX::Epetra::LinearSystemAztecOO::
01185 createMLPreconditioner(Teuchos::ParameterList& p) const
01186 {
01187
01188
01189 if( !Teuchos::is_null(MLPreconditionerPtr) )
01190 throwError("createMLPreconditioner", "ML Prec NOT NULL");
01191
01192
01193 Teuchos::ParameterList& teuchosParams = p.sublist("ML");
01194
01195 if (utils.isPrintType(Utils::Debug))
01196 {
01197 utils.out() << "NOX::Epetra::LinearSolverAztecOO : createMLPreconditioner - \n"
01198 << " using Teuchos parameter : " << endl;
01199 teuchosParams.print(utils.out());
01200 }
01201
01202
01203 if (precType == EpetraVbrMatrix || precType == EpetraCrsMatrix ||
01204 precType == EpetraRowMatrix) {
01205
01206 Epetra_RowMatrix* rowMatrix = 0;
01207
01208 if (precMatrixSource == UseJacobian)
01209 rowMatrix = dynamic_cast<Epetra_RowMatrix*>(jacPtr.get());
01210 else if (precMatrixSource == SeparateMatrix)
01211 rowMatrix = dynamic_cast<Epetra_RowMatrix*>(precPtr.get());
01212
01213 if (rowMatrix == 0)
01214 throwError("createMLPreconditioner",
01215 "Dynamic cast to Epetra_RowMatrix failed!");
01216
01217 MLPreconditionerPtr = Teuchos::rcp( new ML_Epetra::MultiLevelPreconditioner(
01218 *rowMatrix, teuchosParams, true) );
01219 return true;
01220
01221 }
01222
01223
01224
01225
01226 throwError("createMLPreconditioner",
01227 "ML requires an Epetra_RowMatrix !");
01228
01229 return false;
01230 }
01231 #endif
01232
01233
01234 bool NOX::Epetra::LinearSystemAztecOO::
01235 recomputePreconditioner(const NOX::Epetra::Vector& x,
01236 Teuchos::ParameterList& linearSolverParams) const
01237 {
01238 if (utils.isPrintType(Utils::LinearSolverDetails)) {
01239 utils.out() << "\n Recomputing preconditioner" << endl;
01240 }
01241
01242 if (precAlgorithm == None_) {
01243 return true;
01244 }
01245 else if (precAlgorithm == AztecOO_) {
01246
01247 if (precMatrixSource == SeparateMatrix)
01248 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
01249 *precPtr, &linearSolverParams);
01250
01251 }
01252 else if (precAlgorithm == Ifpack_) {
01253
01254 throwError("recomputePreconditioner",
01255 "The preconditioner can only be \"Recomputed\" when using the \"NewIfpack\" or \"ML\" Preconditioner types!");
01256 }
01257 else if (precAlgorithm == NewIfpack_) {
01258
01259 if (precMatrixSource == SeparateMatrix)
01260 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
01261 *precPtr, &linearSolverParams);
01262
01263 int err = newIfpackPreconditionerPtr->Compute();
01264 if (err != 0)
01265 precError(err, "recomputePreconditioner", "Ifpack", "Compute");
01266
01267 }
01268 #ifdef HAVE_NOX_ML_EPETRA
01269 else if (precAlgorithm == ML_) {
01270
01271 if (precMatrixSource == SeparateMatrix)
01272 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
01273 *precPtr, &linearSolverParams);
01274
01275 int err = MLPreconditionerPtr->ReComputePreconditioner();
01276
01277 if (err != 0)
01278 precError(err, "recomputePreconditioner", "ML",
01279 "ReComputePreconditioner");
01280
01281
01282 }
01283 #endif
01284 else if (precAlgorithm == UserDefined_) {
01285
01286 if (precMatrixSource == SeparateMatrix)
01287 precInterfacePtr->computePreconditioner(x.getEpetraVector(),
01288 *precPtr, &linearSolverParams);
01289
01290 }
01291
01292 return true;
01293 }
01294
01295
01296 bool NOX::Epetra::LinearSystemAztecOO::destroyPreconditioner() const
01297 {
01298 if (isPrecConstructed) {
01299
01300 if (precAlgorithm == AztecOO_) {
01301 aztecSolverPtr->DestroyPreconditioner();
01302 }
01303
01304 else if (precAlgorithm == Ifpack_) {
01305 ifpackPreconditionerPtr = Teuchos::null;
01306 ifpackGraphPtr = Teuchos::null;
01307 }
01308 else if (precAlgorithm == NewIfpack_) {
01309 newIfpackPreconditionerPtr = Teuchos::null;
01310 }
01311 #ifdef HAVE_NOX_ML_EPETRA
01312 else if (precAlgorithm == ML_)
01313 MLPreconditionerPtr = Teuchos::null;
01314 #endif
01315 if (utils.isPrintType(Utils::LinearSolverDetails)) {
01316 utils.out() << "\n Destroying preconditioner" << endl;
01317 }
01318 }
01319 isPrecConstructed = false;
01320 solvePrecOpPtr = Teuchos::null;
01321 return true;
01322 }
01323
01324
01325 NOX::Epetra::LinearSystemAztecOO::OperatorType
01326 NOX::Epetra::LinearSystemAztecOO::getOperatorType(const Epetra_Operator& Op)
01327 {
01328
01329
01330
01331
01332 const Epetra_Operator* testOperator = 0;
01333
01334
01335 testOperator = dynamic_cast<const Epetra_CrsMatrix*>(&Op);
01336 if (testOperator != 0)
01337 return EpetraCrsMatrix;
01338
01339
01340 testOperator = dynamic_cast<const Epetra_VbrMatrix*>(&Op);
01341 if (testOperator != 0)
01342 return EpetraVbrMatrix;
01343
01344
01345 testOperator = dynamic_cast<const Epetra_RowMatrix*>(&Op);
01346 if (testOperator != 0)
01347 return EpetraRowMatrix;
01348
01349
01350 return EpetraOperator;
01351 }
01352
01353
01354 bool NOX::Epetra::LinearSystemAztecOO::
01355 computeJacobian(const NOX::Epetra::Vector& x)
01356 {
01357 bool success = jacInterfacePtr->computeJacobian(x.getEpetraVector(),
01358 *jacPtr);
01359 return success;
01360 }
01361
01362
01363 Teuchos::RCP<NOX::Epetra::Scaling>
01364 NOX::Epetra::LinearSystemAztecOO::getScaling()
01365 {
01366 return scaling;
01367 }
01368
01369
01370 void NOX::Epetra::LinearSystemAztecOO::
01371 resetScaling(const Teuchos::RCP<NOX::Epetra::Scaling>& scalingObject)
01372 {
01373 scaling = scalingObject;
01374 return;
01375 }
01376
01377
01378 void NOX::Epetra::LinearSystemAztecOO::
01379 throwError(const string& functionName, const string& errorMsg) const
01380 {
01381 if (utils.isPrintType(Utils::Error)) {
01382 utils.out() << "NOX::Epetra::LinearSystemAztecOO::" << functionName
01383 << " - " << errorMsg << endl;
01384 }
01385 throw "NOX Error";
01386 }
01387
01388
01389 Teuchos::RCP<const NOX::Epetra::Interface::Jacobian>
01390 NOX::Epetra::LinearSystemAztecOO::getJacobianInterface() const
01391 {
01392 return jacInterfacePtr;
01393 }
01394
01395
01396 Teuchos::RCP<const NOX::Epetra::Interface::Preconditioner>
01397 NOX::Epetra::LinearSystemAztecOO::getPrecInterface() const
01398 {
01399 return precInterfacePtr;
01400 }
01401
01402
01403 bool
01404 NOX::Epetra::LinearSystemAztecOO::hasPreconditioner() const
01405 {
01406 return (precAlgorithm != None_);
01407 }
01408
01409
01410 bool
01411 NOX::Epetra::LinearSystemAztecOO::isPreconditionerConstructed() const
01412 {
01413 return isPrecConstructed;
01414 }
01415
01416
01417 Teuchos::RCP<const Epetra_Operator>
01418 NOX::Epetra::LinearSystemAztecOO::getJacobianOperator() const
01419 {
01420 return jacPtr;
01421 }
01422
01423
01424 Teuchos::RCP<Epetra_Operator>
01425 NOX::Epetra::LinearSystemAztecOO::getJacobianOperator()
01426 {
01427 return jacPtr;
01428 }
01429
01430
01431 Teuchos::RCP<const Epetra_Operator>
01432 NOX::Epetra::LinearSystemAztecOO::getPrecOperator() const
01433 {
01434 return precPtr;
01435 }
01436
01437
01438 Teuchos::RCP<const Epetra_Operator>
01439 NOX::Epetra::LinearSystemAztecOO::getGeneratedPrecOperator() const
01440 {
01441 return solvePrecOpPtr;
01442 }
01443
01444
01445 Teuchos::RCP<Epetra_Operator>
01446 NOX::Epetra::LinearSystemAztecOO::getGeneratedPrecOperator()
01447 {
01448 return solvePrecOpPtr;
01449 }
01450
01451
01452 NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
01453 NOX::Epetra::LinearSystemAztecOO::
01454 getPreconditionerPolicy(bool advanceReuseCounter)
01455 {
01456
01457 if (precReusePolicy == PRPT_REBUILD)
01458 return PRPT_REBUILD;
01459
01460 if (precReusePolicy == PRPT_RECOMPUTE) {
01461 if (isPrecConstructed)
01462 return PRPT_RECOMPUTE;
01463 else
01464 return PRPT_REBUILD;
01465 }
01466
01467
01468
01469
01470
01471 if (precReusePolicy == PRPT_REUSE) {
01472
01473
01474 if (!isPrecConstructed) {
01475 if (advanceReuseCounter)
01476 precQueryCounter++;
01477 return PRPT_REBUILD;
01478 }
01479
01480 if (utils.isPrintType(Utils::Details))
01481 if (advanceReuseCounter)
01482 utils.out() << "\n Preconditioner Reuse: Age of Prec --> "
01483 << precQueryCounter << " / " << maxAgeOfPrec << endl;
01484
01485
01486 if( maxAgeOfPrec == -2 ) {
01487 if (advanceReuseCounter)
01488 precQueryCounter++;
01489 return PRPT_REUSE;
01490 }
01491
01492
01493
01494 else if( maxAgeOfPrec == -1 ) {
01495 if (advanceReuseCounter)
01496 precQueryCounter++;
01497 maxAgeOfPrec = -2;
01498 return PRPT_REBUILD;
01499 }
01500
01501
01502 else {
01503 if( precQueryCounter == 0 || precQueryCounter >= maxAgeOfPrec ) {
01504 if (advanceReuseCounter)
01505 precQueryCounter = 1;
01506 return PRPT_REBUILD;
01507 }
01508 else {
01509 if (advanceReuseCounter)
01510 precQueryCounter++;
01511 return PRPT_REUSE;
01512 }
01513 }
01514 }
01515
01516 return PRPT_REBUILD;
01517 }
01518
01519
01520 double
01521 NOX::Epetra::LinearSystemAztecOO::getTimeCreatePreconditioner() const
01522 {
01523 return timeCreatePreconditioner;
01524 }
01525
01526
01527 double
01528 NOX::Epetra::LinearSystemAztecOO::getTimeApplyJacobianInverse() const
01529 {
01530 return timeApplyJacbianInverse;
01531 }
01532
01533
01534 void
01535 NOX::Epetra::LinearSystemAztecOO::setJacobianOperatorForSolve(
01536 const Teuchos::RCP<const Epetra_Operator>& solveJacOp)
01537 {
01538 jacPtr = Teuchos::rcp_const_cast<Epetra_Operator>(solveJacOp);
01539 jacType = getOperatorType(*solveJacOp);
01540 this->setAztecOOJacobian();
01541 }
01542
01543
01544 void
01545 NOX::Epetra::LinearSystemAztecOO::setPrecOperatorForSolve(
01546 const Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
01547 {
01548 solvePrecOpPtr = Teuchos::rcp_const_cast<Epetra_Operator>(solvePrecOp);
01549 this->setAztecOOPreconditioner();
01550 }
01551
01552
01553 void
01554 NOX::Epetra::LinearSystemAztecOO::setAztecOOJacobian() const
01555 {
01556 if ((jacType == EpetraRowMatrix) ||
01557 (jacType == EpetraVbrMatrix) ||
01558 (jacType == EpetraCrsMatrix)) {
01559 aztecSolverPtr->SetUserMatrix(dynamic_cast<Epetra_RowMatrix*>(jacPtr.get()), false);
01560 }
01561 else
01562 aztecSolverPtr->SetUserOperator(jacPtr.get());
01563 }
01564
01565
01566 void
01567 NOX::Epetra::LinearSystemAztecOO::setAztecOOPreconditioner() const
01568 {
01569 if ( !Teuchos::is_null(solvePrecOpPtr) && precAlgorithm != AztecOO_)
01570 aztecSolverPtr->SetPrecOperator(solvePrecOpPtr.get());
01571 }
01572
01573
01574 void NOX::Epetra::LinearSystemAztecOO::
01575 precError(int error_code,
01576 const std::string& nox_function,
01577 const std::string& prec_type,
01578 const std::string& prec_function) const
01579 {
01580 if (error_code != 0) {
01581
01582 std::ostringstream msg;
01583
01584 if (throwErrorOnPrecFailure)
01585 msg << "Error - ";
01586 else
01587 msg << "Warning - ";
01588
01589 msg << "NOX::Epetra::LinearSystemAztecOO::" << nox_function << " - The " << prec_type << " preconditioner has returned a nonzero error code of " << error_code << " for the function \"" << prec_function << "\". Please consult the preconditioner documentation for this error type.";
01590
01591 if (throwErrorOnPrecFailure) {
01592 TEST_FOR_EXCEPTION(true, std::logic_error, msg.str());
01593 }
01594 else
01595 utils.out() << msg.str() << endl;
01596 }
01597 }
01598
01599
01600