45 #include "Teuchos_TimeMonitor.hpp" 46 #include "EpetraExt_BlockUtility.h" 47 #include "Teuchos_Assert.hpp" 49 Stokhos::GaussSeidelPreconditioner::
50 GaussSeidelPreconditioner(
51 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
53 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
54 const Teuchos::RCP<const Epetra_Map>& base_map_,
55 const Teuchos::RCP<const Epetra_Map>& sg_map_,
56 const Teuchos::RCP<NOX::Epetra::LinearSystem>& det_solver_,
57 const Teuchos::RCP<Teuchos::ParameterList>& params_):
58 label(
"Stokhos Gauss-Seidel Preconditioner"),
61 epetraCijk(epetraCijk_),
64 is_stoch_parallel(epetraCijk->isStochasticParallel()),
65 stoch_row_map(epetraCijk->getStochasticRowMap()),
66 det_solver(det_solver_),
71 Cijk(epetraCijk->getParallelCijk()),
72 is_parallel(epetraCijk->isStochasticParallel())
75 Teuchos::RCP<const Epetra_BlockMap> stoch_col_map =
76 epetraCijk->getStochasticColMap();
78 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*base_map,
81 col_exporter = Teuchos::rcp(
new Epetra_Export(*sg_col_map, *sg_map));
85 Stokhos::GaussSeidelPreconditioner::
86 ~GaussSeidelPreconditioner()
91 Stokhos::GaussSeidelPreconditioner::
92 setupPreconditioner(
const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
93 const Epetra_Vector&
x)
96 sg_poly = sg_op->getSGPolynomial();
97 EpetraExt::BlockVector sg_x_block(View, *base_map,
x);
98 det_solver->setJacobianOperatorForSolve(sg_poly->getCoeffPtr(0));
99 det_solver->createPreconditioner(*(sg_x_block.GetBlock(0)),
100 params->sublist(
"Deterministic Solver Parameters"),
105 Stokhos::GaussSeidelPreconditioner::
106 SetUseTranspose(
bool UseTranspose)
108 useTranspose = UseTranspose;
109 TEUCHOS_TEST_FOR_EXCEPTION(
110 UseTranspose ==
true, std::logic_error,
111 "Stokhos::GaussSeidelPreconditioner::SetUseTranspose(): " <<
112 "Preconditioner does not support transpose!" << std::endl);
118 Stokhos::GaussSeidelPreconditioner::
119 Apply(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 121 return sg_op->Apply(Input,Result);
125 Stokhos::GaussSeidelPreconditioner::
126 ApplyInverse(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 128 int max_iter = params->get(
"Max Iterations",100 );
129 double sg_tol = params->get(
"Tolerance", 1e-12);
130 bool scale_op = params->get(
"Scale Operator by Inverse Basis Norms",
true);
131 bool only_use_linear = params->get(
"Only Use Linear Terms",
true);
135 const Epetra_MultiVector *input = &Input;
136 bool made_copy =
false;
137 if (Input.Values() == Result.Values()) {
138 input =
new Epetra_MultiVector(Input);
142 int k_limit = sg_poly->size();
143 int dim = sg_poly->basis()->dimension();
144 if (only_use_linear && sg_poly->size() > dim+1)
146 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
148 int m = input->NumVectors();
149 if (sg_df_block == Teuchos::null || sg_df_block->NumVectors() != m) {
151 Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
153 Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
154 kx = Teuchos::rcp(
new Epetra_MultiVector(*base_map, m));
156 sg_df_col = Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map,
158 sg_df_tmp = Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map,
164 EpetraExt::BlockMultiVector sg_dx_block(View, *base_map, Result);
165 EpetraExt::BlockMultiVector sg_f_block(View, *base_map, *input);
166 sg_dx_block.PutScalar(0.0);
169 double norm_f,norm_df;
170 sg_f_block.Norm2(&norm_f);
171 sg_op->Apply(sg_dx_block, *sg_df_block);
172 sg_df_block->Update(-1.0, sg_f_block, 1.0);
173 sg_df_block->Norm2(&norm_df);
175 Teuchos::RCP<Epetra_MultiVector>
f, df,
dx;
177 sg_df_block->Update(1.0, sg_f_block, 0.0);
180 while (((norm_df/norm_f)>sg_tol) && (iter<max_iter)) {
181 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 182 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total global solve Time");
186 sg_y_block->Update(1.0, sg_f_block, 0.0);
188 sg_df_col->PutScalar(0.0);
190 for (Cijk_type::i_iterator i_it=Cijk->i_begin();
191 i_it!=Cijk->i_end(); ++i_it) {
194 f = sg_f_block.GetBlock(i);
195 df = sg_df_block->GetBlock(i);
196 dx = sg_dx_block.GetBlock(i);
199 Teuchos::ParameterList& det_solver_params =
200 params->sublist(
"Deterministic Solver Parameters");
201 for (
int col=0; col<m; col++) {
202 NOX::Epetra::Vector nox_df(Teuchos::rcp((*df)(col),
false),
203 NOX::Epetra::Vector::CreateView);
204 NOX::Epetra::Vector nox_dx(Teuchos::rcp((*
dx)(col),
false),
205 NOX::Epetra::Vector::CreateView);
208 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 209 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total deterministic solve Time");
211 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
215 df->Update(1.0, *
f, 0.0);
217 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
218 k_it != Cijk->k_end(i_it); ++k_it) {
220 if (k!=0 && k<k_limit) {
221 (*sg_poly)[k].Apply(*
dx, *kx);
222 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
223 j_it != Cijk->j_end(k_it); ++j_it) {
225 int j_gid = epetraCijk->GCID(
j);
226 double c =
value(j_it);
229 bool owned = epetraCijk->myGRID(j_gid);
230 if (!is_parallel || owned) {
231 sg_df_block->GetBlock(
j)->Update(-c, *kx, 1.0);
232 sg_y_block->GetBlock(
j)->Update(-c, *kx, 1.0);
235 sg_df_col->GetBlock(
j)->Update(-c, *kx, 1.0);
240 (*sg_poly)[0].Apply(*
dx, *kx);
241 sg_y_block->GetBlock(i)->Update(-1.0, *kx, 1.0);
247 sg_df_tmp->Export(*sg_df_col, *col_exporter, InsertAdd);
248 sg_df_block->Update(1.0, *sg_df_tmp, 1.0);
249 sg_y_block->Update(1.0, *sg_df_tmp, 1.0);
252 sg_y_block->Norm2(&norm_df);
263 Stokhos::GaussSeidelPreconditioner::
266 return sg_op->NormInf();
271 Stokhos::GaussSeidelPreconditioner::
274 return const_cast<char*
>(label.c_str());
278 Stokhos::GaussSeidelPreconditioner::
285 Stokhos::GaussSeidelPreconditioner::
288 return sg_op->HasNormInf();
292 Stokhos::GaussSeidelPreconditioner::
298 Stokhos::GaussSeidelPreconditioner::
299 OperatorDomainMap()
const 305 Stokhos::GaussSeidelPreconditioner::
306 OperatorRangeMap()
const SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
SparseArrayIterator< index_iterator, value_iterator >::value_reference value(const SparseArrayIterator< index_iterator, value_iterator > &it)