Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GaussSeidelPreconditioner.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
45#include "Teuchos_TimeMonitor.hpp"
46#include "EpetraExt_BlockUtility.h"
47#include "Teuchos_Assert.hpp"
48
49Stokhos::GaussSeidelPreconditioner::
50GaussSeidelPreconditioner(
51 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
52 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
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"),
59 sg_comm(sg_comm_),
60 sg_basis(sg_basis_),
61 epetraCijk(epetraCijk_),
62 base_map(base_map_),
63 sg_map(sg_map_),
64 is_stoch_parallel(epetraCijk->isStochasticParallel()),
65 stoch_row_map(epetraCijk->getStochasticRowMap()),
66 det_solver(det_solver_),
67 params(params_),
68 useTranspose(false),
69 sg_op(),
70 sg_poly(),
71 Cijk(epetraCijk->getParallelCijk()),
72 is_parallel(epetraCijk->isStochasticParallel())
73{
74 if (is_parallel) {
75 Teuchos::RCP<const Epetra_BlockMap> stoch_col_map =
76 epetraCijk->getStochasticColMap();
77 sg_col_map =
78 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*base_map,
79 *stoch_col_map,
80 sg_map->Comm()));
81 col_exporter = Teuchos::rcp(new Epetra_Export(*sg_col_map, *sg_map));
82 }
83}
84
85Stokhos::GaussSeidelPreconditioner::
86~GaussSeidelPreconditioner()
87{
88}
89
90void
91Stokhos::GaussSeidelPreconditioner::
92setupPreconditioner(const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
93 const Epetra_Vector& x)
94{
95 sg_op = sg_op_;
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"),
101 false);
102}
103
104int
105Stokhos::GaussSeidelPreconditioner::
106SetUseTranspose(bool UseTranspose)
107{
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);
113
114 return 0;
115}
116
117int
118Stokhos::GaussSeidelPreconditioner::
119Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
120{
121 return sg_op->Apply(Input,Result);
122}
123
124int
125Stokhos::GaussSeidelPreconditioner::
126ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
127{
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);
132
133 // We have to be careful if Input and Result are the same vector.
134 // If this is the case, the only possible solution is to make a copy
135 const Epetra_MultiVector *input = &Input;
136 bool made_copy = false;
137 if (Input.Values() == Result.Values()) {
138 input = new Epetra_MultiVector(Input);
139 made_copy = true;
140 }
141
142 int k_limit = sg_poly->size();
143 int dim = sg_poly->basis()->dimension();
144 if (only_use_linear && sg_poly->size() > dim+1)
145 k_limit = dim + 1;
146 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
147
148 int m = input->NumVectors();
149 if (sg_df_block == Teuchos::null || sg_df_block->NumVectors() != m) {
150 sg_df_block =
151 Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
152 sg_y_block =
153 Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
154 kx = Teuchos::rcp(new Epetra_MultiVector(*base_map, m));
155 if (is_parallel) {
156 sg_df_col = Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map,
157 *sg_col_map, m));
158 sg_df_tmp = Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map,
159 *sg_map, m));
160 }
161 }
162
163 // Extract blocks
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);
167
168 // Compute initial residual norm
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);
174
175 Teuchos::RCP<Epetra_MultiVector> f, df, dx;
176
177 sg_df_block->Update(1.0, sg_f_block, 0.0);
178
179 int iter = 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");
183#endif
184 iter++;
185
186 sg_y_block->Update(1.0, sg_f_block, 0.0);
187 if (is_parallel)
188 sg_df_col->PutScalar(0.0);
189
190 for (Cijk_type::i_iterator i_it=Cijk->i_begin();
191 i_it!=Cijk->i_end(); ++i_it) {
192 int i = index(i_it);
193
194 f = sg_f_block.GetBlock(i);
195 df = sg_df_block->GetBlock(i);
196 dx = sg_dx_block.GetBlock(i);
197
198 dx->PutScalar(0.0);
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);
206 // Solve linear system
207 {
208#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
209 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total deterministic solve Time");
210#endif
211 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
212 }
213 }
214
215 df->Update(1.0, *f, 0.0);
216
217 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
218 k_it != Cijk->k_end(i_it); ++k_it) {
219 int k = index(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) {
224 int j = index(j_it);
225 int j_gid = epetraCijk->GCID(j);
226 double c = value(j_it);
227 if (scale_op)
228 c /= norms[j_gid];
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);
233 }
234 else
235 sg_df_col->GetBlock(j)->Update(-c, *kx, 1.0);
236 }
237 }
238 }
239
240 (*sg_poly)[0].Apply(*dx, *kx);
241 sg_y_block->GetBlock(i)->Update(-1.0, *kx, 1.0);
242
243 } //End of k loop
244
245 // Add in contributions from off-process
246 if (is_parallel) {
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);
250 }
251
252 sg_y_block->Norm2(&norm_df);
253 //std::cout << "norm_df = " << norm_df << std::endl;
254 } //End of iter loop
255
256 if (made_copy)
257 delete input;
258
259 return 0;
260}
261
262double
263Stokhos::GaussSeidelPreconditioner::
264NormInf() const
265{
266 return sg_op->NormInf();
267}
268
269
270const char*
271Stokhos::GaussSeidelPreconditioner::
272Label () const
273{
274 return const_cast<char*>(label.c_str());
275}
276
277bool
278Stokhos::GaussSeidelPreconditioner::
279UseTranspose() const
280{
281 return useTranspose;
282}
283
284bool
285Stokhos::GaussSeidelPreconditioner::
286HasNormInf() const
287{
288 return sg_op->HasNormInf();
289}
290
291const Epetra_Comm &
292Stokhos::GaussSeidelPreconditioner::
293Comm() const
294{
295 return *sg_comm;
296}
297const Epetra_Map&
298Stokhos::GaussSeidelPreconditioner::
299OperatorDomainMap() const
300{
301 return *sg_map;
302}
303
304const Epetra_Map&
305Stokhos::GaussSeidelPreconditioner::
306OperatorRangeMap() const
307{
308 return *sg_map;
309}
expr expr expr dx(i, j)
int NormInf(double *Result) const
int NumVectors() const
Abstract base class for multivariate orthogonal polynomials.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)