Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_EpetraSparse3Tensor.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
42#include "Stokhos_ConfigDefs.h"
45#include "Teuchos_Assert.hpp"
46#ifdef HAVE_STOKHOS_ISORROPIA
47#include "Isorropia_Epetra.hpp"
48#endif
49#include "Epetra_Map.h"
50
53 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
54 const Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> >& Cijk_,
55 const Teuchos::RCP<const EpetraExt::MultiComm>& globalMultiComm_,
56 int k_begin_, int k_end_) :
57 basis(basis_),
58 Cijk(Cijk_),
59 globalMultiComm(globalMultiComm_),
60 k_begin(k_begin_),
61 k_end(k_end_)
62{
63 // Get stochastic distribution
65 stoch_comm = Teuchos::rcp(&(globalMultiComm->TimeDomainComm()), false);
66 is_parallel = (stoch_comm->NumProc() > 1);
67 if (k_end == -1)
68 k_end = Cijk->num_k();
69
70 // Build stochastic row map -- stochastic blocks evenly distributed
71 // across stochastic procs
73 *stoch_comm));
74
75 // Build Cijk tensor parallel over i
76 if (!is_parallel)
78 else
80
81 // Build stochastic graph
83
84 // Build stochastic column map
85 stoch_col_map = Teuchos::rcp(&(stoch_graph->ColMap()), false);
86}
87
90 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
91 const Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> >& Cijk_,
92 const Teuchos::RCP<const EpetraExt::MultiComm>& globalMultiComm_,
93 const Teuchos::RCP<const Epetra_BlockMap>& stoch_row_map_,
94 const Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> >& Cijk_parallel_,
95 int k_begin_, int k_end_) :
96 basis(basis_),
97 Cijk(Cijk_),
98 globalMultiComm(globalMultiComm_),
99 k_begin(k_begin_),
100 k_end(k_end_),
101 stoch_row_map(stoch_row_map_),
102 Cijk_parallel(Cijk_parallel_)
103{
104 // Get stochastic distribution
106 stoch_comm = Teuchos::rcp(&(globalMultiComm->TimeDomainComm()), false);
107 is_parallel = (stoch_comm->NumProc() > 1);
108 if (k_end == -1)
109 k_end = Cijk->num_k();
110
111 // Build Cijk tensor parallel over i
112 if (Cijk_parallel == Teuchos::null) {
113 if (!is_parallel)
115 else
117 }
118
119 // Build stochastic graph
121
122 // Build stochastic column map
123 stoch_col_map = Teuchos::rcp(&(stoch_graph->ColMap()), false);
124}
125
128 int k_begin_, int k_end_) :
129 basis(epetraCijk.basis),
130 Cijk(epetraCijk.Cijk),
131 globalMultiComm(epetraCijk.globalMultiComm),
132 num_global_stoch_blocks(epetraCijk.num_global_stoch_blocks),
133 k_begin(k_begin_),
134 k_end(k_end_),
135 stoch_comm(epetraCijk.stoch_comm),
136 is_parallel(epetraCijk.is_parallel),
137 stoch_row_map(epetraCijk.stoch_row_map)
138{
139 if (k_end == -1)
140 k_end = Cijk->num_k();
141
142 // Build Cijk tensor parallel over i
143 if (!is_parallel)
145 else
147
148 // Build stochastic graph
150
151 // Build stochastic column map
152 stoch_col_map = Teuchos::rcp(&(stoch_graph->ColMap()), false);
153}
154
155void
157rebalance(Teuchos::ParameterList& isorropia_params)
158{
159#ifdef HAVE_STOKHOS_ISORROPIA
160 // Reblance with Isorropia
161 Teuchos::RCP<const Epetra_CrsGraph> rebalanced_graph =
162 Teuchos::rcp(Isorropia::Epetra::createBalancedCopy(
163 *stoch_graph, isorropia_params));
164 stoch_graph = rebalanced_graph;
165
166 // Get new maps
167 stoch_row_map = Teuchos::rcp(&(stoch_graph->RowMap()), false);
168 stoch_col_map = Teuchos::rcp(&(stoch_graph->ColMap()), false);
169
170 // Build new parallel Cijk
171 Cijk_parallel = buildParallelCijk();
172#else
173 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
174 "Error! Stokhos must be configured with " <<
175 "isorropia support to rebalance the stochastic " <<
176 "graph.");
177#endif
178}
179
180Teuchos::RCP<Stokhos::EpetraSparse3Tensor::Cijk_type>
182buildParallelCijk() const
183{
184 Teuchos::RCP<Cijk_type> Cijk_tmp = Teuchos::rcp(new Cijk_type);
185 int *rowIndices = stoch_row_map->MyGlobalElements();
186 int myBlockRows = stoch_row_map->NumMyElements();
187 for (int i=0; i<myBlockRows; i++) {
188 int myRow = rowIndices[i];
189 Cijk_type::i_iterator i_it = Cijk->find_i(myRow);
190 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
191 k_it != Cijk->k_end(i_it); ++k_it) {
192 int k = index(k_it);
193 if (k >= k_begin && k < k_end) {
194 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
195 j_it != Cijk->j_end(k_it); ++j_it) {
196 int j = index(j_it);
197 double c = value(j_it);
198 Cijk_tmp->add_term(myRow, j, k, c);
199 }
200 }
201 }
202 }
203 Cijk_tmp->fillComplete();
204
205 return Cijk_tmp;
206}
207
208void
211{
212 Teuchos::RCP<Cijk_type> Cijk_tmp = Teuchos::rcp(new Cijk_type);
213 for (Cijk_type::i_iterator i_it = Cijk_parallel->i_begin();
214 i_it != Cijk_parallel->i_end(); ++i_it) {
215 int i = stoch_row_map->LID(index(i_it));
216 TEUCHOS_TEST_FOR_EXCEPTION(i == -1, std::logic_error,
217 "Error! Global row index " << index(i_it) <<
218 " is not on processor " << globalMultiComm->MyPID());
219 for (Cijk_type::ik_iterator k_it = Cijk_parallel->k_begin(i_it);
220 k_it != Cijk_parallel->k_end(i_it); ++k_it) {
221 int k = index(k_it);
222 for (Cijk_type::ikj_iterator j_it = Cijk_parallel->j_begin(k_it);
223 j_it != Cijk_parallel->j_end(k_it); ++j_it) {
224 int j= stoch_col_map->LID(index(j_it));
225 TEUCHOS_TEST_FOR_EXCEPTION(j == -1, std::logic_error,
226 "Error! Global col index " << index(j_it) <<
227 " is not on processor " << globalMultiComm->MyPID());
228 double c = value(j_it);
229 Cijk_tmp->add_term(i, j, k, c);
230 }
231 }
232 }
233 Cijk_tmp->fillComplete();
234
235 Cijk_parallel = Cijk_tmp;
236}
Teuchos::RCP< const Epetra_CrsGraph > stoch_graph
Stochastic operator graph.
void transformToLocal()
Transform Cijk to local i and j indices.
EpetraSparse3Tensor(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &basis, const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &Cijk, const Teuchos::RCP< const EpetraExt::MultiComm > &globalMultiComm, int k_begin=0, int k_end=-1)
Constructor from a full Cijk.
void rebalance(Teuchos::ParameterList &isorropia_params)
Rebalance maps and graph using Isorropia.
Teuchos::RCP< Cijk_type > buildParallelCijk() const
Build parallel Cijk tensor from a parallel row map.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > basis
Basis.
int num_global_stoch_blocks
Number of global stochastic blocks.
Teuchos::RCP< const Cijk_type > Cijk_parallel
Cijk tensor parallel over i.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Stochastic row-map.
Teuchos::RCP< const EpetraExt::MultiComm > globalMultiComm
Multi-comm.
Teuchos::RCP< const Cijk_type > Cijk
Triple product.
Teuchos::RCP< const Epetra_BlockMap > stoch_col_map
Stochastic col-map.
bool is_parallel
Whether stochastic blocks are parallel.
Teuchos::RCP< const Epetra_Comm > stoch_comm
Stochastic comm.
Abstract base class for multivariate orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
j_sparse_array::const_iterator ikj_iterator
Iterator for looping over j entries given i and k.
kj_sparse_array::const_iterator ik_iterator
Iterator for looping over k entries given i.
ikj_sparse_array::const_iterator i_iterator
Iterator for looping over i entries.
Teuchos::RCP< Epetra_CrsGraph > sparse3Tensor2CrsGraph(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm)
Build an Epetra_CrsGraph from a sparse 3 tensor.