Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ApproxSchurComplementPreconditioner.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#include <algorithm>
45#include <iostream>
46#include <iterator>
47
50 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
51 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
52 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
53 const Teuchos::RCP<const Epetra_Map>& base_map_,
54 const Teuchos::RCP<const Epetra_Map>& sg_map_,
55 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& prec_factory_,
56 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
57 label("Stokhos Approximate Schur Complement Preconditioner"),
58 sg_comm(sg_comm_),
59 sg_basis(sg_basis_),
60 epetraCijk(epetraCijk_),
61 base_map(base_map_),
62 sg_map(sg_map_),
63 prec_factory(prec_factory_),
64 mean_prec(),
65 useTranspose(false),
66 sg_op(),
67 sg_poly(),
68 Cijk(epetraCijk->getParallelCijk()),
69 P(sg_basis->order()),
70 block_indices(P+2),
71 upper_block_Cijk(P+1),
72 lower_block_Cijk(P+1),
73 scale_op(true),
74 symmetric(false),
75 only_use_linear(true),
76 rhs_block()
77{
78 // Check if parallel, which we don't support
79 TEUCHOS_TEST_FOR_EXCEPTION(
80 epetraCijk->isStochasticParallel(), std::logic_error,
81 "Stokhos::ApproxSchurComplementPreconditioner does not support " <<
82 "a parallel stochastic distribution.");
83
84 scale_op = params_->get("Scale Operator by Inverse Basis Norms", true);
85 symmetric = params_->get("Symmetric Gauss-Seidel", false);
86 only_use_linear = params_->get("Only Use Linear Terms", true);
87
88 Cijk_type::k_iterator k_begin = Cijk->k_begin();
89 Cijk_type::k_iterator k_end = Cijk->k_end();
91 k_end = Cijk->find_k(sg_basis()->dimension() + 1);
92
94 for (Cijk_type::k_iterator k=k_begin; k!=k_end; ++k) {
95 int nj = Cijk->num_j(k);
96 if (max_num_mat_vec < nj)
97 max_num_mat_vec = nj;
98 }
99
100 // Get indices for each block
101 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > prod_basis =
102 Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<int,double> >(sg_basis, true);
103 int d = prod_basis->dimension();
104 MultiIndex<int> term(d);
105 for (int p=0; p<=P; p++) {
106 term[0] = p;
107 block_indices[p] = prod_basis->index(term);
108 upper_block_Cijk[p] = Teuchos::rcp(new Cijk_type);
109 lower_block_Cijk[p] = Teuchos::rcp(new Cijk_type);
110 }
111 block_indices[P+1] = sg_basis->size();
112
113 // std::cout << "block_indices = [";
114 // std::copy(block_indices.begin(), block_indices.end(),
115 // std::ostream_iterator<int>(std::cout, " "));
116 // std::cout << "]" << std::endl;
117
118 // Build Cijk tensors for each order block
119 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
120 int k = index(k_it);
121 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
122 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
123 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
124 int j = index(j_it);
125 Teuchos::Array<int>::iterator col_it =
126 std::upper_bound(block_indices.begin(), block_indices.end(), j);
127 int p_col = col_it - block_indices.begin() - 1;
128 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
129 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
130 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
131 int i = index(i_it);
132 double c = value(i_it);
133 Teuchos::Array<int>::iterator row_it =
134 std::upper_bound(block_indices.begin(), block_indices.end(), i);
135 int p_row = row_it - block_indices.begin() - 1;
136 //std::cout << "i = " << i << ", p_row = " << p_row << ", j = " << j << ", p_col = " << p_col;
137 if (p_col > p_row) {
138 upper_block_Cijk[p_col]->add_term(i,j,k,c);
139 //std::cout << " upper" << std::endl;
140 }
141 else if (p_col < p_row) {
142 lower_block_Cijk[p_row]->add_term(i,j,k,c);
143 //std::cout << " lower" << std::endl;
144 }
145 }
146 }
147 }
148 for (int p=0; p<=P; p++) {
149 upper_block_Cijk[p]->fillComplete();
150 lower_block_Cijk[p]->fillComplete();
151 }
152}
153
158
159void
161setupPreconditioner(const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
162 const Epetra_Vector& x)
163{
164 sg_op = sg_op_;
165 sg_poly = sg_op->getSGPolynomial();
166 mean_prec = prec_factory->compute(sg_poly->getCoeffPtr(0));
167 label = std::string("Stokhos Approximate Schur Complement Preconditioner:\n")
168 + std::string(" ***** ") + std::string(mean_prec->Label());
169}
170
171int
173SetUseTranspose(bool UseTranspose)
174{
175 useTranspose = UseTranspose;
176 sg_op->SetUseTranspose(UseTranspose);
177 mean_prec->SetUseTranspose(UseTranspose);
178
179 return 0;
180}
181
182int
184Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
185{
186 return sg_op->Apply(Input, Result);
187}
188
189int
191ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
192{
193#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
194 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Approximate Schur Complement Time");
195#endif
196
197 // We have to be careful if Input and Result are the same vector.
198 // If this is the case, the only possible solution is to make a copy
199 const Epetra_MultiVector *input = &Input;
200 bool made_copy = false;
201 if (Input.Values() == Result.Values()) {
202 input = new Epetra_MultiVector(Input);
203 made_copy = true;
204 }
205
206 // Allocate temporary storage
207 int m = input->NumVectors();
208 if (rhs_block == Teuchos::null || rhs_block->NumVectors() != m)
209 rhs_block =
210 Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
211 if (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec)
212 tmp = Teuchos::rcp(new Epetra_MultiVector(*base_map,
213 m*max_num_mat_vec));
214 j_ptr.resize(m*max_num_mat_vec);
215 mj_indices.resize(m*max_num_mat_vec);
216
217 // Extract blocks
218 EpetraExt::BlockMultiVector input_block(View, *base_map, *input);
219 EpetraExt::BlockMultiVector result_block(View, *base_map, Result);
220
221 result_block.PutScalar(0.0);
222
223 // Set right-hand-side to input_block
224 rhs_block->Update(1.0, input_block, 0.0);
225
226 // At level l, linear system has the structure
227 // [ A_{l-1} B_l ][ u_l^{l-1} ] = [ r_l^{l-1} ]
228 // [ C_l D_l ][ u_l^l ] [ r_l^l ]
229
230 for (int l=P; l>=1; l--) {
231 // Compute D_l^{-1} r_l^l
232 divide_diagonal_block(block_indices[l], block_indices[l+1],
233 *rhs_block, result_block);
234
235 // Compute r_l^{l-1} = r_l^{l-1} - B_l D_l^{-1} r_l^l
236 multiply_block(upper_block_Cijk[l], -1.0, result_block, *rhs_block);
237 }
238
239 // Solve A_0 u_0 = r_0
240 divide_diagonal_block(0, 1, *rhs_block, result_block);
241
242 for (int l=1; l<=P; l++) {
243 // Compute r_l^l - C_l*u_l^{l-1}
244 multiply_block(lower_block_Cijk[l], -1.0, result_block, *rhs_block);
245
246 // Compute D_l^{-1} (r_l^l - C_l*u_l^{l-1})
247 divide_diagonal_block(block_indices[l], block_indices[l+1],
248 *rhs_block, result_block);
249 }
250
251 if (made_copy)
252 delete input;
253
254 return 0;
255}
256
257double
259NormInf() const
260{
261 return sg_op->NormInf();
262}
263
264
265const char*
267Label() const
268{
269 return const_cast<char*>(label.c_str());
270}
271
272bool
274UseTranspose() const
275{
276 return useTranspose;
277}
278
279bool
281HasNormInf() const
282{
283 return sg_op->HasNormInf();
284}
285
286const Epetra_Comm &
288Comm() const
289{
290 return *sg_comm;
291}
292const Epetra_Map&
298
299const Epetra_Map&
305
306void
309 const Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> >& cijk,
310 double alpha,
311 const EpetraExt::BlockMultiVector& Input,
312 EpetraExt::BlockMultiVector& Result) const
313{
314 // Input and Result are the whole vector/multi-vector, not just the portion
315 // needed for the particular sub-block
316 int m = Input.NumVectors();
317 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
318 Cijk_type::k_iterator k_begin = cijk->k_begin();
319 Cijk_type::k_iterator k_end = cijk->k_end();
320 if (only_use_linear)
321 k_end = cijk->find_k(sg_basis()->dimension() + 1);
322 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
323 int k = index(k_it);
324 Cijk_type::kj_iterator j_begin = cijk->j_begin(k_it);
325 Cijk_type::kj_iterator j_end = cijk->j_end(k_it);
326 int nj = cijk->num_j(k_it);
327 if (nj > 0) {
328 int l = 0;
329 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
330 int j = index(j_it);
331 for (int mm=0; mm<m; mm++) {
332 j_ptr[l*m+mm] = (*(Input.GetBlock(j)))[mm];
333 mj_indices[l*m+mm] = l*m+mm;
334 }
335 l++;
336 }
337 Epetra_MultiVector input_tmp(View, *base_map, &j_ptr[0], nj*m);
338 Epetra_MultiVector result_tmp(View, *tmp, &mj_indices[0], nj*m);
339 (*sg_poly)[k].Apply(input_tmp, result_tmp);
340 l = 0;
341 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
342 Cijk_type::kji_iterator i_begin = cijk->i_begin(j_it);
343 Cijk_type::kji_iterator i_end = cijk->i_end(j_it);
344 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
345 int i = index(i_it);
346 double c = value(i_it);
347 if (scale_op)
348 c /= norms[i];
349 for (int mm=0; mm<m; mm++)
350 (*Result.GetBlock(i))(mm)->Update(alpha*c, *result_tmp(l*m+mm), 1.0);
351 }
352 l++;
353 }
354 }
355 }
356}
357
358void
360divide_diagonal_block(int row_begin, int row_end,
361 const EpetraExt::BlockMultiVector& Input,
362 EpetraExt::BlockMultiVector& Result) const
363{
364 for (int i=row_begin; i<row_end; i++) {
365#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
366 TEUCHOS_FUNC_TIME_MONITOR(
367 "Stokhos: ASC Deterministic Preconditioner Time");
368#endif
369 mean_prec->ApplyInverse(*(Input.GetBlock(i)), *(Result.GetBlock(i)));
370 }
371}
int NumVectors() const
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 bool UseTranspose() const
Returns the current UseTranspose setting.
void multiply_block(const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &cijk, double alpha, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
virtual const char * Label() const
Returns a character string describing the operator.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
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.
Teuchos::RCP< const Cijk_type > Cijk
Pointer to triple product.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
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 ...
void divide_diagonal_block(int row_begin, int row_end, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
ApproxSchurComplementPreconditioner(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 bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Teuchos::Array< Teuchos::RCP< Cijk_type > > upper_block_Cijk
Triple product tensor for each sub-block.
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.
Abstract base class for multivariate orthogonal polynomials.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
ji_sparse_array::const_iterator kj_iterator
Iterator for looping over j entries given k.