Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
82 void
84 setupPreconditioner(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 
95 int
97 SetUseTranspose(bool UseTheTranspose)
98 {
99  useTranspose = UseTheTranspose;
100  sg_op->SetUseTranspose(UseTheTranspose);
101  mean_prec->SetUseTranspose(UseTheTranspose);
102 
103  return 0;
104 }
105 
106 int
108 Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
109 {
110  return sg_op->Apply(Input, Result);
111 }
112 
113 int
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 
268 double
270 NormInf() const
271 {
272  return sg_op->NormInf();
273 }
274 
275 
276 const char*
278 Label() const
279 {
280  return const_cast<char*>(label.c_str());
281 }
282 
283 bool
286 {
287  return useTranspose;
288 }
289 
290 bool
292 HasNormInf() const
293 {
294  return sg_op->HasNormInf();
295 }
296 
297 const Epetra_Comm &
299 Comm() const
300 {
301  return *sg_comm;
302 }
303 const Epetra_Map&
306 {
307  return *sg_map;
308 }
309 
310 const Epetra_Map&
313 {
314  return *sg_map;
315 }
double * Values() const
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
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.
Bi-directional iterator for traversing a sparse array.
Bi-directional reverse iterator for traversing a sparse array.
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 ...
virtual const char * Label() const
Returns a character string describing the operator.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
bool scale_op
Flag indicating whether operator be scaled with <^2>
int NumVectors() const
bool only_use_linear
Limit Gauss-Seidel loop to linear terms.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
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 ...