Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_AdaptivityManager.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
47
48#include "EpetraExt_BlockVector.h"
49#include "EpetraExt_RowMatrixOut.h"
50
51#ifdef HAVE_STOKHOS_BOOST
52#include <boost/functional/hash.hpp>
53#endif
54
55#ifdef HAVE_STOKHOS_BOOST // we have boost, use the hash Stokhos, use the hash!
56std::size_t Stokhos::AdaptivityManager::Sparse3TensorHash::IJKHash::
57operator()(const Stokhos::AdaptivityManager::Sparse3TensorHash::IJK & ijk) const
58{
59 std::size_t seed=0;
60 boost::hash_combine(seed,ijk.i_);
61 boost::hash_combine(seed,ijk.j_);
62 boost::hash_combine(seed,ijk.k_);
63 return seed;
64}
65
67{
71
72 for(k_iterator k_it = Cijk.k_begin();k_it!=Cijk.k_end();k_it++) {
73 int k = *k_it;
74 for(kj_iterator j_it = Cijk.j_begin(k_it);j_it!=Cijk.j_end(k_it);j_it++) {
75 int j = *j_it;
76 for(kji_iterator i_it = Cijk.i_begin(j_it);i_it!=Cijk.i_end(j_it);i_it++) {
77 int i = *i_it;
78 hashMap_[IJK(i,j,k)] = i_it.value();
79 }
80 }
81 }
82}
83
85{
86 boost::unordered_map<IJK,double>::const_iterator itr;
87 itr = hashMap_.find(IJK(i,j,k));
88
89 if(itr==hashMap_.end()) return 0.0;
90
91 return itr->second;
92}
93#else // no BOOST, just default to the slow thing
95 : Cijk_(Cijk)
96{ }
97
99{
100 return Cijk_.getValue(i,j,k);
101}
102#endif
103
105 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> >& sg_master_basis,
106 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & sg_basis_row_dof,
107 const Epetra_CrsGraph & determ_graph,
108 bool onlyUseLinear,int kExpOrder,
109 bool scaleOp)
110 : sg_master_basis_(sg_master_basis), sg_basis_row_dof_(sg_basis_row_dof), scaleOp_(scaleOp)
111{
113
114 setupWithGraph(determ_graph,onlyUseLinear,kExpOrder);
115}
116
118 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> >& sg_master_basis,
119 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & sg_basis_row_dof,
120 const Epetra_Comm & comm,
121 bool scaleOp)
122 : sg_master_basis_(sg_master_basis), sg_basis_row_dof_(sg_basis_row_dof), scaleOp_(scaleOp)
123{
125}
126
127Teuchos::RCP<Epetra_CrsMatrix>
129{
130 return Teuchos::rcp(new Epetra_CrsMatrix(Copy,*graph_));
131}
132
133void Stokhos::AdaptivityManager::setupWithGraph(const Epetra_CrsGraph & determGraph,bool onlyUseLinear,int kExpOrder)
134{
135 graph_ = adapt_utils::buildAdaptedGraph(determGraph, sg_master_basis_, sg_basis_row_dof_, onlyUseLinear, kExpOrder);
136
137 adapt_utils::buildAdaptedColOffsets(determGraph,myRowGidOffsets_,myColGidOffsets_);
138 adapt_utils::buildColBasisFunctions(determGraph,sg_master_basis_,sg_basis_row_dof_,sg_basis_col_dof_);
139}
140
143void
146 bool onlyUseLinear,bool includeMean) const
147{
148 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
149
150 // build the sparse hash only once
151 Sparse3TensorHash hashLookup(Cijk);
152
153 // Zero out matrix
154 A.PutScalar(0.0);
155
156 // Compute loop bounds
157 Cijk_type::k_iterator k_begin = Cijk.k_begin();
158 Cijk_type::k_iterator k_end = Cijk.k_end();
159 if (!includeMean && index(k_begin) == 0)
160 ++k_begin;
161 if (onlyUseLinear) {
162 int dim = sg_master_basis_->dimension();
163 k_end = Cijk.find_k(dim+1);
164 }
165
166 // Assemble matrix
167 // const Teuchos::Array<double>& norms = sg_basis->norm_squared();
168 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
169 int k = index(k_it);
170 Teuchos::RCP<Epetra_CrsMatrix> block =
171 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(poly.getCoeffPtr(k),
172 true);
173
174 // add in matrix k
175 sumInOperator(A,hashLookup,k,*block);
176 }
177}
178
179void
182{
183 // this allows the simple interface of taking a Sparse3Tensor but immediately computes
184 // the sparse hash (if boost is enabled)
185
186 Sparse3TensorHash hashLookup(Cijk);
187 sumInOperator(A,hashLookup,k,J_k);
188}
189
190void
193{
194 TEUCHOS_ASSERT(J_k.NumMyRows() == int(sg_basis_row_dof_.size()));
195 TEUCHOS_ASSERT(J_k.NumMyCols() == int(sg_basis_col_dof_.size()));
196
197 const Teuchos::Array<double> & normValues = sg_master_basis_->norm_squared();
198
199 // loop over deterministic rows
200 for(int localM=0;localM<J_k.NumMyRows();localM++) {
201 // int m = J_k.GRID(localM); // unused
202
203 // grab row basis
204 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > rowStochBasis
205 = sg_basis_row_dof_[localM];
206
207 // grab row from deterministic system
208 int d_numEntries;
209 int * d_Indices;
210 double * d_Values;
211
212 J_k.ExtractMyRowView(localM,d_numEntries,d_Values,d_Indices);
213
214 // loop over stochastic degrees of freedom of this row
215 for(int rb_i=0;rb_i<rowStochBasis->size();rb_i++) {
216 int i = sg_master_basis_->index(rowStochBasis->term(rb_i));
217
218 double normValue = normValues[i]; // sg_master_basis->norm_squared(i);
219
220 int sg_m = getGlobalRowId(localM,rb_i);
221
222 // we wipe out old values, capacity should gurantee
223 // we don't allocate more often than neccessary!
224 std::vector<int> sg_indices;
225 std::vector<double> sg_values;
226
227 // sg_indices.resize(0);
228 // sg_values.resize(0);
229
230 // loop over each column
231 for(int colInd=0;colInd<d_numEntries;colInd++) {
232 int localN = d_Indices[colInd]; // grab local deterministic column id
233
234 // grab row basis
235 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > colStochBasis
236 = sg_basis_col_dof_[localN];
237
238 // build values array
239 for(int cb_j=0;cb_j<colStochBasis->size();cb_j++) {
240 int j = sg_master_basis_->index(colStochBasis->term(cb_j));
241 int sg_n = getGlobalColId(localN,cb_j);
242 double cijk = Cijk.getValue(i,j,k);
243
244 // no reason to work it in!
245 if(cijk==0) continue;
246
247 if(scaleOp_)
248 cijk = cijk/normValue;
249
250 sg_indices.push_back(sg_n);
251 sg_values.push_back(cijk*d_Values[colInd]);
252 }
253 }
254
255 // add in matrix values
256 A.SumIntoGlobalValues(sg_m,sg_indices.size(),&sg_values[0],&sg_indices[0]);
257 }
258 }
259}
260
264{
265 Teuchos::RCP<const EpetraExt::BlockVector> x_sg_bv = x_sg.getBlockVector();
266
267 // copy from adapted vector to deterministic
268 for(std::size_t i=0;i<sg_basis_row_dof_.size();i++) {
269 int P_i = getRowStochasticBasisSize(i);
270 int localId = rowMap_->LID(getGlobalRowId(i,0));
271
272 for(int j=0;j<P_i;j++,localId++) {
273 int blk = sg_master_basis_->index(sg_basis_row_dof_[i]->term(j));
274 x[localId] = x_sg_bv->GetBlock(blk)->operator[](i);
275 }
276 }
277}
278
282{
283 int numBlocks = x_sg.size();
284 Teuchos::RCP<EpetraExt::BlockVector> x_sg_bv = x_sg.getBlockVector();
285
286 // zero out determinstic vectors
287 for(int blk=0;blk<numBlocks;blk++)
288 x_sg_bv->GetBlock(blk)->PutScalar(0.0);
289
290 // copy from adapted vector to deterministic
291 for(std::size_t i=0;i<sg_basis_row_dof_.size();i++) {
292 int P_i = getRowStochasticBasisSize(i);
293 int localId = rowMap_->LID(getGlobalRowId(i,0));
294
295 for(int j=0;j<P_i;j++,localId++) {
296 int blk = sg_master_basis_->index(sg_basis_row_dof_[i]->term(j));
297 x_sg_bv->GetBlock(blk)->operator[](i) = x[localId];
298 }
299 }
300}
Copy
const Epetra_Comm & Comm() const
int NumMyCols() const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int PutScalar(double ScalarConstant)
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int NumMyRows() const
Sparse3TensorHash(const Stokhos::Sparse3Tensor< int, double > &Cijk)
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > sg_master_basis_
Teuchos::RCP< Epetra_CrsMatrix > buildMatrixFromGraph() const
Teuchos::RCP< Epetra_Map > rowMap_
void setupOperator(Epetra_CrsMatrix &A, const Sparse3Tensor< int, double > &Cijk, Stokhos::EpetraOperatorOrthogPoly &poly, bool onlyUseLinear=false, bool includeMean=true) const
void copyToAdaptiveVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x) const
void copyFromAdaptiveVector(const Epetra_Vector &x, Stokhos::EpetraVectorOrthogPoly &x_sg) const
void setupWithGraph(const Epetra_CrsGraph &graph, bool onlyUseLinear, int kExpOrder)
void sumInOperator(Epetra_CrsMatrix &A, const Stokhos::Sparse3Tensor< int, double > &Cijk, int k, const Epetra_CrsMatrix &J_k) const
AdaptivityManager(const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &sg_master_basis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &sg_basis_row_dof, const Epetra_CrsGraph &determ_graph, bool onlyUseLinear, int kExpOrder, bool scaleOp=true)
std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > sg_basis_row_dof_
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Teuchos::RCP< coeff_type > getCoeffPtr(ordinal_type i)
Return ref-count pointer to coefficient i.
ordinal_type size() const
Return size.
Teuchos::RCP< EpetraExt::BlockVector > getBlockVector()
Get block vector.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
k_iterator find_k(ordinal_type k) const
Return k iterator for given index k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
k_iterator k_end() const
Iterator pointing to last k entry.
Teuchos::RCP< Epetra_Map > buildAdaptedRowMapAndOffsets(const Epetra_Comm &Comm, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, std::vector< int > &myRowGidOffsets)
void buildColBasisFunctions(const Epetra_CrsGraph &determGraph, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &masterBasis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_col_basis)
Teuchos::RCP< Epetra_CrsGraph > buildAdaptedGraph(const Epetra_CrsGraph &determGraph, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &masterBasis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, bool onlyUseLinear=false, int kExpOrder=-1)
void buildAdaptedColOffsets(const Epetra_CrsGraph &determGraph, const std::vector< int > &myRowGidOffsets, std::vector< int > &myColGidOffsets)
Bi-directional iterator for traversing a sparse array.