Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
sparsity_example.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 
45 #ifdef HAVE_MPI
46 #include "Epetra_MpiComm.h"
47 #else
48 #include "Epetra_SerialComm.h"
49 #endif
50 
51 #include "Ifpack_RCMReordering.h"
52 
53 // sparsity_example
54 //
55 // usage:
56 // sparsity_example [options]
57 //
58 // output:
59 // prints the sparsity of the sparse 3 tensor specified by the basis,
60 // dimension, order, given by summing over the third index, to a matrix
61 // market file. This sparsity pattern yields the sparsity of the block
62 // stochastic Galerkin matrix which can be visualized, e.g., by matlab.
63 // The full/linear flag determines whether the third index ranges over
64 // the full polynomial dimension, or only over the zeroth and first order
65 // terms.
66 
67 // Basis types
69 const int num_basis_types = 6;
72 const char *basis_type_names[] = {
73  "hermite", "legendre", "clenshaw-curtis", "gauss-patterson", "rys", "jacobi" };
74 
75 // Growth policies
76 const int num_growth_types = 2;
79 const char *growth_type_names[] = { "slow", "moderate" };
80 
81 // Product Basis types
83 const int num_prod_basis_types = 4;
86 const char *prod_basis_type_names[] = {
87  "complete", "tensor", "total", "smolyak" };
88 
89 // Ordering types
91 const int num_ordering_types = 3;
94 const char *ordering_type_names[] = {
95  "total", "lexicographic", "morton-z" };
96 
97 int main(int argc, char **argv)
98 {
99  try {
100 
101  // Initialize MPI
102 #ifdef HAVE_MPI
103  MPI_Init(&argc,&argv);
104 #endif
105 
106  // Setup command line options
107  Teuchos::CommandLineProcessor CLP;
108  CLP.setDocString(
109  "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
110  int d = 3;
111  CLP.setOption("dimension", &d, "Stochastic dimension");
112  int p = 5;
113  CLP.setOption("order", &p, "Polynomial order");
114  double drop = 1.0e-12;
115  CLP.setOption("drop", &drop, "Drop tolerance");
116  std::string file = "A.mm";
117  CLP.setOption("filename", &file, "Matrix Market filename");
119  CLP.setOption("basis", &basis_type,
121  "Basis type");
123  CLP.setOption("growth", &growth_type,
125  "Growth type");
126  ProductBasisType prod_basis_type = COMPLETE;
127  CLP.setOption("product_basis", &prod_basis_type,
130  "Product basis type");
131  OrderingType ordering_type = TOTAL_ORDERING;
132  CLP.setOption("ordering", &ordering_type,
135  "Product basis ordering");
136  double alpha = 1.0;
137  CLP.setOption("alpha", &alpha, "Jacobi alpha index");
138  double beta = 1.0;
139  CLP.setOption("beta", &beta, "Jacobi beta index");
140  bool full = true;
141  CLP.setOption("full", "linear", &full, "Use full or linear expansion");
142  bool use_old = false;
143  CLP.setOption("old", "new", &use_old, "Use old or new Cijk algorithm");
144  bool print = false;
145  CLP.setOption("print", "no-print", &print, "Print Cijk to screen");
146  bool save_3tensor = false;
147  CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
148  "Save full 3tensor to file");
149  std::string file_3tensor = "Cijk.dat";
150  CLP.setOption("filename_3tensor", &file_3tensor,
151  "Filename to store full 3-tensor");
152  bool unique = false;
153  CLP.setOption("unique", "no-unique", &unique,
154  "Only save the unique non-zeros");
155  bool rcm = false;
156  CLP.setOption("rcm", "no-rcm", &rcm, "Reorder using RCM");
157 
158  // Parse arguments
159  CLP.parse( argc, argv );
160 
161  // Basis
162  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
163  for (int i=0; i<d; i++) {
164  if (basis_type == HERMITE)
165  bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(
166  p, true, growth_type));
167  else if (basis_type == LEGENDRE)
168  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(
169  p, true, growth_type));
170  else if (basis_type == CC_LEGENDRE)
171  bases[i] =
173  p, true));
174  else if (basis_type == GP_LEGENDRE)
175  bases[i] =
177  p, true));
178  else if (basis_type == RYS)
179  bases[i] = Teuchos::rcp(new Stokhos::RysBasis<int,double>(
180  p, 1.0, true, growth_type));
181  else if (basis_type == JACOBI)
182  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
183  p, alpha, beta, true, growth_type));
184  }
185  Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
189  if (prod_basis_type == COMPLETE)
190  basis =
192  bases, drop, use_old));
193  else if (prod_basis_type == TENSOR) {
194  if (ordering_type == TOTAL_ORDERING)
195  basis =
197  bases, drop));
198  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
199  basis =
201  bases, drop));
202  else if (ordering_type == MORTON_Z_ORDERING)
203  basis =
205  bases, drop));
206  }
207  else if (prod_basis_type == TOTAL) {
208  if (ordering_type == TOTAL_ORDERING)
209  basis =
211  bases, drop));
212  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
213  basis =
215  bases, drop));
216  else if (ordering_type == MORTON_Z_ORDERING)
217  basis =
219  bases, drop));
220  }
221  else if (prod_basis_type == SMOLYAK) {
222  Stokhos::TotalOrderIndexSet<int> index_set(d, p);
223  if (ordering_type == TOTAL_ORDERING)
224  basis =
226  bases, index_set, drop));
227  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
228  basis =
230  bases, index_set, drop));
231  else if (ordering_type == MORTON_Z_ORDERING)
232  basis =
234  bases, index_set, drop));
235  }
236 
237  // Triple product tensor
239  Teuchos::RCP<Cijk_type> Cijk;
240  if (full)
241  Cijk = basis->computeTripleProductTensor();
242  else
243  Cijk = basis->computeLinearTripleProductTensor();
244 
245  std::cout << "basis size = " << basis->size()
246  << " num nonzero Cijk entries = " << Cijk->num_entries()
247  << std::endl;
248 
249 #ifdef HAVE_MPI
250  Epetra_MpiComm comm(MPI_COMM_WORLD);
251 #else
252  Epetra_SerialComm comm;
253 #endif
254 
255  if (rcm) {
256  Teuchos::RCP<Cijk_type> Cijk3 = Teuchos::rcp(new Cijk_type);
257  {
258  Cijk_type::k_iterator k_begin = Cijk->k_begin();
259  Cijk_type::k_iterator k_end = Cijk->k_end();
260  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
261  int k = index(k_it);
262  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
263  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
264  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
265  int j = index(j_it);
266  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
267  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
268  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
269  int i = index(i_it);
270  double cijk = value(i_it);
271  if (i != 0 || (i == 0 && j == 0 && k == 0))
272  Cijk3->add_term(i, j, k, cijk);
273  }
274  }
275  }
276  }
277  Cijk3->fillComplete();
278 
279  Teuchos::RCP<Epetra_CrsGraph> graph =
280  Stokhos::sparse3Tensor2CrsGraph(*basis, *Cijk3, comm);
281  Epetra_CrsMatrix mat(Copy, *graph);
282  mat.FillComplete();
283  mat.PutScalar(1.0);
284  Ifpack_RCMReordering ifpack_rcm;
285  ifpack_rcm.SetParameter("reorder: root node", basis->size()-1);
286  ifpack_rcm.Compute(mat);
287 
288  Teuchos::RCP<Cijk_type> Cijk2 = Teuchos::rcp(new Cijk_type);
289  Cijk_type::k_iterator k_begin = Cijk->k_begin();
290  Cijk_type::k_iterator k_end = Cijk->k_end();
291  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
292  int k = ifpack_rcm.Reorder(index(k_it));
293  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
294  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
295  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
296  int j = ifpack_rcm.Reorder(index(j_it));
297  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
298  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
299  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
300  int i = ifpack_rcm.Reorder(index(i_it));
301  double cijk = value(i_it);
302  Cijk2->add_term(i, j, k, cijk);
303  }
304  }
305  }
306  Cijk2->fillComplete();
307  Cijk = Cijk2;
308  }
309 
310  if (print) {
311  std::cout << *Cijk << std::endl;
312  }
313 
314  // Print triple product sparsity to matrix market file
315  Stokhos::sparse3Tensor2MatrixMarket(*basis, *Cijk, comm, file);
316 
317  // Print full 3-tensor to file
318  if (save_3tensor) {
319  std::ofstream cijk_file(file_3tensor.c_str());
320  cijk_file.precision(14);
321  cijk_file.setf(std::ios::scientific);
322  cijk_file << "i, j, k, cijk" << std::endl;
323  Cijk_type::k_iterator k_begin = Cijk->k_begin();
324  Cijk_type::k_iterator k_end = Cijk->k_end();
325  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
326  int k = index(k_it);
327  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
328  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
329  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
330  int j = index(j_it);
331  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
332  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
333  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
334  int i = index(i_it);
335  double cijk = value(i_it);
336  if (!unique || ( i >= j && j >= k ))
337  cijk_file << i << ", "
338  << j << ", "
339  << k << ", "
340  << cijk << std::endl;
341  }
342  }
343  }
344  cijk_file.close();
345  }
346 
347  Teuchos::TimeMonitor::summarize(std::cout);
348 
349  }
350  catch (std::exception& e) {
351  std::cout << e.what() << std::endl;
352  }
353 
354  return 0;
355 }
Hermite polynomial basis.
SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
const ProductBasisType prod_basis_type_values[]
const char * basis_type_names[]
void sparse3Tensor2MatrixMarket(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm, const std::string &file)
const char * ordering_type_names[]
const char * growth_type_names[]
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
const BasisType basis_type_values[]
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const int num_growth_types
OrderingType
const int num_prod_basis_types
Legendre polynomial basis using Gauss-Patterson quadrature points.
const int num_ordering_types
ProductBasisType
Rys polynomial basis.
Jacobi polynomial basis.
BasisType
ordinal_type num_entries() const
Return number of non-zero entries.
expr expr expr expr j
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Multivariate orthogonal polynomial basis generated from a Smolyak sparse grid.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
int main(int argc, char **argv)
An isotropic total order index set.
A comparison functor implementing a strict weak ordering based Morton Z-ordering. ...
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
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.
SparseArrayIterator< index_iterator, value_iterator >::value_reference value(const SparseArrayIterator< index_iterator, value_iterator > &it)
const char * prod_basis_type_names[]
A comparison functor implementing a strict weak ordering based lexographic ordering.
const Stokhos::GrowthPolicy growth_type_values[]
const OrderingType ordering_type_values[]
const int num_basis_types