Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cijk_partition.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_Epetra.hpp"
43#include "Teuchos_CommandLineProcessor.hpp"
44#include "Teuchos_ParameterList.hpp"
45#include "Isorropia_Epetra.hpp"
46#include "Isorropia_EpetraPartitioner.hpp"
47
48#ifdef HAVE_MPI
49#include "Epetra_MpiComm.h"
50#else
51#include "Epetra_SerialComm.h"
52#endif
53
54// sparsity_example
55//
56// usage:
57// sparsity_example [options]
58//
59// output:
60// prints the sparsity of the sparse 3 tensor specified by the basis,
61// dimension, order, given by summing over the third index, to a matrix
62// market file. This sparsity pattern yields the sparsity of the block
63// stochastic Galerkin matrix which can be visualized, e.g., by matlab.
64// The full/linear flag determines whether the third index ranges over
65// the full polynomial dimension, or only over the zeroth and first order
66// terms.
67
68// Basis types
70const int num_basis_types = 6;
73const char *basis_type_names[] = {
74 "hermite", "legendre", "clenshaw-curtis", "gauss-patterson", "rys", "jacobi" };
75
76// Growth policies
77const int num_growth_types = 2;
80const char *growth_type_names[] = { "slow", "moderate" };
81
82// Product Basis types
87const char *prod_basis_type_names[] = {
88 "complete", "tensor", "total", "smolyak" };
89
90// Ordering types
92const int num_ordering_types = 2;
95const char *ordering_type_names[] = {
96 "total", "lexicographic" };
97
98// Partitioning types
103const char *partitioning_type_names[] = {
104 "rcb", "hg_matrix" };
105
106using Teuchos::rcp;
107using Teuchos::RCP;
108using Teuchos::ParameterList;
109using Teuchos::Array;
110
111int main(int argc, char **argv)
112{
113 try {
114
115 // Initialize MPI
116#ifdef HAVE_MPI
117 MPI_Init(&argc,&argv);
118#endif
119
120 // Setup command line options
121 Teuchos::CommandLineProcessor CLP;
122 CLP.setDocString(
123 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
124 int d = 3;
125 CLP.setOption("dimension", &d, "Stochastic dimension");
126 int p = 5;
127 CLP.setOption("order", &p, "Polynomial order");
128 double drop = 1.0e-12;
129 CLP.setOption("drop", &drop, "Drop tolerance");
130 std::string file = "A.mm";
131 CLP.setOption("filename", &file, "Matrix Market filename");
133 CLP.setOption("basis", &basis_type,
135 "Basis type");
137 CLP.setOption("growth", &growth_type,
139 "Growth type");
140 ProductBasisType prod_basis_type = TOTAL;
141 CLP.setOption("product_basis", &prod_basis_type,
144 "Product basis type");
145 OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
146 CLP.setOption("ordering", &ordering_type,
149 "Product basis ordering");
150 PartitioningType partitioning_type = RCB;
151 CLP.setOption("partitioning", &partitioning_type,
154 "Partitioning Method");
155 double imbalance_tol = 1.0;
156 CLP.setOption("imbalance", &imbalance_tol, "Imbalance tolerance");
157 double alpha = 1.0;
158 CLP.setOption("alpha", &alpha, "Jacobi alpha index");
159 double beta = 1.0;
160 CLP.setOption("beta", &beta, "Jacobi beta index");
161 bool full = true;
162 CLP.setOption("full", "linear", &full, "Use full or linear expansion");
163 int tile_size = 128;
164 CLP.setOption("tile_size", &tile_size, "Tile size");
165 bool save_3tensor = false;
166 CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
167 "Save full 3tensor to file");
168 std::string file_3tensor = "Cijk.dat";
169 CLP.setOption("filename_3tensor", &file_3tensor,
170 "Filename to store full 3-tensor");
171
172 // Parse arguments
173 CLP.parse( argc, argv );
174
175 // Basis
176 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
177 for (int i=0; i<d; i++) {
178 if (basis_type == HERMITE)
179 bases[i] = rcp(new Stokhos::HermiteBasis<int,double>(
180 p, true, growth_type));
181 else if (basis_type == LEGENDRE)
182 bases[i] = rcp(new Stokhos::LegendreBasis<int,double>(
183 p, true, growth_type));
184 else if (basis_type == CC_LEGENDRE)
185 bases[i] =
187 p, true));
188 else if (basis_type == GP_LEGENDRE)
189 bases[i] =
191 p, true));
192 else if (basis_type == RYS)
193 bases[i] = rcp(new Stokhos::RysBasis<int,double>(
194 p, 1.0, true, growth_type));
195 else if (basis_type == JACOBI)
196 bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
197 p, alpha, beta, true, growth_type));
198 }
199 RCP<const Stokhos::ProductBasis<int,double> > basis;
202 if (prod_basis_type == COMPLETE)
203 basis =
205 bases, drop));
206 else if (prod_basis_type == TENSOR) {
207 if (ordering_type == TOTAL_ORDERING)
208 basis =
210 bases, drop));
211 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
212 basis =
214 bases, drop));
215 }
216
217 else if (prod_basis_type == TOTAL) {
218 if (ordering_type == TOTAL_ORDERING)
219 basis =
221 bases, drop));
222 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
223 basis =
225 bases, drop));
226 }
227 else if (prod_basis_type == SMOLYAK) {
228 Stokhos::TotalOrderIndexSet<int> index_set(d, p);
229 if (ordering_type == TOTAL_ORDERING)
230 basis =
232 bases, index_set, drop));
233 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
234 basis =
236 bases, index_set, drop));
237 }
238
239 // Triple product tensor
240 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
241 RCP<Cijk_type> Cijk;
242 if (full)
243 Cijk = basis->computeTripleProductTensor();
244 else
245 Cijk = basis->computeLinearTripleProductTensor();
246
247 int basis_size = basis->size();
248 std::cout << "basis size = " << basis_size
249 << " num nonzero Cijk entries = " << Cijk->num_entries()
250 << std::endl;
251
252#ifdef HAVE_MPI
253 Epetra_MpiComm comm(MPI_COMM_WORLD);
254#else
256#endif
257
258 // File for saving sparse Cijk tensor and parts
259 std::ofstream cijk_file;
260 if (save_3tensor) {
261 cijk_file.open(file_3tensor.c_str());
262 cijk_file.precision(14);
263 cijk_file.setf(std::ios::scientific);
264 cijk_file << "i, j, k, part" << std::endl;
265 }
266
267 Teuchos::Array<int> parts;
268 if (partitioning_type == RCB) {
269 // Store Cijk (i,j,k) triples in Epetra_MultiVector
270 int num_cijk_entries = Cijk->num_entries();
271 Epetra_LocalMap cijk_map(num_cijk_entries, 0, comm);
272 Epetra_MultiVector ijk_triples(cijk_map, 3);
273 int idx = 0;
274 Cijk_type::k_iterator k_begin = Cijk->k_begin();
275 Cijk_type::k_iterator k_end = Cijk->k_end();
276 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
277 int k = index(k_it);
278 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
279 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
280 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
281 int j = index(j_it);
282 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
283 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
284 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
285 int i = index(i_it);
286 ijk_triples[0][idx] = i;
287 ijk_triples[1][idx] = j;
288 ijk_triples[2][idx] = k;
289 ++idx;
290 }
291 }
292 }
293
294 // Partition ijk_triples using isorropia
295 ParameterList params;
296 params.set("partitioning method", "rcb");
297 int num_parts = num_cijk_entries / tile_size;
298 if (num_cijk_entries % tile_size > 0)
299 ++num_parts;
300 std::stringstream ss;
301 ss << num_parts;
302 params.set<std::string>("num parts", ss.str());
303 RCP<const Epetra_MultiVector> ijk_triples_rcp =
304 rcp(&ijk_triples,false);
305 Isorropia::Epetra::Partitioner partitioner(ijk_triples_rcp, params);
306 partitioner.compute();
307
308 std::cout << "num parts requested = " << num_parts
309 << " num parts computed = " << partitioner.numProperties() << std::endl;
310
311 parts.resize(num_cijk_entries);
312 int sz;
313 partitioner.extractPartsCopy(num_cijk_entries, sz, &parts[0]);
314 TEUCHOS_ASSERT(sz == num_cijk_entries);
315
316 // Print full 3-tensor to file
317 if (save_3tensor) {
318 for (int i=0; i<num_cijk_entries; ++i) {
319 cijk_file << ijk_triples[0][i] << ", "
320 << ijk_triples[1][i] << ", "
321 << ijk_triples[2][i] << ", "
322 << parts[i] << std::endl;
323 }
324 }
325 }
326
327 else if (partitioning_type == HG_MATRIX) {
328 // Build CRS graph from Cijk tensor
329 RCP<const EpetraExt::MultiComm> multiComm =
330 Stokhos::buildMultiComm(comm, basis_size, 1);
331 Stokhos::EpetraSparse3Tensor epetra_Cijk(basis, Cijk, multiComm);
332 RCP<const Epetra_CrsGraph> graph = epetra_Cijk.getStochasticGraph();
333 // RCP<const Epetra_CrsGraph> graph =
334 // Stokhos::sparse3Tensor2CrsGraph(*basis, *Cijk, comm);
335
336 // Partition graph using isorropia/zoltan's hypergraph partitioner
337 ParameterList params;
338 params.set("partitioning method", "hypergraph");
339 //int num_parts = basis_size / tile_size;
340 int num_parts = comm.NumProc();
341 std::stringstream ss;
342 ss << num_parts;
343 params.set<std::string>("num parts", ss.str());
344 std::stringstream ss2;
345 ss2 << imbalance_tol;
346 params.set<std::string>("imbalance tol", ss2.str());
347
348 /*
349 Isorropia::Epetra::Partitioner partitioner(graph, params);
350 partitioner.compute();
351
352 std::cout << "num parts requested = " << num_parts
353 << " num parts computed = " << partitioner.numProperties() << std::endl;
354
355 parts.resize(partitioner.numProperties());
356 int sz;
357 partitioner.extractPartsCopy(partitioner.numProperties(), sz, &parts[0]);
358 TEUCHOS_ASSERT(sz == partitioner.numProperties());
359
360 for (int i=0; i<sz; ++i)
361 std::cout << i << " " << parts[i] << std::endl;
362
363 // Create the permuted map
364 RCP<Epetra_Map> permuted_map = partitioner.createNewMap();
365 std::cout << *permuted_map << std::endl;
366 */
367 Teuchos::RCP<const Epetra_CrsGraph> rebalanced_graph =
368 Teuchos::rcp(Isorropia::Epetra::createBalancedCopy(
369 *graph, params));
370 Epetra_BlockMap permuted_map = rebalanced_graph->RowMap();
371 std::cout << permuted_map << std::endl;
372
373 //Epetra_CrsMatrix mat(Copy, *rebalanced_graph);
374 Epetra_CrsMatrix mat(Copy, *graph);
375 mat.FillComplete();
376 mat.PutScalar(1.0);
377 EpetraExt::RowMatrixToMatrixMarketFile(file.c_str(), mat);
378
379 // Print full 3-tensor to file
380 if (save_3tensor) {
381 Cijk_type::k_iterator k_begin = Cijk->k_begin();
382 Cijk_type::k_iterator k_end = Cijk->k_end();
383 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
384 int k = index(k_it);
385 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
386 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
387 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
388 int j = index(j_it);
389 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
390 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
391 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it){
392 int i = index(i_it);
393 cijk_file << i << ", " << j << ", " << k << ", " << parts[i]
394 << std::endl;
395 }
396 }
397 }
398 }
399 }
400
401 if (save_3tensor) {
402 cijk_file.close();
403 }
404
405 //Teuchos::TimeMonitor::summarize(std::cout);
406
407 }
408 catch (std::exception& e) {
409 std::cout << e.what() << std::endl;
410 }
411
412#ifdef HAVE_MPI
413 MPI_Finalize() ;
414#endif
415
416 return 0;
417}
Copy
BasisType
ProductBasisType
PartitioningType
@ HG_MATRIX
@ RCB
const PartitioningType partitioning_type_values[]
int main(int argc, char **argv)
OrderingType
@ LEXICOGRAPHIC_ORDERING
@ TOTAL_ORDERING
const int num_ordering_types
const BasisType basis_type_values[]
const char * ordering_type_names[]
const OrderingType ordering_type_values[]
const int num_basis_types
const int num_partitioning_types
const int num_growth_types
const char * basis_type_names[]
const int num_prod_basis_types
const char * prod_basis_type_names[]
const ProductBasisType prod_basis_type_values[]
const char * growth_type_names[]
const Stokhos::GrowthPolicy growth_type_values[]
@ JACOBI
@ CC_LEGENDRE
@ LEGENDRE
@ HERMITE
@ RYS
@ GP_LEGENDRE
@ COMPLETE
@ SMOLYAK
@ TOTAL
@ TENSOR
const char * partitioning_type_names[]
int FillComplete(bool OptimizeDataStorage=true)
int PutScalar(double ScalarConstant)
int NumProc() const
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis using Gauss-Patterson quadrature points.
Hermite polynomial basis.
Jacobi polynomial basis.
Legendre polynomial basis.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Rys polynomial basis.
Multivariate orthogonal polynomial basis generated from a Smolyak sparse grid.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
An isotropic total order index set.
A comparison functor implementing a strict weak ordering based total-order ordering,...
Teuchos::RCP< const EpetraExt::MultiComm > buildMultiComm(const Epetra_Comm &globalComm, int num_global_stochastic_blocks, int num_spatial_procs=-1)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.