48 #include "NOX_Epetra_Interface_Jacobian.H" 49 #include "EpetraExt_BlockVector.h" 50 #include "Teuchos_ParameterList.hpp" 51 #include "Teuchos_TimeMonitor.hpp" 54 NOX::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),
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);
110 NOX::Epetra::LinearSystemSGJacobi::
111 ~LinearSystemSGJacobi()
115 bool NOX::Epetra::LinearSystemSGJacobi::
116 applyJacobian(
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);
125 bool NOX::Epetra::LinearSystemSGJacobi::
126 applyJacobianTranspose(
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);
136 bool NOX::Epetra::LinearSystemSGJacobi::
137 applyJacobianInverse(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;
201 bool NOX::Epetra::LinearSystemSGJacobi::
202 applyRightPreconditioning(
bool useTranspose,
203 Teuchos::ParameterList& params,
204 const NOX::Epetra::Vector& input,
205 NOX::Epetra::Vector& result)
const 210 Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemSGJacobi::
216 void NOX::Epetra::LinearSystemSGJacobi::
217 resetScaling(
const Teuchos::RCP<NOX::Epetra::Scaling>& s)
223 bool NOX::Epetra::LinearSystemSGJacobi::
224 computeJacobian(
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);
234 bool NOX::Epetra::LinearSystemSGJacobi::
235 createPreconditioner(
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"),
248 bool NOX::Epetra::LinearSystemSGJacobi::
249 destroyPreconditioner()
const 251 return det_solver->destroyPreconditioner();
254 bool NOX::Epetra::LinearSystemSGJacobi::
255 recomputePreconditioner(
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"));
267 NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
268 NOX::Epetra::LinearSystemSGJacobi::
269 getPreconditionerPolicy(
bool advanceReuseCounter)
271 return det_solver->getPreconditionerPolicy(advanceReuseCounter);
274 bool NOX::Epetra::LinearSystemSGJacobi::
275 isPreconditionerConstructed()
const 277 return det_solver->isPreconditionerConstructed();
280 bool NOX::Epetra::LinearSystemSGJacobi::
281 hasPreconditioner()
const 283 return det_solver->hasPreconditioner();
286 Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
287 getJacobianOperator()
const 292 Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
293 getJacobianOperator()
298 Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
299 getGeneratedPrecOperator()
const 301 return Teuchos::null;
304 Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
305 getGeneratedPrecOperator()
307 return Teuchos::null;
310 void NOX::Epetra::LinearSystemSGJacobi::
311 setJacobianOperatorForSolve(
const 312 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
314 Teuchos::RCP<const Stokhos::SGOperator> const_sg_op =
320 void NOX::Epetra::LinearSystemSGJacobi::
321 setPrecOperatorForSolve(
const 322 Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
An abstract class to represent a generic stochastic Galerkin operator as an Epetra_Operator.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()=0
Get SG polynomial.
Factory for generating stochastic Galerkin preconditioners.