Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
95 void
97 setupPreconditioner(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 
153 int
155 SetUseTranspose(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 
167 int
169 Apply(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 
215 int
217 ApplyInverse(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 
279 double
281 NormInf() const
282 {
283  // I think this is only an upper bound
284  return mean_prec->NormInf() * G_prec->NormInf();
285 }
286 
287 
288 const char*
290 Label() const
291 {
292  return const_cast<char*>(label.c_str());
293 }
294 
295 bool
298 {
299  return useTranspose;
300 }
301 
302 bool
304 HasNormInf() const
305 {
306  return mean_prec->NormInf() && G_prec->NormInf();
307 }
308 
309 const Epetra_Comm &
311 Comm() const
312 {
313  return *sg_comm;
314 }
315 const Epetra_Map&
318 {
319  return *sg_map;
320 }
321 
322 const Epetra_Map&
325 {
326  return *sg_map;
327 }
328 
329 double
331 MatrixTrace(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 
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
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.
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.
Bi-directional iterator for traversing a sparse array.
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>
Teuchos::RCP< Teuchos::ParameterList > params
Preconditioner parameters.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
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.
double MatrixTrace(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B) const
Compute trace of matrix A&#39;B.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this 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 ...
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
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.