Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MatrixFreeOperator.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
43#include "EpetraExt_BlockMultiVector.h"
44#include "EpetraExt_BlockUtility.h"
45#include "Teuchos_TimeMonitor.hpp"
46
49 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
50 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
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"),
58 sg_comm(sg_comm_),
59 sg_basis(sg_basis_),
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()),
66 global_col_map(),
67 global_col_map_trans(),
68 stoch_col_map(epetraCijk->getStochasticColMap()),
69 col_importer(),
70 col_importer_trans(),
71 Cijk(epetraCijk->getParallelCijk()),
72 block_ops(),
73 scale_op(true),
74 include_mean(true),
75 only_use_linear(false),
76 useTranspose(false),
77 use_block_apply(true),
78 expansion_size(sg_basis->size()),
79 num_blocks(0),
80 input_col(),
81 input_col_trans(),
82 input_block(),
83 result_block(),
84 tmp(),
85 tmp_trans(),
86 k_begin(Cijk->k_begin()),
87 k_end(Cijk->k_end())
88{
89 scale_op = params->get("Scale Operator by Inverse Basis Norms", true);
90 include_mean = params->get("Include Mean", true);
91 only_use_linear = params->get("Only Use Linear Terms", false);
92 use_block_apply = params->get("Use Block Apply", true);
93
94 // Compute maximum number of mat-vec's needed
95 if (!include_mean && index(k_begin) == 0)
96 ++k_begin;
97 if (only_use_linear) {
98 int dim = sg_basis->dimension();
99 k_end = Cijk->find_k(dim+1);
100 }
101 max_num_mat_vec = 0;
102 for (Cijk_type::k_iterator k=k_begin; k!=k_end; ++k) {
103 int nj = Cijk->num_j(k);
104 if (max_num_mat_vec < nj)
105 max_num_mat_vec = nj;
106 }
107
108 // Build up column map of SG operator
109 int num_col_blocks = expansion_size;
110 int num_row_blocks = expansion_size;
111 if (is_stoch_parallel) {
112
113 // Build column map from base domain map. This will communicate
114 // stochastic components to column map, but not deterministic. It would
115 // be more efficient to do both, but the Epetra_Operator interface
116 // doesn't have the concept of a column map.
118 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*domain_base_map,
120 *sg_comm));
122 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*range_base_map,
124 *sg_comm));
125
126 // Build importer from Domain Map to Column Map
128 Teuchos::rcp(new Epetra_Import(*global_col_map, *domain_sg_map));
130 Teuchos::rcp(new Epetra_Import(*global_col_map_trans, *range_sg_map));
131
132 num_col_blocks = epetraCijk->numMyCols();
133 num_row_blocks = epetraCijk->numMyRows();
134 }
135
136 input_block.resize(num_col_blocks);
137 result_block.resize(num_row_blocks);
138}
139
144
145double
147countApplyFlops() const
148{
149 int n_apply = 0;
150 int n_add = 0;
151 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
152 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
153 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
154 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
155 ++n_apply;
156 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
157 i_it != Cijk->i_end(j_it); ++i_it) {
158 ++n_add;
159 }
160 }
161 }
162
163 const Epetra_CrsMatrix& mat =
164 dynamic_cast<const Epetra_CrsMatrix&>((*block_ops)[0]);
165 const int nnz = mat.NumGlobalNonzeros();
166 const int nrow = mat.NumGlobalRows();
167 return 2.0 * static_cast<double>(n_apply) * static_cast<double>(nnz) +
168 static_cast<double>(n_add) * static_cast<double>(nrow);
169}
170
171void
174 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
175{
176 block_ops = ops;
177 num_blocks = block_ops->size();
178 if (num_blocks < Cijk->num_k())
179 k_end = Cijk->find_k(num_blocks);
180}
181
182Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
185{
186 return block_ops;
187}
188
189Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
191getSGPolynomial() const
192{
193 return block_ops;
194}
195
196int
198SetUseTranspose (bool UseTheTranspose)
199{
200 useTranspose = UseTheTranspose;
201 for (int i=0; i<num_blocks; i++)
202 (*block_ops)[i].SetUseTranspose(useTranspose);
203
204 return 0;
205}
206
207int
209Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
210{
211#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
212 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: SG Operator Apply()");
213#endif
214
215 // Note for transpose:
216 // The stochastic matrix is symmetric, however the matrix blocks may not
217 // be. So the algorithm here is the same whether we are using the transpose
218 // or not. We just apply the transpose of the blocks in the case of
219 // applying the global transpose, and make sure the imported Input
220 // vectors use the right map.
221
222 // We have to be careful if Input and Result are the same vector.
223 // If this is the case, the only possible solution is to make a copy
224 const Epetra_MultiVector *input = &Input;
225 bool made_copy = false;
226 if (Input.Values() == Result.Values() && !is_stoch_parallel) {
227 input = new Epetra_MultiVector(Input);
228 made_copy = true;
229 }
230
231 // Initialize
232 Result.PutScalar(0.0);
233
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();
239 }
240
241 // Allocate temporary storage
242 int m = Input.NumVectors();
243 if (useTranspose == false &&
244 (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec))
245 tmp = Teuchos::rcp(new Epetra_MultiVector(*result_base_map,
246 m*max_num_mat_vec));
247 else if (useTranspose == true &&
248 (tmp_trans == Teuchos::null ||
249 tmp_trans->NumVectors() != m*max_num_mat_vec))
250 tmp_trans = Teuchos::rcp(new Epetra_MultiVector(*result_base_map,
251 m*max_num_mat_vec));
252 Epetra_MultiVector *tmp_result;
253 if (useTranspose == false)
254 tmp_result = tmp.get();
255 else
256 tmp_result = tmp_trans.get();
257
258 // Map input into column map
259 const Epetra_MultiVector *tmp_col;
260 if (!is_stoch_parallel)
261 tmp_col = input;
262 else {
263 if (useTranspose == false) {
264 if (input_col == Teuchos::null || input_col->NumVectors() != m)
265 input_col = Teuchos::rcp(new Epetra_MultiVector(*global_col_map, m));
266 input_col->Import(*input, *col_importer, Insert);
267 tmp_col = input_col.get();
268 }
269 else {
270 if (input_col_trans == Teuchos::null ||
271 input_col_trans->NumVectors() != m)
272 input_col_trans =
273 Teuchos::rcp(new Epetra_MultiVector(*global_col_map_trans, m));
274 input_col_trans->Import(*input, *col_importer_trans, Insert);
275 tmp_col = input_col_trans.get();
276 }
277 }
278
279 // Extract blocks
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);
286
287 // Apply block SG operator via
288 // w_i =
289 // \sum_{j=0}^P \sum_{k=0}^L J_k v_j < \psi_i \psi_j \psi_k > / <\psi_i^2>
290 // for i=0,...,P where P = expansion_size, L = num_blocks, w_j is the jth
291 // input block, w_i is the ith result block, and J_k is the kth block operator
292
293 // k_begin and k_end are initialized in the constructor
294 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
295 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
296 int k = index(k_it);
297 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
298 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
299 int nj = Cijk->num_j(k_it);
300 if (nj > 0) {
301 Teuchos::Array<double*> j_ptr(nj*m);
302 Teuchos::Array<int> mj_indices(nj*m);
303 int l = 0;
304 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
305 int j = index(j_it);
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;
309 }
310 l++;
311 }
312 Epetra_MultiVector input_tmp(View, *input_base_map, &j_ptr[0], nj*m);
313 Epetra_MultiVector result_tmp(View, *tmp_result, &mj_indices[0], nj*m);
314 if (use_block_apply) {
315 (*block_ops)[k].Apply(input_tmp, result_tmp);
316 }
317 else {
318 for (int jj=0; jj<nj*m; jj++)
319 (*block_ops)[k].Apply(*(input_tmp(jj)), *(result_tmp(jj)));
320 }
321 l = 0;
322 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
323 int j = index(j_it);
324 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
325 i_it != Cijk->i_end(j_it); ++i_it) {
326 int i = index(i_it);
327 double c = value(i_it);
328 if (scale_op) {
329 int i_gid;
330 if (useTranspose)
331 i_gid = epetraCijk->GCID(j);
332 else
333 i_gid = epetraCijk->GRID(i);
334 c /= norms[i_gid];
335 }
336 for (int mm=0; mm<m; mm++)
337 (*result_block[i])(mm)->Update(c, *result_tmp(l*m+mm), 1.0);
338 }
339 l++;
340 }
341 }
342 }
343
344 // Destroy blocks
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;
349
350 if (made_copy)
351 delete input;
352
353 return 0;
354}
355
356int
358 Epetra_MultiVector& Result) const
359{
360 throw "MatrixFreeOperator::ApplyInverse not defined!";
361 return -1;
362}
363
364double
366{
367 return 1.0;
368}
369
370
371const char*
373{
374 return const_cast<char*>(label.c_str());
375}
376
377bool
379{
380 return useTranspose;
381}
382
383bool
385{
386 return false;
387}
388
389const Epetra_Comm &
391{
392 return *sg_comm;
393}
394const Epetra_Map&
396{
397 if (useTranspose)
398 return *range_sg_map;
399 return *domain_sg_map;
400}
401
402const Epetra_Map&
404{
405 if (useTranspose)
406 return *domain_sg_map;
407 return *range_sg_map;
408}
Insert
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.
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 > &params)
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.