43#include "EpetraExt_BlockMultiVector.h"
44#include "EpetraExt_BlockUtility.h"
48 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
50 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
51 const Teuchos::RCP<const Epetra_Map>& domain_base_map_,
52 const Teuchos::RCP<const Epetra_Map>& range_base_map_,
53 const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
54 const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
55 const Teuchos::RCP<Teuchos::ParameterList>& params) :
56 label(
"Stokhos KL MatrixFree Operator"),
59 epetraCijk(epetraCijk_),
60 domain_base_map(domain_base_map_),
61 range_base_map(range_base_map_),
62 domain_sg_map(domain_sg_map_),
63 range_sg_map(range_sg_map_),
64 is_stoch_parallel(epetraCijk->isStochasticParallel()),
66 global_col_map_trans(),
67 stoch_col_map(epetraCijk->getStochasticColMap()),
70 Cijk(epetraCijk->getParallelCijk()),
75 expansion_size(sg_basis->size()),
83 k_begin(Cijk->k_begin()),
86 scale_op = params->get(
"Scale Operator by Inverse Basis Norms",
true);
96 int nj =
Cijk->num_j(k);
111 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*
domain_base_map,
115 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*
range_base_map,
140 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
143 num_blocks = sg_basis->dimension() + 1;
146Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
153Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
164 useTranspose = UseTheTranspose;
165 for (
int i=0; i<num_blocks; i++)
166 (*block_ops)[i].SetUseTranspose(useTranspose);
178 bool made_copy =
false;
179 if (Input.
Values() == Result.
Values() && !is_stoch_parallel) {
187 const Epetra_Map* input_base_map = domain_base_map.get();
188 const Epetra_Map* result_base_map = range_base_map.get();
189 if (useTranspose ==
true) {
190 input_base_map = range_base_map.get();
191 result_base_map = domain_base_map.get();
196 if (useTranspose ==
false &&
197 (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec))
200 else if (useTranspose ==
true &&
201 (tmp_trans == Teuchos::null ||
202 tmp_trans->NumVectors() != m*max_num_mat_vec))
206 if (useTranspose ==
false)
207 tmp_result = tmp.get();
209 tmp_result = tmp_trans.get();
213 if (!is_stoch_parallel)
216 if (useTranspose ==
false) {
217 if (input_col == Teuchos::null || input_col->NumVectors() != m)
219 input_col->Import(*input, *col_importer,
Insert);
220 tmp_col = input_col.get();
223 if (input_col_trans == Teuchos::null ||
224 input_col_trans->NumVectors() != m)
227 input_col_trans->Import(*input, *col_importer_trans,
Insert);
228 tmp_col = input_col_trans.get();
233 EpetraExt::BlockMultiVector sg_input(
View, *input_base_map, *tmp_col);
234 EpetraExt::BlockMultiVector sg_result(
View, *result_base_map, Result);
235 for (
int i=0; i<input_block.size(); i++)
236 input_block[i] = sg_input.GetBlock(i);
237 for (
int i=0; i<result_block.size(); i++)
238 result_block[i] = sg_result.GetBlock(i);
239 int N = result_block[0]->MyLength();
241 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
242 int d = sg_basis->dimension();
243 Teuchos::Array<double> zero(d), one(d);
244 for(
int j = 0;
j<d;
j++) {
248 Teuchos::Array< double > phi_0(expansion_size), phi_1(expansion_size);
249 sg_basis->evaluateBases(zero, phi_0);
250 sg_basis->evaluateBases(one, phi_1);
257 int nj = Cijk->num_j(k_it);
259 Teuchos::Array<double*> j_ptr(nj*m);
260 Teuchos::Array<int> mj_indices(nj*m);
264 for (
int mm=0; mm<m; mm++) {
265 j_ptr[l*m+mm] = input_block[
j]->Values()+mm*N;
266 mj_indices[l*m+mm] = l*m+mm;
272 (*block_ops)[k].Apply(input_tmp, result_tmp);
276 int j_gid = epetraCijk->GCID(
j);
278 i_it != Cijk->i_end(j_it); ++i_it) {
280 int i_gid = epetraCijk->GRID(i);
281 double c = value(i_it);
287 c -= phi_0[k]/(phi_1[k]*phi_0[0])*norms[i_gid];
295 for (
int mm=0; mm<m; mm++)
296 (*result_block[i])(mm)->
Update(c, *result_tmp(l*m+mm), 1.0);
304 for (
int i=0; i<input_block.size(); i++)
305 input_block[i] = Teuchos::null;
306 for (
int i=0; i<result_block.size(); i++)
307 result_block[i] = Teuchos::null;
319 throw "KLMatrixFreeOperator::ApplyInverse not defined!";
335 return const_cast<char*
>(label.c_str());
363 return *range_sg_map;
364 return *domain_sg_map;
372 return *domain_sg_map;
373 return *range_sg_map;
int PutScalar(double ScalarConstant)
Cijk_type::k_iterator k_begin
Starting k iterator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
KLMatrixFreeOperator(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< Epetra_Import > col_importer
Importer from domain map to column map.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Stores SG parallel communicator.
Teuchos::RCP< const Epetra_Map > domain_base_map
Stores domain base map.
bool include_mean
Flag indicating whether to include mean term.
Teuchos::RCP< const Epetra_Map > domain_sg_map
Stores domain SG map.
Teuchos::RCP< Epetra_Map > global_col_map
Stores operator column SG map.
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.
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 ...
Teuchos::RCP< const Epetra_Map > range_base_map
Stores range base map.
int max_num_mat_vec
Maximum number of matvecs in Apply.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
Teuchos::RCP< const Cijk_type > Cijk
Stores triple product tensor.
virtual ~KLMatrixFreeOperator()
Destructor.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
Cijk_type::k_iterator k_end
Ending k iterator.
bool is_stoch_parallel
Whether we have parallelism over stochastic blocks.
Teuchos::RCP< Epetra_Import > col_importer_trans
Importer from range map to column map.
Teuchos::RCP< const Epetra_Map > range_sg_map
Stores range SG map.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const char * Label() const
Returns a character std::string describing the operator.
Teuchos::Array< Teuchos::RCP< Epetra_MultiVector > > result_block
MultiVectors for each block for Apply() result.
int expansion_size
Number of terms in expansion.
Teuchos::RCP< Epetra_Map > global_col_map_trans
Stores operator column SG map for transpose.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
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 const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.