Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_KLMatrixFreeOperator.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
48 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
49 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
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"),
57 sg_comm(sg_comm_),
58 sg_basis(sg_basis_),
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()),
65 global_col_map(),
66 global_col_map_trans(),
67 stoch_col_map(epetraCijk->getStochasticColMap()),
68 col_importer(),
69 col_importer_trans(),
70 Cijk(epetraCijk->getParallelCijk()),
71 block_ops(),
72 scale_op(true),
73 include_mean(true),
74 useTranspose(false),
75 expansion_size(sg_basis->size()),
76 num_blocks(0),
77 input_col(),
78 input_col_trans(),
79 input_block(),
80 result_block(),
81 tmp(),
82 tmp_trans(),
83 k_begin(Cijk->k_begin()),
84 k_end(Cijk->k_end())
85{
86 scale_op = params->get("Scale Operator by Inverse Basis Norms", true);
87 include_mean = params->get("Include Mean", true);
88
89 // Compute maximum number of mat-vec's needed
90 if (!include_mean && index(k_begin) == 0)
91 ++k_begin;
92 int dim = sg_basis->dimension();
93 k_end = Cijk->find_k(dim+1);
95 for (Cijk_type::k_iterator k=k_begin; k!=k_end; ++k) {
96 int nj = Cijk->num_j(k);
97 if (max_num_mat_vec < nj)
98 max_num_mat_vec = nj;
99 }
100
101 // Build up column map of SG operator
102 int num_col_blocks = expansion_size;
103 int num_row_blocks = expansion_size;
104 if (is_stoch_parallel) {
105
106 // Build column map from base domain map. This will communicate
107 // stochastic components to column map, but not deterministic. It would
108 // be more efficient to do both, but the Epetra_Operator interface
109 // doesn't have the concept of a column map.
111 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*domain_base_map,
113 *sg_comm));
115 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*range_base_map,
117 *sg_comm));
118
119 // Build importer from Domain Map to Column Map
121 Teuchos::rcp(new Epetra_Import(*global_col_map, *domain_sg_map));
123 Teuchos::rcp(new Epetra_Import(*global_col_map_trans, *range_sg_map));
124
125 num_col_blocks = epetraCijk->numMyCols();
126 num_row_blocks = epetraCijk->numMyRows();
127 }
128
129 input_block.resize(num_col_blocks);
130 result_block.resize(num_row_blocks);
131}
132
136
137void
140 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
141{
142 block_ops = ops;
143 num_blocks = sg_basis->dimension() + 1;
144}
145
146Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
149{
150 return block_ops;
151}
152
153Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
155getSGPolynomial() const
156{
157 return block_ops;
158}
159
160int
162SetUseTranspose(bool UseTheTranspose)
163{
164 useTranspose = UseTheTranspose;
165 for (int i=0; i<num_blocks; i++)
166 (*block_ops)[i].SetUseTranspose(useTranspose);
167
168 return 0;
169}
170
171int
173Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
174{
175 // We have to be careful if Input and Result are the same vector.
176 // If this is the case, the only possible solution is to make a copy
177 const Epetra_MultiVector *input = &Input;
178 bool made_copy = false;
179 if (Input.Values() == Result.Values() && !is_stoch_parallel) {
180 input = new Epetra_MultiVector(Input);
181 made_copy = true;
182 }
183
184 // Initialize
185 Result.PutScalar(0.0);
186
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();
192 }
193
194 // Allocate temporary storage
195 int m = Input.NumVectors();
196 if (useTranspose == false &&
197 (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec))
198 tmp = Teuchos::rcp(new Epetra_MultiVector(*result_base_map,
199 m*max_num_mat_vec));
200 else if (useTranspose == true &&
201 (tmp_trans == Teuchos::null ||
202 tmp_trans->NumVectors() != m*max_num_mat_vec))
203 tmp_trans = Teuchos::rcp(new Epetra_MultiVector(*result_base_map,
204 m*max_num_mat_vec));
205 Epetra_MultiVector *tmp_result;
206 if (useTranspose == false)
207 tmp_result = tmp.get();
208 else
209 tmp_result = tmp_trans.get();
210
211 // Map input into column map
212 const Epetra_MultiVector *tmp_col;
213 if (!is_stoch_parallel)
214 tmp_col = input;
215 else {
216 if (useTranspose == false) {
217 if (input_col == Teuchos::null || input_col->NumVectors() != m)
218 input_col = Teuchos::rcp(new Epetra_MultiVector(*global_col_map, m));
219 input_col->Import(*input, *col_importer, Insert);
220 tmp_col = input_col.get();
221 }
222 else {
223 if (input_col_trans == Teuchos::null ||
224 input_col_trans->NumVectors() != m)
225 input_col_trans =
226 Teuchos::rcp(new Epetra_MultiVector(*global_col_map_trans, m));
227 input_col_trans->Import(*input, *col_importer_trans, Insert);
228 tmp_col = input_col_trans.get();
229 }
230 }
231
232 // Extract blocks
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();
240
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++) {
245 zero[j] = 0.0;
246 one[j] = 1.0;
247 }
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);
251
252 // k_begin and k_end are initialized in the constructor
253 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
254 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
255 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
256 int k = index(k_it);
257 int nj = Cijk->num_j(k_it);
258 if (nj > 0) {
259 Teuchos::Array<double*> j_ptr(nj*m);
260 Teuchos::Array<int> mj_indices(nj*m);
261 int l = 0;
262 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
263 int j = index(j_it);
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;
267 }
268 l++;
269 }
270 Epetra_MultiVector input_tmp(View, *input_base_map, &j_ptr[0], nj*m);
271 Epetra_MultiVector result_tmp(View, *tmp_result, &mj_indices[0], nj*m);
272 (*block_ops)[k].Apply(input_tmp, result_tmp);
273 l = 0;
274 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
275 int j = index(j_it);
276 int j_gid = epetraCijk->GCID(j);
277 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
278 i_it != Cijk->i_end(j_it); ++i_it) {
279 int i = index(i_it);
280 int i_gid = epetraCijk->GRID(i);
281 double c = value(i_it);
282 if (k == 0)
283 c /= phi_0[0];
284 else {
285 c /= phi_1[k];
286 if (i_gid == j_gid)
287 c -= phi_0[k]/(phi_1[k]*phi_0[0])*norms[i_gid];
288 }
289 if (scale_op) {
290 if (useTranspose)
291 c /= norms[j_gid];
292 else
293 c /= norms[i_gid];
294 }
295 for (int mm=0; mm<m; mm++)
296 (*result_block[i])(mm)->Update(c, *result_tmp(l*m+mm), 1.0);
297 }
298 l++;
299 }
300 }
301 }
302
303 // Destroy blocks
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;
308
309 if (made_copy)
310 delete input;
311
312 return 0;
313}
314
315int
317ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
318{
319 throw "KLMatrixFreeOperator::ApplyInverse not defined!";
320 return -1;
321}
322
323double
325NormInf() const
326{
327 return 1.0;
328}
329
330
331const char*
333Label () const
334{
335 return const_cast<char*>(label.c_str());
336}
337
338bool
340UseTranspose() const
341{
342 return useTranspose;
343}
344
345bool
347HasNormInf() const
348{
349 return false;
350}
351
352const Epetra_Comm &
354Comm() const
355{
356 return *sg_comm;
357}
358const Epetra_Map&
360OperatorDomainMap() const
361{
362 if (useTranspose)
363 return *range_sg_map;
364 return *domain_sg_map;
365}
366
367const Epetra_Map&
369OperatorRangeMap() const
370{
371 if (useTranspose)
372 return *domain_sg_map;
373 return *range_sg_map;
374}
Insert
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 > &params)
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 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.
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.