Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_KroneckerProductPreconditioner.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 "Epetra_LocalMap.h"
45#include "EpetraExt_BlockMultiVector.h"
46#include "Teuchos_Assert.hpp"
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>& mean_prec_factory_,
56 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& G_prec_factory_,
57 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
58 label("Stokhos Kronecker Product 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 mean_prec_factory(mean_prec_factory_),
65 G_prec_factory(G_prec_factory_),
66 params(params_),
67 mean_prec(),
68 useTranspose(false),
69 sg_op(),
70 sg_poly(),
71 Cijk(epetraCijk->getParallelCijk()),
72 scale_op(true),
73 only_use_linear(false)
74{
75 scale_op = params->get("Scale Operator by Inverse Basis Norms", true);
76 only_use_linear = params_->get("Only Use Linear Terms", false);
77 // Build new parallel Cijk if we are only using the linear terms, Cijk
78 // is distributed over proc's, and Cijk includes more than just the linear
79 // terms (so we have the right column map; otherwise we will be importing
80 // much more than necessary)
81 if (only_use_linear && epetraCijk->isStochasticParallel()) {
82 int dim = sg_basis->dimension();
83 if (epetraCijk->getKEnd() > dim+1)
84 epetraCijk =
85 Teuchos::rcp(new EpetraSparse3Tensor(*epetraCijk, 1, dim+1));
86
87 }
88}
89
92{
93}
94
95void
97setupPreconditioner(const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
98 const Epetra_Vector& x)
99{
100 sg_op = sg_op_;
101 sg_poly = sg_op->getSGPolynomial();
102 mean_prec = mean_prec_factory->compute(sg_poly->getCoeffPtr(0));
103 label = std::string("Stokhos Kronecker Product Preconditioner:\n") +
104 std::string(" ***** ") +
105 std::string(mean_prec->Label());
106
107 // Build graph of G matrix
108 Teuchos::RCP<const Epetra_CrsGraph> graph = epetraCijk->getStochasticGraph();
109
110 // Construct G matrix: G_{ij} = \sum tr(A'B)/tr(A'A)*<psi_alpha,psi_i,psi_j>.
111 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
112 G = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
113 Teuchos::RCP<Epetra_CrsMatrix> A0 =
114 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_poly->getCoeffPtr(0), true);
115 double traceAB0 = MatrixTrace(*A0, *A0);
116 Cijk_type::k_iterator k_begin = Cijk->k_begin();
117 Cijk_type::k_iterator k_end = Cijk->k_end();
118 if (only_use_linear) {
119 int dim = sg_basis->dimension();
120 k_end = Cijk->find_k(dim+1);
121 }
122 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
123 int k = index(k_it);
124 Teuchos::RCP<Epetra_CrsMatrix> A_k =
125 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_poly->getCoeffPtr(k),
126 true);
127 double traceAB = MatrixTrace(*A_k, *A0);
128 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
129 j_it != Cijk->j_end(k_it); ++j_it) {
130 int j = epetraCijk->GCID(index(j_it));
131 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
132 i_it != Cijk->i_end(j_it); ++i_it) {
133 int i = epetraCijk->GRID(index(i_it));
134 double c = value(i_it)*traceAB/traceAB0;
135 if (scale_op)
136 c /= norms[i];
137 G->SumIntoGlobalValues(i, 1, &c, &j);
138 }
139 }
140 }
141 G->FillComplete();
142
143 // Build G preconditioner
144 G_prec = G_prec_factory->compute(G);
145
146 label = std::string("Stokhos Kronecker Product Preconditioner:\n") +
147 std::string(" ***** ") +
148 std::string(mean_prec->Label()) + std::string("\n") +
149 std::string(" ***** ") +
150 std::string(G_prec->Label());
151}
152
153int
155SetUseTranspose(bool UseTranspose)
156{
157 useTranspose = UseTranspose;
158 mean_prec->SetUseTranspose(useTranspose);
159 TEUCHOS_TEST_FOR_EXCEPTION(
160 UseTranspose == true, std::logic_error,
161 "Stokhos::KroneckerProductPreconditioner::SetUseTranspose(): " <<
162 "Preconditioner does not support transpose!" << std::endl);
163
164 return 0;
165}
166
167int
169Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
170{
171 EpetraExt::BlockMultiVector sg_input(View, *base_map, Input);
172 EpetraExt::BlockMultiVector sg_result(View, *base_map, Result);
173
174 const Epetra_Map& G_map = G->RowMap();
175 int NumMyElements = G_map.NumMyElements();
176 int vecLen = sg_input.GetBlock(0)->MyLength(); // Global length of vector.
177 int m = sg_input.NumVectors();
178
179 if (result_MVT == Teuchos::null || result_MVT->NumVectors() != vecLen*m) {
180 result_MVT = Teuchos::rcp(new Epetra_MultiVector(G_map, vecLen*m));
181 }
182
183 // Preconditioner is P = (G x I)(I x A_0)
184
185 // Apply I x A_0
186 for (int i=0; i<NumMyElements; i++) {
187 mean_prec->Apply(*(sg_input.GetBlock(i)), *(sg_result.GetBlock(i)));
188 }
189
190 Teuchos::RCP<Epetra_MultiVector> x;
191 for (int irow=0 ; irow<NumMyElements; irow++) {
192 x = sg_result.GetBlock(irow);
193 for (int vcol=0; vcol<m; vcol++) {
194 for (int icol=0; icol<vecLen; icol++) {
195 (*result_MVT)[m*vcol+icol][irow] = (*x)[vcol][icol];
196 }
197 }
198 }
199
200 // Apply G x I
201 G_prec->Apply(*result_MVT, *result_MVT);
202
203 for (int irow=0; irow<NumMyElements; irow++) {
204 x = sg_result.GetBlock(irow);
205 for (int vcol=0; vcol<m; vcol++) {
206 for (int icol=0; icol<vecLen; icol++) {
207 (*x)[vcol][icol] = (*result_MVT)[m*vcol+icol][irow];
208 }
209 }
210 }
211
212 return 0;
213}
214
215int
217ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
218{
219#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
220 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Kronecker Product Prec Time");
221#endif
222
223 EpetraExt::BlockMultiVector sg_input(View, *base_map, Input);
224 EpetraExt::BlockMultiVector sg_result(View, *base_map, Result);
225
226 const Epetra_Map& G_map = G->RowMap();
227 int NumMyElements = G_map.NumMyElements();
228 int vecLen = sg_input.GetBlock(0)->MyLength(); // Global length of vector.
229 int m = sg_input.NumVectors();
230
231 if (result_MVT == Teuchos::null || result_MVT->NumVectors() != vecLen*m) {
232 result_MVT = Teuchos::rcp(new Epetra_MultiVector(G_map, vecLen*m));
233 }
234
235 // Preconditioner is P^{-1} = (I x A_0^{-1})(G^{-1} x I)
236 // => y = P^{-1}x => Y = A_0^{-1}(G^{-1}X^T)^T where X = multivec(x)
237
238 Teuchos::RCP<Epetra_MultiVector> x;
239 for (int irow=0 ; irow<NumMyElements; irow++) {
240 x = sg_input.GetBlock(irow);
241 for (int vcol=0; vcol<m; vcol++) {
242 for (int icol=0; icol<vecLen; icol++) {
243 (*result_MVT)[m*vcol+icol][irow] = (*x)[vcol][icol];
244 }
245 }
246 }
247
248 // Apply (G^{-1} x I)
249 {
250#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
251 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: G Preconditioner Apply Inverse");
252#endif
253 G_prec->ApplyInverse(*result_MVT, *result_MVT);
254 }
255
256 for (int irow=0; irow<NumMyElements; irow++) {
257 x = sg_result.GetBlock(irow);
258 for (int vcol=0; vcol<m; vcol++) {
259 for (int icol=0; icol<vecLen; icol++) {
260 (*x)[vcol][icol] = (*result_MVT)[m*vcol+icol][irow];
261 }
262 }
263 }
264
265 // Apply (I x A_0^{-1})
266 {
267#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
268 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Mean Preconditioner Apply Inverse");
269#endif
270 for (int i=0; i<NumMyElements; i++) {
271 mean_prec->ApplyInverse(*(sg_result.GetBlock(i)),
272 *(sg_result.GetBlock(i)));
273 }
274 }
275
276 return 0;
277}
278
279double
281NormInf() const
282{
283 // I think this is only an upper bound
284 return mean_prec->NormInf() * G_prec->NormInf();
285}
286
287
288const char*
290Label() const
291{
292 return const_cast<char*>(label.c_str());
293}
294
295bool
297UseTranspose() const
298{
299 return useTranspose;
300}
301
302bool
304HasNormInf() const
305{
306 return mean_prec->NormInf() && G_prec->NormInf();
307}
308
309const Epetra_Comm &
311Comm() const
312{
313 return *sg_comm;
314}
315const Epetra_Map&
317OperatorDomainMap() const
318{
319 return *sg_map;
320}
321
322const Epetra_Map&
324OperatorRangeMap() const
325{
326 return *sg_map;
327}
328
329double
331MatrixTrace(const Epetra_CrsMatrix& A, const Epetra_CrsMatrix& B) const {
332 int n = A.NumMyRows(); // # of rows on the processor.
333 double traceAB = 0.0;
334 for (int i=0; i<n; i++) {
335 int m = A.NumMyEntries(i); // # of non zero entries in row i.
336 for (int j=0; j<m; j++) {
337 traceAB += A[i][j]*B[i][j];
338 }
339 }
340
341 return traceAB;
342}
343
Copy
int NumMyElements() const
int NumMyEntries(int Row) const
int NumMyRows() const
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.
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.
bool only_use_linear
Limit construction of G to linear terms.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
Teuchos::RCP< Teuchos::ParameterList > params
Preconditioner parameters.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
double MatrixTrace(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B) const
Compute trace of matrix A'B.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
KroneckerProductPreconditioner(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 > &mean_prec_factory, const Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > &G_prec_factory, const Teuchos::RCP< Teuchos::ParameterList > &params)
Constructor.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
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.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.