49 #include "NOX_Epetra_Interface_Jacobian.H" 50 #include "EpetraExt_BlockVector.h" 51 #include "EpetraExt_BlockUtility.h" 52 #include "Teuchos_ParameterList.hpp" 53 #include "Teuchos_TimeMonitor.hpp" 55 NOX::Epetra::LinearSystemSGGS::
57 Teuchos::ParameterList& printingParams,
58 Teuchos::ParameterList& linearSolverParams,
59 const Teuchos::RCP<NOX::Epetra::LinearSystem>& det_solver_,
60 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
61 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
63 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data,
64 const Teuchos::RCP<Epetra_Operator>& J,
65 const Teuchos::RCP<const Epetra_Map>& base_map_,
66 const Teuchos::RCP<const Epetra_Map>& sg_map_,
67 const Teuchos::RCP<NOX::Epetra::Scaling> s):
68 det_solver(det_solver_),
70 epetraCijk(sg_parallel_data->getEpetraCijk()),
71 is_stoch_parallel(epetraCijk->isStochasticParallel()),
72 stoch_row_map(epetraCijk->getStochasticRowMap()),
73 Cijk(epetraCijk->getParallelCijk()),
74 jacInterfacePtr(iJac),
78 utils(printingParams),
79 is_parallel(epetraCijk->isStochasticParallel())
84 sg_df_block = Teuchos::rcp(
new EpetraExt::BlockVector(*base_map, *sg_map));
85 sg_y_block = Teuchos::rcp(
new EpetraExt::BlockVector(*base_map, *sg_map));
88 only_use_linear = linearSolverParams.get(
"Only Use Linear Terms",
false);
89 k_limit = sg_poly->size();
90 int dim = sg_basis->dimension();
91 if (only_use_linear && sg_poly->size() > dim+1)
95 Teuchos::RCP<const Epetra_BlockMap> stoch_col_map =
96 epetraCijk->getStochasticColMap();
98 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*base_map,
101 col_exporter = Teuchos::rcp(
new Epetra_Export(*sg_col_map, *sg_map));
102 sg_df_col = Teuchos::rcp(
new EpetraExt::BlockVector(*base_map,
104 sg_df_tmp = Teuchos::rcp(
new EpetraExt::BlockVector(*base_map, *sg_map));
108 NOX::Epetra::LinearSystemSGGS::
113 bool NOX::Epetra::LinearSystemSGGS::
114 applyJacobian(
const NOX::Epetra::Vector& input,
115 NOX::Epetra::Vector& result)
const 117 sg_op->SetUseTranspose(
false);
118 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
120 return (status == 0);
123 bool NOX::Epetra::LinearSystemSGGS::
124 applyJacobianTranspose(
const NOX::Epetra::Vector& input,
125 NOX::Epetra::Vector& result)
const 127 sg_op->SetUseTranspose(
true);
128 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
129 sg_op->SetUseTranspose(
false);
131 return (status == 0);
134 bool NOX::Epetra::LinearSystemSGGS::
135 applyJacobianInverse(Teuchos::ParameterList ¶ms,
136 const NOX::Epetra::Vector &input,
137 NOX::Epetra::Vector &result)
139 int max_iter = params.get(
"Max Iterations",100 );
140 double sg_tol = params.get(
"Tolerance", 1e-12);
141 bool scale_op = params.get(
"Scale Operator by Inverse Basis Norms",
true);
144 EpetraExt::BlockVector sg_dx_block(View, *base_map, result.getEpetraVector());
145 EpetraExt::BlockVector sg_f_block(View, *base_map, input.getEpetraVector());
146 sg_dx_block.PutScalar(0.0);
149 double norm_f,norm_df;
150 sg_f_block.Norm2(&norm_f);
151 sg_op->Apply(sg_dx_block, *sg_df_block);
152 sg_df_block->Update(-1.0, sg_f_block, 1.0);
153 sg_df_block->Norm2(&norm_df);
155 Teuchos::RCP<Epetra_Vector>
f, df,
dx;
156 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
165 while (((norm_df/norm_f)>sg_tol) && (iter<max_iter)) {
166 TEUCHOS_FUNC_TIME_MONITOR(
"Total global solve Time");
169 sg_y_block->Update(1.0, sg_f_block, 0.0);
171 sg_df_col->PutScalar(0.0);
173 for (Cijk_type::i_iterator i_it=Cijk->i_begin();
174 i_it!=Cijk->i_end(); ++i_it) {
175 int i = Stokhos::index(i_it);
176 f = sg_f_block.GetBlock(i);
177 df = sg_df_block->GetBlock(i);
178 dx = sg_dx_block.GetBlock(i);
181 Teuchos::ParameterList& det_solver_params =
182 params.sublist(
"Deterministic Solver Parameters");
183 NOX::Epetra::Vector nox_df(df, NOX::Epetra::Vector::CreateView);
184 NOX::Epetra::Vector nox_dx(
dx, NOX::Epetra::Vector::CreateView);
187 TEUCHOS_FUNC_TIME_MONITOR(
"Total deterministic solve Time");
188 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
191 df->Update(1.0, *
f, 0.0);
193 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
194 k_it != Cijk->k_end(i_it); ++k_it) {
195 int k = Stokhos::index(k_it);
196 if (k!=0 && k<k_limit) {
197 (*sg_poly)[k].Apply(*
dx, *kx);
198 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
199 j_it != Cijk->j_end(k_it); ++j_it) {
200 int j = Stokhos::index(j_it);
201 int j_gid = epetraCijk->GCID(
j);
202 double c = Stokhos::value(j_it);
205 bool owned = epetraCijk->myGRID(j_gid);
206 if (!is_parallel || owned) {
207 sg_df_block->GetBlock(
j)->Update(-c, *kx, 1.0);
208 sg_y_block->GetBlock(
j)->Update(-c, *kx, 1.0);
211 sg_df_col->GetBlock(
j)->Update(-c, *kx, 1.0);
216 (*sg_poly)[0].Apply(*
dx, *kx);
217 sg_y_block->GetBlock(i)->Update(-1.0, *kx, 1.0);
223 sg_df_tmp->Export(*sg_df_col, *col_exporter, InsertAdd);
224 sg_df_block->Update(1.0, *sg_df_tmp, 1.0);
225 sg_y_block->Update(1.0, *sg_df_tmp, 1.0);
228 sg_y_block->Norm2(&norm_df);
229 utils.out() <<
"\t Gauss-Seidel relative residual norm at iteration " 230 << iter <<
" is " << norm_df/norm_f << std::endl;
237 bool NOX::Epetra::LinearSystemSGGS::
238 applyRightPreconditioning(
bool useTranspose,
239 Teuchos::ParameterList& params,
240 const NOX::Epetra::Vector& input,
241 NOX::Epetra::Vector& result)
const 246 Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemSGGS::
252 void NOX::Epetra::LinearSystemSGGS::
253 resetScaling(
const Teuchos::RCP<NOX::Epetra::Scaling>& s)
259 bool NOX::Epetra::LinearSystemSGGS::
260 computeJacobian(
const NOX::Epetra::Vector&
x)
262 bool success = jacInterfacePtr->computeJacobian(
x.getEpetraVector(),
264 sg_poly = sg_op->getSGPolynomial();
265 det_solver->setJacobianOperatorForSolve(sg_poly->getCoeffPtr(0));
269 bool NOX::Epetra::LinearSystemSGGS::
270 createPreconditioner(
const NOX::Epetra::Vector&
x,
271 Teuchos::ParameterList& p,
272 bool recomputeGraph)
const 274 EpetraExt::BlockVector sg_x_block(View, *base_map,
x.getEpetraVector());
276 det_solver->createPreconditioner(*(sg_x_block.GetBlock(0)),
277 p.sublist(
"Deterministic Solver Parameters"),
283 bool NOX::Epetra::LinearSystemSGGS::
284 destroyPreconditioner()
const 286 return det_solver->destroyPreconditioner();
289 bool NOX::Epetra::LinearSystemSGGS::
290 recomputePreconditioner(
const NOX::Epetra::Vector&
x,
291 Teuchos::ParameterList& linearSolverParams)
const 293 EpetraExt::BlockVector sg_x_block(View, *base_map,
x.getEpetraVector());
295 det_solver->recomputePreconditioner(
296 *(sg_x_block.GetBlock(0)),
297 linearSolverParams.sublist(
"Deterministic Solver Parameters"));
303 NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
304 NOX::Epetra::LinearSystemSGGS::
305 getPreconditionerPolicy(
bool advanceReuseCounter)
307 return det_solver->getPreconditionerPolicy(advanceReuseCounter);
310 bool NOX::Epetra::LinearSystemSGGS::
311 isPreconditionerConstructed()
const 313 return det_solver->isPreconditionerConstructed();
316 bool NOX::Epetra::LinearSystemSGGS::
317 hasPreconditioner()
const 319 return det_solver->hasPreconditioner();
322 Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
323 getJacobianOperator()
const 328 Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
329 getJacobianOperator()
334 Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
335 getGeneratedPrecOperator()
const 337 return Teuchos::null;
340 Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
341 getGeneratedPrecOperator()
343 return Teuchos::null;
346 void NOX::Epetra::LinearSystemSGGS::
347 setJacobianOperatorForSolve(
const 348 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
350 Teuchos::RCP<const Stokhos::SGOperator> const_sg_op =
356 void NOX::Epetra::LinearSystemSGGS::
357 setPrecOperatorForSolve(
const 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.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)