48#include "NOX_Epetra_Interface_Jacobian.H"
49#include "EpetraExt_BlockVector.h"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_TimeMonitor.hpp"
54NOX::Epetra::LinearSystemSGJacobi::
56 Teuchos::ParameterList& printingParams,
57 Teuchos::ParameterList& linearSolverParams_,
58 const Teuchos::RCP<NOX::Epetra::LinearSystem>& det_solver_,
59 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
60 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
62 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data,
63 const Teuchos::RCP<Epetra_Operator>& J,
64 const Teuchos::RCP<const Epetra_Map>& base_map_,
65 const Teuchos::RCP<const Epetra_Map>& sg_map_,
66 const Teuchos::RCP<NOX::Epetra::Scaling> s):
67 det_solver(det_solver_),
68 epetraCijk(sg_parallel_data->getEpetraCijk()),
69 jacInterfacePtr(iJac),
75 sg_op = Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(J,
true);
76 sg_poly = sg_op->getSGPolynomial();
79 Teuchos::rcp(
new EpetraExt::BlockVector(*base_map, *sg_map));
82 Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
83 Teuchos::rcp(&(linearSolverParams_.sublist(
"Jacobi SG Operator")),
false);
84 sgOpParams->set(
"Include Mean",
false);
85 if (!sgOpParams->isParameter(
"Only Use Linear Terms"))
86 sgOpParams->set(
"Only Use Linear Terms",
false);
92 if (sgOpParams->get<
bool>(
"Only Use Linear Terms") &&
93 epetraCijk->isStochasticParallel()) {
94 int dim = sg_basis->dimension();
95 if (epetraCijk->getKEnd() > dim+1)
98 *epetraCijk, 1, dim+1));
103 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
104 sg_parallel_data->getMultiComm();
106 sg_op_factory.build(sg_comm, sg_basis, epetraCijk,
107 base_map, base_map, sg_map, sg_map);
110NOX::Epetra::LinearSystemSGJacobi::
111~LinearSystemSGJacobi()
115bool NOX::Epetra::LinearSystemSGJacobi::
116applyJacobian(
const NOX::Epetra::Vector& input,
117 NOX::Epetra::Vector& result)
const
119 sg_op->SetUseTranspose(
false);
120 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
122 return (status == 0);
125bool NOX::Epetra::LinearSystemSGJacobi::
126applyJacobianTranspose(
const NOX::Epetra::Vector& input,
127 NOX::Epetra::Vector& result)
const
129 sg_op->SetUseTranspose(
true);
130 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
131 sg_op->SetUseTranspose(
false);
133 return (status == 0);
136bool NOX::Epetra::LinearSystemSGJacobi::
137applyJacobianInverse(Teuchos::ParameterList ¶ms,
138 const NOX::Epetra::Vector &input,
139 NOX::Epetra::Vector &result)
141 int max_iter = params.get(
"Max Iterations",100 );
142 double sg_tol = params.get(
"Tolerance", 1e-12);
145 EpetraExt::BlockVector sg_dx_block(View, *base_map, result.getEpetraVector());
146 EpetraExt::BlockVector sg_f_block(View, *base_map, input.getEpetraVector());
147 sg_dx_block.PutScalar(0.0);
150 double norm_f,norm_df;
151 sg_f_block.Norm2(&norm_f);
152 sg_op->Apply(sg_dx_block, *sg_df_block);
153 sg_df_block->Update(-1.0, sg_f_block, 1.0);
154 sg_df_block->Norm2(&norm_df);
156 Teuchos::RCP<Epetra_Vector> df,
dx;
157 Teuchos::ParameterList& det_solver_params =
158 params.sublist(
"Deterministic Solver Parameters");
160 int myBlockRows = epetraCijk->numMyRows();
162 while (((norm_df/norm_f)>sg_tol) && (iter<max_iter)) {
163 TEUCHOS_FUNC_TIME_MONITOR(
"Total global solve Time");
168 sg_df_block->Update(1.0, sg_f_block, 0.0);
170 mat_free_op->Apply(sg_dx_block, *sg_df_block);
171 sg_df_block->Update(1.0, sg_f_block, -1.0);
174 for (
int i=0; i<myBlockRows; i++) {
175 df = sg_df_block->GetBlock(i);
176 dx = sg_dx_block.GetBlock(i);
177 NOX::Epetra::Vector nox_df(df, NOX::Epetra::Vector::CreateView);
178 NOX::Epetra::Vector nox_dx(
dx, NOX::Epetra::Vector::CreateView);
180 (*sg_poly)[0].Apply(*
dx, *kx);
185 TEUCHOS_FUNC_TIME_MONITOR(
"Total deterministic solve Time");
186 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
189 df->Update(-1.0, *kx, 1.0);
192 sg_df_block->Norm2(&norm_df);
193 utils.out() <<
"\t Jacobi relative residual norm at iteration "
194 << iter <<
" is " << norm_df/norm_f << std::endl;
201bool NOX::Epetra::LinearSystemSGJacobi::
202applyRightPreconditioning(
bool useTranspose,
203 Teuchos::ParameterList& params,
204 const NOX::Epetra::Vector& input,
205 NOX::Epetra::Vector& result)
const
210Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemSGJacobi::
216void NOX::Epetra::LinearSystemSGJacobi::
217resetScaling(
const Teuchos::RCP<NOX::Epetra::Scaling>& s)
223bool NOX::Epetra::LinearSystemSGJacobi::
224computeJacobian(
const NOX::Epetra::Vector& x)
226 bool success = jacInterfacePtr->computeJacobian(
x.getEpetraVector(),
228 sg_poly = sg_op->getSGPolynomial();
229 det_solver->setJacobianOperatorForSolve(sg_poly->getCoeffPtr(0));
230 mat_free_op->setupOperator(sg_poly);
234bool NOX::Epetra::LinearSystemSGJacobi::
235createPreconditioner(
const NOX::Epetra::Vector& x,
236 Teuchos::ParameterList& p,
237 bool recomputeGraph)
const
239 EpetraExt::BlockVector sg_x_block(View, *base_map,
x.getEpetraVector());
241 det_solver->createPreconditioner(*(sg_x_block.GetBlock(0)),
242 p.sublist(
"Deterministic Solver Parameters"),
248bool NOX::Epetra::LinearSystemSGJacobi::
249destroyPreconditioner()
const
251 return det_solver->destroyPreconditioner();
254bool NOX::Epetra::LinearSystemSGJacobi::
255recomputePreconditioner(
const NOX::Epetra::Vector& x,
256 Teuchos::ParameterList& linearSolverParams)
const
258 EpetraExt::BlockVector sg_x_block(View, *base_map,
x.getEpetraVector());
260 det_solver->recomputePreconditioner(
261 *(sg_x_block.GetBlock(0)),
262 linearSolverParams.sublist(
"Deterministic Solver Parameters"));
267NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
268NOX::Epetra::LinearSystemSGJacobi::
269getPreconditionerPolicy(
bool advanceReuseCounter)
271 return det_solver->getPreconditionerPolicy(advanceReuseCounter);
274bool NOX::Epetra::LinearSystemSGJacobi::
275isPreconditionerConstructed()
const
277 return det_solver->isPreconditionerConstructed();
280bool NOX::Epetra::LinearSystemSGJacobi::
281hasPreconditioner()
const
283 return det_solver->hasPreconditioner();
286Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
287getJacobianOperator()
const
292Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
298Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
299getGeneratedPrecOperator()
const
301 return Teuchos::null;
304Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
305getGeneratedPrecOperator()
307 return Teuchos::null;
310void NOX::Epetra::LinearSystemSGJacobi::
311setJacobianOperatorForSolve(
const
312 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
314 Teuchos::RCP<const Stokhos::SGOperator> const_sg_op =
315 Teuchos::rcp_dynamic_cast<const Stokhos::SGOperator>(solveJacOp,
true);
316 sg_op = Teuchos::rcp_const_cast<Stokhos::SGOperator>(const_sg_op);
317 sg_poly = sg_op->getSGPolynomial();
320void NOX::Epetra::LinearSystemSGJacobi::
321setPrecOperatorForSolve(
const
322 Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
Abstract base class for multivariate orthogonal polynomials.
Factory for generating stochastic Galerkin preconditioners.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x