49#include "NOX_Epetra_Interface_Jacobian.H"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_TimeMonitor.hpp"
52#include "Teuchos_Assert.hpp"
54NOX::Epetra::LinearSystemMPBD::
56 Teuchos::ParameterList& printingParams,
57 Teuchos::ParameterList& linearSolverParams,
58 const Teuchos::RCP<NOX::Epetra::LinearSystem>& block_solver_,
59 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
60 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
61 const Teuchos::RCP<Epetra_Operator>& J,
62 const Teuchos::RCP<const Epetra_Map>& base_map_,
63 const Teuchos::RCP<NOX::Epetra::Scaling> s):
64 block_solver(block_solver_),
65 jacInterfacePtr(iJac),
70 mp_op = Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(J,
true);
71 block_ops = mp_op->getMPOps();
72 num_mp_blocks = block_ops->size();
74 std::string prec_strategy = linearSolverParams.get(
"Preconditioner Strategy",
76 if (prec_strategy ==
"Standard")
77 precStrategy = STANDARD;
78 else if (prec_strategy ==
"Mean")
80 else if (prec_strategy ==
"On the fly")
81 precStrategy = ON_THE_FLY;
83 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
84 "Invalid preconditioner strategy " << prec_strategy);
86 if (precStrategy == STANDARD)
87 precs.resize(num_mp_blocks);
90NOX::Epetra::LinearSystemMPBD::
95bool NOX::Epetra::LinearSystemMPBD::
96applyJacobian(
const NOX::Epetra::Vector& input,
97 NOX::Epetra::Vector& result)
const
99 mp_op->SetUseTranspose(
false);
100 int status = mp_op->Apply(input.getEpetraVector(), result.getEpetraVector());
102 return (status == 0);
105bool NOX::Epetra::LinearSystemMPBD::
106applyJacobianTranspose(
const NOX::Epetra::Vector& input,
107 NOX::Epetra::Vector& result)
const
109 mp_op->SetUseTranspose(
true);
110 int status = mp_op->Apply(input.getEpetraVector(), result.getEpetraVector());
111 mp_op->SetUseTranspose(
false);
113 return (status == 0);
116bool NOX::Epetra::LinearSystemMPBD::
117applyJacobianInverse(Teuchos::ParameterList ¶ms,
118 const NOX::Epetra::Vector &input,
119 NOX::Epetra::Vector &result)
121 TEUCHOS_FUNC_TIME_MONITOR(
"Total deterministic solve Time");
124 EpetraExt::BlockVector input_block(View, *base_map,
125 input.getEpetraVector());
126 EpetraExt::BlockVector result_block(View, *base_map,
127 result.getEpetraVector());
128 result_block.PutScalar(0.0);
131 Teuchos::ParameterList& block_solver_params =
132 params.sublist(
"Deterministic Solver Parameters");
135 bool final_status =
true;
137 for (
int i=0; i<num_mp_blocks; i++) {
138 NOX::Epetra::Vector nox_input(input_block.GetBlock(i),
139 NOX::Epetra::Vector::CreateView);
140 NOX::Epetra::Vector nox_result(result_block.GetBlock(i),
141 NOX::Epetra::Vector::CreateView);
143 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
145 if (precStrategy == STANDARD)
146 block_solver->setPrecOperatorForSolve(precs[i]);
147 else if (precStrategy == ON_THE_FLY) {
148 block_solver->createPreconditioner(*(prec_x->GetBlock(i)),
149 block_solver_params,
false);
152 status = block_solver->applyJacobianInverse(block_solver_params, nox_input,
154 final_status = final_status && status;
160bool NOX::Epetra::LinearSystemMPBD::
161applyRightPreconditioning(
bool useTranspose,
162 Teuchos::ParameterList& params,
163 const NOX::Epetra::Vector& input,
164 NOX::Epetra::Vector& result)
const
169Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemMPBD::
175void NOX::Epetra::LinearSystemMPBD::
176resetScaling(
const Teuchos::RCP<NOX::Epetra::Scaling>& s)
182bool NOX::Epetra::LinearSystemMPBD::
183computeJacobian(
const NOX::Epetra::Vector& x)
185 bool success = jacInterfacePtr->computeJacobian(
x.getEpetraVector(),
187 block_ops = mp_op->getMPOps();
192bool NOX::Epetra::LinearSystemMPBD::
193createPreconditioner(
const NOX::Epetra::Vector& x,
194 Teuchos::ParameterList& p,
195 bool recomputeGraph)
const
197 EpetraExt::BlockVector mp_x_block(View, *base_map,
x.getEpetraVector());
198 Teuchos::ParameterList& solverParams =
199 p.sublist(
"Deterministic Solver Parameters");
200 bool total_success =
true;
201 if (precStrategy == STANDARD) {
202 for (
int i=0; i<num_mp_blocks; i++) {
203 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
204 if (precs[i] != Teuchos::null)
205 block_solver->setPrecOperatorForSolve(precs[i]);
207 block_solver->createPreconditioner(*(mp_x_block.GetBlock(i)),
208 solverParams, recomputeGraph);
209 precs[i] = block_solver->getGeneratedPrecOperator();
210 total_success = total_success && success;
214 else if (precStrategy ==
MEAN) {
215 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
217 block_solver->createPreconditioner(*(mp_x_block.GetBlock(0)),
218 solverParams, recomputeGraph);
219 total_success = total_success && success;
222 else if (precStrategy == ON_THE_FLY) {
223 if (prec_x == Teuchos::null)
224 prec_x = Teuchos::rcp(
new EpetraExt::BlockVector(mp_x_block));
226 *prec_x = mp_x_block;
229 return total_success;
232bool NOX::Epetra::LinearSystemMPBD::
233destroyPreconditioner()
const
235 return block_solver->destroyPreconditioner();
238bool NOX::Epetra::LinearSystemMPBD::
239recomputePreconditioner(
const NOX::Epetra::Vector& x,
240 Teuchos::ParameterList& p)
const
242 EpetraExt::BlockVector mp_x_block(View, *base_map,
x.getEpetraVector());
243 Teuchos::ParameterList& solverParams =
244 p.sublist(
"Deterministic Solver Parameters");
245 bool total_success =
true;
247 if (precStrategy == STANDARD) {
248 for (
int i=0; i<num_mp_blocks; i++) {
249 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
250 if (precs[i] != Teuchos::null)
251 block_solver->setPrecOperatorForSolve(precs[i]);
253 block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(i)),
255 precs[i] = block_solver->getGeneratedPrecOperator();
256 total_success = total_success && success;
260 else if (precStrategy ==
MEAN) {
261 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
263 block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(0)),
265 total_success = total_success && success;
268 else if (precStrategy == ON_THE_FLY) {
269 if (prec_x == Teuchos::null)
270 prec_x = Teuchos::rcp(
new EpetraExt::BlockVector(mp_x_block));
272 *prec_x = mp_x_block;
275 return total_success;
278NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
279NOX::Epetra::LinearSystemMPBD::
280getPreconditionerPolicy(
bool advanceReuseCounter)
282 return block_solver->getPreconditionerPolicy(advanceReuseCounter);
285bool NOX::Epetra::LinearSystemMPBD::
286isPreconditionerConstructed()
const
288 return block_solver->isPreconditionerConstructed();
291bool NOX::Epetra::LinearSystemMPBD::
292hasPreconditioner()
const
294 return block_solver->hasPreconditioner();
297Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
298getJacobianOperator()
const
303Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
309Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
310getGeneratedPrecOperator()
const
312 if (precStrategy ==
MEAN)
313 return block_solver->getGeneratedPrecOperator();
315 TEUCHOS_TEST_FOR_EXCEPTION(
316 true, std::logic_error,
317 "Cannot call getGeneratedPrecOperator() unless prec strategy is Mean");
318 return Teuchos::null;
321Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
322getGeneratedPrecOperator()
324 if (precStrategy ==
MEAN)
325 return block_solver->getGeneratedPrecOperator();
327 TEUCHOS_TEST_FOR_EXCEPTION(
328 true, std::logic_error,
329 "Cannot call getGeneratedPrecOperator() unless prec strategy is Mean");
330 return Teuchos::null;
333void NOX::Epetra::LinearSystemMPBD::
334setJacobianOperatorForSolve(
const
335 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
337 Teuchos::RCP<const Stokhos::BlockDiagonalOperator> const_mp_op =
338 Teuchos::rcp_dynamic_cast<const Stokhos::BlockDiagonalOperator>(solveJacOp,
340 mp_op = Teuchos::rcp_const_cast<Stokhos::BlockDiagonalOperator>(const_mp_op);
341 block_ops = mp_op->getMPOps();
344void NOX::Epetra::LinearSystemMPBD::
345setPrecOperatorForSolve(
const Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
347 if (precStrategy ==
MEAN)
348 block_solver->setPrecOperatorForSolve(solvePrecOp);
350 TEUCHOS_TEST_FOR_EXCEPTION(
351 true, std::logic_error,
352 "Cannot call setPrecOperatorForSolve() unless prec strategy is Mean");
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x