Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_KLReducedMatrixFreeOperator.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 
42 #include "EpetraExt_BlockMultiVector.h"
44 #include "Stokhos_PCEAnasaziKL.hpp"
45 #include "Teuchos_Assert.hpp"
46 #include "Teuchos_TimeMonitor.hpp"
50 #include <sstream>
51 
54  const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
55  const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
56  const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
57  const Teuchos::RCP<const Epetra_Map>& domain_base_map_,
58  const Teuchos::RCP<const Epetra_Map>& range_base_map_,
59  const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
60  const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
61  const Teuchos::RCP<Teuchos::ParameterList>& params_) :
62  label("Stokhos KL Reduced Matrix Free Operator"),
63  sg_comm(sg_comm_),
64  sg_basis(sg_basis_),
65  epetraCijk(epetraCijk_),
66  domain_base_map(domain_base_map_),
67  range_base_map(range_base_map_),
68  domain_sg_map(domain_sg_map_),
69  range_sg_map(range_sg_map_),
70  Cijk(epetraCijk->getParallelCijk()),
71  block_ops(),
72  params(params_),
73  useTranspose(false),
74  expansion_size(sg_basis->size()),
75  num_blocks(0),
76  num_KL(0),
77  num_KL_computed(0),
78  mean(),
79  block_vec_map(),
80  block_vec_poly(),
81  dot_products(),
82  sparse_kl_coeffs(),
83  kl_blocks()
84 {
85  num_KL = params->get("Number of KL Terms", 5);
86  drop_tolerance = params->get("Sparse 3 Tensor Drop Tolerance", 1e-6);
87  do_error_tests = params->get("Do Error Tests", false);
88 }
89 
90 void
93  const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
94 {
95  block_ops = ops;
96  num_blocks = block_ops->size();
97 
98  // Build a vector polynomial out of matrix nonzeros
99  mean = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
100  block_ops->getCoeffPtr(0));
101  block_vec_map =
102  Teuchos::rcp(new Epetra_Map(-1, mean->NumMyNonzeros(), 0,
103  domain_base_map->Comm()));
104  block_vec_poly =
105  Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
106  sg_basis, block_ops->map(), block_vec_map, sg_comm));
107 
108  // Setup KL blocks
109  setup();
110 }
111 
112 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
115 {
116  return block_ops;
117 }
118 
119 Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
122 {
123  return block_ops;
124 }
125 
128 {
129 }
130 
131 int
133 SetUseTranspose(bool UseTheTranspose)
134 {
135  useTranspose = UseTheTranspose;
136  kl_mat_free_op->SetUseTranspose(useTranspose);
137  for (int i=0; i<num_blocks; i++)
138  (*block_ops)[i].SetUseTranspose(useTranspose);
139 
140  return 0;
141 }
142 
143 int
145 Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
146 {
147  return kl_mat_free_op->Apply(Input, Result);
148 }
149 
150 int
152 ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
153 {
154  throw "KLReducedMatrixFreeOperator::ApplyInverse not defined!";
155  return -1;
156 }
157 
158 double
160 NormInf() const
161 {
162  return 1.0;
163 }
164 
165 
166 const char*
168 Label () const
169 {
170  return const_cast<char*>(label.c_str());
171 }
172 
173 bool
176 {
177  return useTranspose;
178 }
179 
180 bool
182 HasNormInf() const
183 {
184  return false;
185 }
186 
187 const Epetra_Comm &
189 Comm() const
190 {
191  return *sg_comm;
192 }
193 const Epetra_Map&
196 {
197  if (useTranspose)
198  return *range_sg_map;
199  return *domain_sg_map;
200 }
201 
202 const Epetra_Map&
205 {
206  if (useTranspose)
207  return *domain_sg_map;
208  return *range_sg_map;
209 }
210 
211 void
214 {
215 #ifdef HAVE_STOKHOS_ANASAZI
216 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
217  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::KLReducedMatrixFreeOperator -- Calculation/setup of KL opeator");
218 #endif
219 
220  mean = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
221  block_ops->getCoeffPtr(0));
222 
223  // Copy matrix coefficients into vectors
224  for (int coeff=0; coeff<num_blocks; coeff++) {
225  Teuchos::RCP<const Epetra_CrsMatrix> block_coeff =
226  Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>
227  (block_ops->getCoeffPtr(coeff));
228  int row = 0;
229  for (int i=0; i<mean->NumMyRows(); i++) {
230  int num_col;
231  mean->NumMyRowEntries(i, num_col);
232  for (int j=0; j<num_col; j++)
233  (*block_vec_poly)[coeff][row++] = (*block_coeff)[i][j];
234  }
235  }
236 
237  int myPID = sg_comm->MyPID();
238 
239  // Compute KL expansion of solution sg_J_vec_poly
240  Stokhos::PCEAnasaziKL pceKL(*block_vec_poly, num_KL);
241  Teuchos::ParameterList anasazi_params = pceKL.getDefaultParams();
242  bool result = pceKL.computeKL(anasazi_params);
243  if (!result && myPID == 0)
244  std::cout << "KL Eigensolver did not converge!" << std::endl;
245  Teuchos::RCP<Epetra_MultiVector> evecs = pceKL.getEigenvectors();
246  Teuchos::Array<double> evals = pceKL.getEigenvalues();
247  //num_KL_computed = evecs->NumVectors();
248  if (myPID == 0)
249  std::cout << "num computed eigenvectors = "
250  << evecs->NumVectors() << std::endl;
251  double kl_tol = params->get("KL Tolerance", 1e-6);
252  num_KL_computed = 0;
253  while (num_KL_computed < evals.size() &&
254  std::sqrt(evals[num_KL_computed]/evals[0]) > kl_tol)
255  num_KL_computed++;
256  if (num_KL_computed == evals.size() && myPID == 0)
257  std::cout << "Can't achieve KL tolerance " << kl_tol
258  << ". Smallest eigenvalue / largest eigenvalue = "
259  << std::sqrt(evals[num_KL_computed-1]/evals[0]) << std::endl;
260  if (myPID == 0)
261  std::cout << "num KL eigenvectors = " << num_KL_computed << std::endl;
262 
263  // Compute dot products of Jacobian blocks and KL eigenvectors
264  dot_products.resize(num_KL_computed);
265  for (int rv=0; rv < num_KL_computed; rv++) {
266  dot_products[rv].resize(num_blocks-1);
267  for (int coeff=1; coeff < num_blocks; coeff++) {
268  double dot;
269  (*block_vec_poly)[coeff].Dot(*((*evecs)(rv)), &dot);
270  dot_products[rv][coeff-1] = dot;
271  }
272  }
273 
274  // Compute KL coefficients
275  const Teuchos::Array<double>& norms = sg_basis->norm_squared();
276  sparse_kl_coeffs =
277  Teuchos::rcp(new Stokhos::Sparse3Tensor<int,double>);
278  for (Cijk_type::i_iterator i_it=Cijk->i_begin();
279  i_it!=Cijk->i_end(); ++i_it) {
280  int i = epetraCijk->GRID(index(i_it));
281  sparse_kl_coeffs->sum_term(i, i, 0, norms[i]);
282  }
283  Cijk_type::k_iterator l_begin = ++(Cijk->k_begin());
284  Cijk_type::k_iterator l_end = Cijk->k_end();
285  for (Cijk_type::k_iterator l_it=l_begin; l_it!=l_end; ++l_it) {
286  int l = index(l_it);
287  for (Cijk_type::kj_iterator j_it = Cijk->j_begin(l_it);
288  j_it != Cijk->j_end(l_it); ++j_it) {
289  int j = epetraCijk->GCID(index(j_it));
290  for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
291  i_it != Cijk->i_end(j_it); ++i_it) {
292  int i = epetraCijk->GRID(index(i_it));
293  double c = value(i_it);
294  for (int k=1; k<num_KL_computed+1; k++) {
295  double dp = dot_products[k-1][l-1];
296  double v = dp*c;
297  if (std::abs(v) > drop_tolerance)
298  sparse_kl_coeffs->sum_term(i,j,k,v);
299  }
300  }
301  }
302  }
303  sparse_kl_coeffs->fillComplete();
304 
305  bool save_tensor = params->get("Save Sparse 3 Tensor To File", false);
306  if (save_tensor) {
307  static int idx = 0;
308  std::string basename = params->get("Sparse 3 Tensor Base Filename",
309  "sparse_KL_coeffs");
310  std::stringstream ss;
311  ss << basename << "_" << idx++ << ".mm";
312  sparse3Tensor2MatrixMarket(*sparse_kl_coeffs,
313  *(epetraCijk->getStochasticRowMap()), ss.str());
314  }
315 
316  // Transform eigenvectors back to matrices
317  kl_blocks.resize(num_KL_computed);
318  Teuchos::RCP<Epetra_BlockMap> kl_map =
319  Teuchos::rcp(new Epetra_LocalMap(num_KL_computed+1, 0,
320  sg_comm->TimeDomainComm()));
321  kl_ops =
322  Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
323  sg_basis, kl_map, domain_base_map, range_base_map,
324  range_sg_map, sg_comm));
325  kl_ops->setCoeffPtr(0, mean);
326  for (int rv=0; rv<num_KL_computed; rv++) {
327  if (kl_blocks[rv] == Teuchos::null)
328  kl_blocks[rv] = Teuchos::rcp(new Epetra_CrsMatrix(*mean));
329  int row = 0;
330  for (int i=0; i<mean->NumMyRows(); i++) {
331  int num_col;
332  mean->NumMyRowEntries(i, num_col);
333  for (int j=0; j<num_col; j++)
334  (*kl_blocks[rv])[i][j] = (*evecs)[rv][row++];
335  }
336  kl_ops->setCoeffPtr(rv+1, kl_blocks[rv]);
337  }
338 
339  Teuchos::RCP<Stokhos::EpetraSparse3Tensor> reducedEpetraCijk =
340  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(
341  sg_basis, sparse_kl_coeffs, sg_comm,
342  epetraCijk->getStochasticRowMap(), sparse_kl_coeffs,
343  0, -1));
344  reducedEpetraCijk->transformToLocal();
345 
346  // Create matrix-free op
347  kl_mat_free_op = Teuchos::rcp(new Stokhos::MatrixFreeOperator(
348  sg_comm, sg_basis, reducedEpetraCijk,
349  domain_base_map, range_base_map,
350  domain_sg_map, range_sg_map, params));
351  kl_mat_free_op->setupOperator(kl_ops);
352 
353  // Check accuracy of KL expansion
354  if (do_error_tests) {
355  Teuchos::Array<double> point(sg_basis->dimension());
356  for (int i=0; i<sg_basis->dimension(); i++)
357  point[i] = 0.5;
358  Teuchos::Array<double> basis_vals(sg_basis->size());
359  sg_basis->evaluateBases(point, basis_vals);
360 
361  Epetra_Vector val(*block_vec_map);
362  Epetra_Vector val_kl(*block_vec_map);
363  block_vec_poly->evaluate(basis_vals, val);
364  val_kl.Update(1.0, (*block_vec_poly)[0], 0.0);
365  Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > rvs(num_KL_computed);
366  Teuchos::Array<double> val_rvs(num_KL_computed);
367  for (int rv=0; rv<num_KL_computed; rv++) {
368  rvs[rv].reset(sg_basis);
369  rvs[rv][0] = 0.0;
370  for (int coeff=1; coeff<num_blocks; coeff++)
371  rvs[rv][coeff] = dot_products[rv][coeff-1];
372  val_rvs[rv] = rvs[rv].evaluate(point, basis_vals);
373  val_kl.Update(val_rvs[rv], *((*evecs)(rv)), 1.0);
374  }
375  double nrm;
376  val.NormInf(&nrm);
377  val.Update(-1.0, val_kl, 1.0);
378  double diff;
379  val.NormInf(&diff);
380  if (myPID == 0)
381  std::cout << "Infinity norm of random field difference = " << diff/nrm
382  << std::endl;
383 
384  // Check accuracy of operator
385  Epetra_Vector op_input(*domain_sg_map), op_result(*range_sg_map), op_kl_result(*range_sg_map);
386  op_input.PutScalar(1.0);
387  Stokhos::MatrixFreeOperator op(sg_comm, sg_basis, epetraCijk,
388  domain_base_map, range_base_map,
389  domain_sg_map, range_sg_map, params);
390  op.setupOperator(block_ops);
391  op.Apply(op_input, op_result);
392  this->Apply(op_input, op_kl_result);
393  op_result.NormInf(&nrm);
394  op_result.Update(-1.0, op_kl_result, 1.0);
395  op_result.NormInf(&diff);
396  if (myPID == 0)
397  std::cout << "Infinity norm of operator difference = " << diff/nrm
398  << std::endl;
399  }
400 
401 #else
402  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
403  "Stokhos::KLReducedMatrixFreeOperator is available " <<
404  "only when configured with Anasazi support!")
405 #endif
406 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
An Epetra operator representing the block stochastic Galerkin operator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
expr val()
void sparse3Tensor2MatrixMarket(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm, const std::string &file)
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
Bi-directional iterator for traversing a sparse array.
bool do_error_tests
Whether to do KL error tests (can be expensive)
KLReducedMatrixFreeOperator(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 > &domain_base_map, const Teuchos::RCP< const Epetra_Map > &range_base_map, const Teuchos::RCP< const Epetra_Map > &domain_sg_map, const Teuchos::RCP< const Epetra_Map > &range_sg_map, const Teuchos::RCP< Teuchos::ParameterList > &params)
Constructor.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
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 ...
double drop_tolerance
Tolerance for dropping entries in sparse 3 tensor.
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j
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 ...
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< YD, YP... > >::value, typename Kokkos::Details::InnerProductSpaceTraits< typename Kokkos::View< XD, XP... >::non_const_value_type >::dot_type >::type dot(const Kokkos::View< XD, XP... > &x, const Kokkos::View< YD, YP... > &y)
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
virtual const char * Label() const
Returns a character string describing the operator.
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 ...
virtual bool UseTranspose() const
Returns the current UseTranspose setting.