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();
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));
112 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
119 Teuchos::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");
221 block_ops->getCoeffPtr(0));
224 for (
int coeff=0; coeff<num_blocks; coeff++) {
225 Teuchos::RCP<const Epetra_CrsMatrix> block_coeff =
227 (block_ops->getCoeffPtr(coeff));
229 for (
int i=0; i<mean->NumMyRows(); i++) {
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];
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!")
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.
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.
void setup()
Setup KL blocks.
int NumMyRowEntries(int MyRow, int &NumEntries) const
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.
int num_KL
Number of KL terms.
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 > ¶ms)
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.
UnitTestSetup< int, double > setup
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)
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 Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
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 ~KLReducedMatrixFreeOperator()
Destructor.
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.