Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ApproxGaussSeidelPreconditioner.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
43#include "Teuchos_TimeMonitor.hpp"
44
47 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
48 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
49 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
50 const Teuchos::RCP<const Epetra_Map>& base_map_,
51 const Teuchos::RCP<const Epetra_Map>& sg_map_,
52 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& prec_factory_,
53 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
54 label("Stokhos Approximate Gauss-Seidel Preconditioner"),
55 sg_comm(sg_comm_),
56 sg_basis(sg_basis_),
57 epetraCijk(epetraCijk_),
58 base_map(base_map_),
59 sg_map(sg_map_),
60 prec_factory(prec_factory_),
61 mean_prec(),
62 useTranspose(false),
63 sg_op(),
64 sg_poly(),
65 Cijk(epetraCijk->getParallelCijk()),
66 scale_op(true),
67 symmetric(false),
68 only_use_linear(true),
69 mat_vec_tmp(),
70 rhs_block()
71{
72 scale_op = params_->get("Scale Operator by Inverse Basis Norms", true);
73 symmetric = params_->get("Symmetric Gauss-Seidel", false);
74 only_use_linear = params_->get("Only Use Linear Terms", true);
75}
76
79{
80}
81
82void
84setupPreconditioner(const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
85 const Epetra_Vector& x)
86{
87 sg_op = sg_op_;
88 sg_poly = sg_op->getSGPolynomial();
89 mean_prec = prec_factory->compute(sg_poly->getCoeffPtr(0));
90 label = std::string("Stokhos Approximate Gauss-Seidel Preconditioner:\n") +
91 std::string(" ***** ") +
92 std::string(mean_prec->Label());
93}
94
95int
97SetUseTranspose(bool UseTheTranspose)
98{
99 useTranspose = UseTheTranspose;
100 sg_op->SetUseTranspose(UseTheTranspose);
101 mean_prec->SetUseTranspose(UseTheTranspose);
102
103 return 0;
104}
105
106int
108Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
109{
110 return sg_op->Apply(Input, Result);
111}
112
113int
115ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
116{
117#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
118 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Approximate Gauss-Seidel Time");
119#endif
120
121 // We have to be careful if Input and Result are the same vector.
122 // If this is the case, the only possible solution is to make a copy
123 const Epetra_MultiVector *input = &Input;
124 bool made_copy = false;
125 if (Input.Values() == Result.Values()) {
126 input = new Epetra_MultiVector(Input);
127 made_copy = true;
128 }
129
130 int m = input->NumVectors();
131 if (mat_vec_tmp == Teuchos::null || mat_vec_tmp->NumVectors() != m)
132 mat_vec_tmp = Teuchos::rcp(new Epetra_MultiVector(*base_map, m));
133 if (rhs_block == Teuchos::null || rhs_block->NumVectors() != m)
134 rhs_block =
135 Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
136
137 // Extract blocks
138 EpetraExt::BlockMultiVector input_block(View, *base_map, *input);
139 EpetraExt::BlockMultiVector result_block(View, *base_map, Result);
140
141 result_block.PutScalar(0.0);
142
143 int k_limit = sg_poly->size();
144 if (only_use_linear)
145 k_limit = sg_poly->basis()->dimension() + 1;
146 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
147
148 rhs_block->Update(1.0, input_block, 0.0);
149
150 for (Cijk_type::i_iterator i_it=Cijk->i_begin();
151 i_it!=Cijk->i_end(); ++i_it) {
152 int i = index(i_it);
153
154 Teuchos::RCP<Epetra_MultiVector> res_i = result_block.GetBlock(i);
155 {
156 // Apply deterministic preconditioner
157#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
158 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total AGS Deterministic Preconditioner Time");
159#endif
160 mean_prec->ApplyInverse(*(rhs_block->GetBlock(i)), *res_i);
161 }
162
163 int i_gid = epetraCijk->GRID(i);
164 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
165 k_it != Cijk->k_end(i_it); ++k_it) {
166 int k = index(k_it);
167 if (k!=0 && k<k_limit) {
168 bool do_mat_vec = false;
169 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
170 j_it != Cijk->j_end(k_it); ++j_it) {
171 int j = index(j_it);
172 int j_gid = epetraCijk->GCID(j);
173 if (j_gid > i_gid) {
174 bool on_proc = epetraCijk->myGRID(j_gid);
175 if (on_proc) {
176 do_mat_vec = true;
177 break;
178 }
179 }
180 }
181 if (do_mat_vec) {
182 (*sg_poly)[k].Apply(*res_i, *mat_vec_tmp);
183 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
184 j_it != Cijk->j_end(k_it); ++j_it) {
185 int j = index(j_it);
186 int j_gid = epetraCijk->GCID(j);
187 double c = value(j_it);
188 if (scale_op) {
189 if (useTranspose)
190 c /= norms[i_gid];
191 else
192 c /= norms[j_gid];
193 }
194 if (j_gid > i_gid) {
195 bool on_proc = epetraCijk->myGRID(j_gid);
196 if (on_proc) {
197 rhs_block->GetBlock(j)->Update(-c, *mat_vec_tmp, 1.0);
198 }
199 }
200 }
201 }
202 }
203 }
204 }
205
206 // For symmetric Gauss-Seidel
207 if (symmetric) {
208
209 for (Cijk_type::i_reverse_iterator i_it= Cijk->i_rbegin();
210 i_it!=Cijk->i_rend(); ++i_it) {
211 int i = index(i_it);
212
213 Teuchos::RCP<Epetra_MultiVector> res_i = result_block.GetBlock(i);
214 {
215 // Apply deterministic preconditioner
216#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
217 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total AGS Deterministic Preconditioner Time");
218#endif
219 mean_prec->ApplyInverse(*(rhs_block->GetBlock(i)), *res_i);
220 }
221
222 int i_gid = epetraCijk->GRID(i);
223 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
224 k_it != Cijk->k_end(i_it); ++k_it) {
225 int k = index(k_it);
226 if (k!=0 && k<k_limit) {
227 bool do_mat_vec = false;
228 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
229 j_it != Cijk->j_end(k_it); ++j_it) {
230 int j = index(j_it);
231 int j_gid = epetraCijk->GCID(j);
232 if (j_gid < i_gid) {
233 bool on_proc = epetraCijk->myGRID(j_gid);
234 if (on_proc) {
235 do_mat_vec = true;
236 break;
237 }
238 }
239 }
240 if (do_mat_vec) {
241 (*sg_poly)[k].Apply(*res_i, *mat_vec_tmp);
242 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
243 j_it != Cijk->j_end(k_it); ++j_it) {
244 int j = index(j_it);
245 int j_gid = epetraCijk->GCID(j);
246 double c = value(j_it);
247 if (scale_op)
248 c /= norms[j_gid];
249 if (j_gid < i_gid) {
250 bool on_proc = epetraCijk->myGRID(j_gid);
251 if (on_proc) {
252 rhs_block->GetBlock(j)->Update(-c, *mat_vec_tmp, 1.0);
253 }
254 }
255 }
256 }
257 }
258 }
259 }
260 }
261
262 if (made_copy)
263 delete input;
264
265 return 0;
266}
267
268double
270NormInf() const
271{
272 return sg_op->NormInf();
273}
274
275
276const char*
278Label() const
279{
280 return const_cast<char*>(label.c_str());
281}
282
283bool
285UseTranspose() const
286{
287 return useTranspose;
288}
289
290bool
292HasNormInf() const
293{
294 return sg_op->HasNormInf();
295}
296
297const Epetra_Comm &
299Comm() const
300{
301 return *sg_comm;
302}
303const Epetra_Map&
305OperatorDomainMap() const
306{
307 return *sg_map;
308}
309
310const Epetra_Map&
312OperatorRangeMap() const
313{
314 return *sg_map;
315}
int NumVectors() const
double * Values() const
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
bool only_use_linear
Limit Gauss-Seidel loop to linear terms.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of the inverse of the operator applied to a Epetra_MultiVector Input in Result as ...
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
ApproxGaussSeidelPreconditioner(const Teuchos::RCP< const EpetraExt::MultiComm > &sg_comm, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > &epetraCijk, const Teuchos::RCP< const Epetra_Map > &base_map, const Teuchos::RCP< const Epetra_Map > &sg_map, const Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > &prec_factory, const Teuchos::RCP< Teuchos::ParameterList > &params)
Constructor.
virtual const char * Label() const
Returns a character string describing the operator.
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector Input in Result as described ...
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.
Bi-directional reverse iterator for traversing a sparse array.