42#include "EpetraExt_BlockMultiVector.h"
45#include "Teuchos_Assert.hpp"
46#include "Teuchos_TimeMonitor.hpp"
54 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
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"),
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()),
74 expansion_size(sg_basis->size()),
93 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
96 num_blocks = block_ops->size();
99 mean = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
100 block_ops->getCoeffPtr(0));
102 Teuchos::rcp(
new Epetra_Map(-1, mean->NumMyNonzeros(), 0,
103 domain_base_map->Comm()));
106 sg_basis, block_ops->map(), block_vec_map, sg_comm));
112Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
119Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
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);
147 return kl_mat_free_op->Apply(Input, Result);
154 throw "KLReducedMatrixFreeOperator::ApplyInverse not defined!";
170 return const_cast<char*
>(label.c_str());
198 return *range_sg_map;
199 return *domain_sg_map;
207 return *domain_sg_map;
208 return *range_sg_map;
215#ifdef HAVE_STOKHOS_ANASAZI
216#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
217 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::KLReducedMatrixFreeOperator -- Calculation/setup of KL opeator");
220 mean = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
221 block_ops->getCoeffPtr(0));
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));
229 for (
int i=0; i<mean->NumMyRows(); i++) {
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];
237 int myPID = sg_comm->MyPID();
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();
249 std::cout <<
"num computed eigenvectors = "
250 << evecs->NumVectors() << std::endl;
251 double kl_tol = params->get(
"KL Tolerance", 1e-6);
253 while (num_KL_computed < evals.size() &&
254 std::sqrt(evals[num_KL_computed]/evals[0]) > kl_tol)
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;
261 std::cout <<
"num KL eigenvectors = " << num_KL_computed << std::endl;
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++) {
269 (*block_vec_poly)[coeff].Dot(*((*evecs)(rv)), &dot);
270 dot_products[rv][coeff-1] = dot;
275 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
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]);
288 j_it != Cijk->j_end(l_it); ++j_it) {
289 int j = epetraCijk->GCID(index(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];
297 if (std::abs(v) > drop_tolerance)
298 sparse_kl_coeffs->sum_term(i,
j,k,v);
303 sparse_kl_coeffs->fillComplete();
305 bool save_tensor = params->get(
"Save Sparse 3 Tensor To File",
false);
308 std::string basename = params->get(
"Sparse 3 Tensor Base Filename",
310 std::stringstream ss;
311 ss << basename <<
"_" << idx++ <<
".mm";
313 *(epetraCijk->getStochasticRowMap()), ss.str());
317 kl_blocks.resize(num_KL_computed);
318 Teuchos::RCP<Epetra_BlockMap> kl_map =
320 sg_comm->TimeDomainComm()));
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)
330 for (
int i=0; i<mean->NumMyRows(); i++) {
332 mean->NumMyRowEntries(i, num_col);
333 for (
int j=0;
j<num_col;
j++)
334 (*kl_blocks[rv])[i][
j] = (*evecs)[rv][row++];
336 kl_ops->setCoeffPtr(rv+1, kl_blocks[rv]);
339 Teuchos::RCP<Stokhos::EpetraSparse3Tensor> reducedEpetraCijk =
341 sg_basis, sparse_kl_coeffs, sg_comm,
342 epetraCijk->getStochasticRowMap(), sparse_kl_coeffs,
344 reducedEpetraCijk->transformToLocal();
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);
354 if (do_error_tests) {
355 Teuchos::Array<double> point(sg_basis->dimension());
356 for (
int i=0; i<sg_basis->dimension(); i++)
358 Teuchos::Array<double> basis_vals(sg_basis->size());
359 sg_basis->evaluateBases(point, basis_vals);
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);
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);
377 val.Update(-1.0, val_kl, 1.0);
381 std::cout <<
"Infinity norm of random field difference = " << diff/nrm
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);
388 domain_base_map, range_base_map,
389 domain_sg_map, range_sg_map, params);
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);
397 std::cout <<
"Infinity norm of operator difference = " << diff/nrm
402 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
403 "Stokhos::KLReducedMatrixFreeOperator is available " <<
404 "only when configured with Anasazi support!")
UnitTestSetup< Kokkos::Cuda > setup
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
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 ~KLReducedMatrixFreeOperator()
Destructor.
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 ...
int num_KL
Number of KL terms.
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.
void setup()
Setup KL blocks.
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 > ¶ms)
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.
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)
Bi-directional iterator for traversing a sparse array.