Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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"
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
90void
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
112Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
115{
116 return block_ops;
117}
118
119Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
121getSGPolynomial() const
122{
123 return block_ops;
124}
125
130
131int
133SetUseTranspose(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
143int
145Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
146{
147 return kl_mat_free_op->Apply(Input, Result);
148}
149
150int
152ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
153{
154 throw "KLReducedMatrixFreeOperator::ApplyInverse not defined!";
155 return -1;
156}
157
158double
160NormInf() const
161{
162 return 1.0;
163}
164
165
166const char*
168Label () const
169{
170 return const_cast<char*>(label.c_str());
171}
172
173bool
175UseTranspose() const
176{
177 return useTranspose;
178}
179
180bool
182HasNormInf() const
183{
184 return false;
185}
186
187const Epetra_Comm &
189Comm() const
190{
191 return *sg_comm;
192}
193const Epetra_Map&
195OperatorDomainMap() const
196{
197 if (useTranspose)
198 return *range_sg_map;
199 return *domain_sg_map;
200}
201
202const Epetra_Map&
204OperatorRangeMap() const
205{
206 if (useTranspose)
207 return *domain_sg_map;
208 return *range_sg_map;
209}
210
211void
213setup()
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}
expr val()
UnitTestSetup< Kokkos::Cuda > setup
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
double drop_tolerance
Tolerance for dropping entries in sparse 3 tensor.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
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 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 ...
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
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.
bool do_error_tests
Whether to do KL error tests (can be expensive)
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
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 const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual const char * Label() const
Returns a character string describing the operator.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
An Epetra operator representing the block stochastic Galerkin operator.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
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 ...
Abstract base class for multivariate orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
ikj_sparse_array::const_iterator i_iterator
Iterator for looping over i 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.
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)