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,
178 bool made_copy =
false;
179 if (Input.Values() == Result.Values() && !is_stoch_parallel) {
185 Result.PutScalar(0.0);
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();
195 int m = Input.NumVectors();
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;
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.