43#include "EpetraExt_BlockMultiVector.h"
44#include "EpetraExt_BlockUtility.h"
45#include "Teuchos_TimeMonitor.hpp"
49 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
51 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
52 const Teuchos::RCP<const Epetra_Map>& domain_base_map_,
53 const Teuchos::RCP<const Epetra_Map>& range_base_map_,
54 const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
55 const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
56 const Teuchos::RCP<Teuchos::ParameterList>& params) :
57 label(
"Stokhos Matrix Free Operator"),
60 epetraCijk(epetraCijk_),
61 domain_base_map(domain_base_map_),
62 range_base_map(range_base_map_),
63 domain_sg_map(domain_sg_map_),
64 range_sg_map(range_sg_map_),
65 is_stoch_parallel(epetraCijk->isStochasticParallel()),
67 global_col_map_trans(),
68 stoch_col_map(epetraCijk->getStochasticColMap()),
71 Cijk(epetraCijk->getParallelCijk()),
75 only_use_linear(false),
77 use_block_apply(
true),
78 expansion_size(sg_basis->size()),
86 k_begin(Cijk->k_begin()),
89 scale_op = params->get(
"Scale Operator by Inverse Basis Norms",
true);
103 int nj =
Cijk->num_j(k);
118 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*
domain_base_map,
122 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*
range_base_map,
157 i_it != Cijk->i_end(j_it); ++i_it) {
167 return 2.0 *
static_cast<double>(n_apply) *
static_cast<double>(nnz) +
168 static_cast<double>(n_add) *
static_cast<double>(nrow);
174 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
177 num_blocks = block_ops->size();
178 if (num_blocks < Cijk->num_k())
179 k_end = Cijk->find_k(num_blocks);
182Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
189Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
200 useTranspose = UseTheTranspose;
201 for (
int i=0; i<num_blocks; i++)
202 (*block_ops)[i].SetUseTranspose(useTranspose);
211#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
212 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: SG Operator Apply()");
225 bool made_copy =
false;
226 if (Input.Values() == Result.Values() && !is_stoch_parallel) {
232 Result.PutScalar(0.0);
234 const Epetra_Map* input_base_map = domain_base_map.get();
235 const Epetra_Map* result_base_map = range_base_map.get();
236 if (useTranspose ==
true) {
237 input_base_map = range_base_map.get();
238 result_base_map = domain_base_map.get();
242 int m = Input.NumVectors();
243 if (useTranspose ==
false &&
244 (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec))
247 else if (useTranspose ==
true &&
248 (tmp_trans == Teuchos::null ||
249 tmp_trans->NumVectors() != m*max_num_mat_vec))
253 if (useTranspose ==
false)
254 tmp_result = tmp.get();
256 tmp_result = tmp_trans.get();
260 if (!is_stoch_parallel)
263 if (useTranspose ==
false) {
264 if (input_col == Teuchos::null || input_col->NumVectors() != m)
266 input_col->Import(*input, *col_importer,
Insert);
267 tmp_col = input_col.get();
270 if (input_col_trans == Teuchos::null ||
271 input_col_trans->NumVectors() != m)
274 input_col_trans->Import(*input, *col_importer_trans,
Insert);
275 tmp_col = input_col_trans.get();
280 EpetraExt::BlockMultiVector sg_input(
View, *input_base_map, *tmp_col);
281 EpetraExt::BlockMultiVector sg_result(
View, *result_base_map, Result);
282 for (
int i=0; i<input_block.size(); i++)
283 input_block[i] = sg_input.GetBlock(i);
284 for (
int i=0; i<result_block.size(); i++)
285 result_block[i] = sg_result.GetBlock(i);
294 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
299 int nj = Cijk->num_j(k_it);
301 Teuchos::Array<double*> j_ptr(nj*m);
302 Teuchos::Array<int> mj_indices(nj*m);
306 for (
int mm=0; mm<m; mm++) {
307 j_ptr[l*m+mm] = (*input_block[
j])[mm];
308 mj_indices[l*m+mm] = l*m+mm;
314 if (use_block_apply) {
315 (*block_ops)[k].Apply(input_tmp, result_tmp);
318 for (
int jj=0; jj<nj*m; jj++)
319 (*block_ops)[k].Apply(*(input_tmp(jj)), *(result_tmp(jj)));
325 i_it != Cijk->i_end(j_it); ++i_it) {
327 double c = value(i_it);
331 i_gid = epetraCijk->GCID(
j);
333 i_gid = epetraCijk->GRID(i);
336 for (
int mm=0; mm<m; mm++)
337 (*result_block[i])(mm)->
Update(c, *result_tmp(l*m+mm), 1.0);
345 for (
int i=0; i<input_block.size(); i++)
346 input_block[i] = Teuchos::null;
347 for (
int i=0; i<result_block.size(); i++)
348 result_block[i] = Teuchos::null;
360 throw "MatrixFreeOperator::ApplyInverse not defined!";
374 return const_cast<char*
>(label.c_str());
398 return *range_sg_map;
399 return *domain_sg_map;
406 return *domain_sg_map;
407 return *range_sg_map;
int NumGlobalNonzeros() const
int NumGlobalRows() const
virtual const char * Label() const
Returns a character string describing the operator.
Teuchos::RCP< const Epetra_Map > range_base_map
Stores range base map.
Teuchos::RCP< Epetra_Import > col_importer
Importer from domain map to column map.
Teuchos::RCP< const Epetra_Map > range_sg_map
Stores range SG map.
Cijk_type::k_iterator k_end
Ending k iterator.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual ~MatrixFreeOperator()
Destructor.
MatrixFreeOperator(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.
Teuchos::RCP< const Epetra_BlockMap > stoch_col_map
Stores stochastic part of column map.
Teuchos::RCP< const Epetra_Map > domain_sg_map
Stores domain SG map.
Teuchos::Array< Teuchos::RCP< Epetra_MultiVector > > result_block
MultiVectors for each block for Apply() result.
bool include_mean
Flag indicating whether to include mean term.
Teuchos::RCP< const Cijk_type > Cijk
Stores triple product tensor.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Teuchos::RCP< Epetra_Map > global_col_map_trans
Stores operator column SG map for transpose.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
Teuchos::Array< Teuchos::RCP< const Epetra_MultiVector > > input_block
MultiVectors for each block for Apply() input.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
bool use_block_apply
Flag indicating whether to use block Epetra_MultiVector apply.
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 ...
bool is_stoch_parallel
Whether we have parallelism over stochastic blocks.
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 ...
Cijk_type::k_iterator k_begin
Starting k iterator.
bool only_use_linear
Flag indicating whether to only use linear terms.
Teuchos::RCP< Epetra_Map > global_col_map
Stores operator column SG map.
int expansion_size
Number of terms in expansion.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
Teuchos::RCP< const Epetra_Map > domain_base_map
Stores domain base map.
double countApplyFlops() const
Return number of FLOPS for each call to Apply()
Teuchos::RCP< Epetra_Import > col_importer_trans
Importer from range map to column map.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
int max_num_mat_vec
Maximum number of matvecs in Apply.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Stores SG parallel communicator.
Abstract base class for multivariate orthogonal polynomials.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k 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.