Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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#include <fstream>
54#include <iostream>
55
56// sparsity_example
57//
58// usage:
59// sparsity_example [options]
60//
61// output:
62// prints the sparsity of the sparse 3 tensor specified by the basis,
63// dimension, order, given by summing over the third index, to a matrix
64// market file. This sparsity pattern yields the sparsity of the block
65// stochastic Galerkin matrix which can be visualized, e.g., by matlab.
66// The full/linear flag determines whether the third index ranges over
67// the full polynomial dimension, or only over the zeroth and first order
68// terms.
69
70// Basis types
72const int num_basis_types = 6;
75const char *basis_type_names[] = {
76 "hermite", "legendre", "clenshaw-curtis", "gauss-patterson", "rys", "jacobi" };
77
78// Growth policies
79const int num_growth_types = 2;
82const char *growth_type_names[] = { "slow", "moderate" };
83
84// Product Basis types
89const char *prod_basis_type_names[] = {
90 "complete", "tensor", "total", "smolyak" };
91
92// Ordering types
94const int num_ordering_types = 3;
97const char *ordering_type_names[] = {
98 "total", "lexicographic", "morton-z" };
99
100int main(int argc, char **argv)
101{
102 try {
103
104 // Initialize MPI
105#ifdef HAVE_MPI
106 MPI_Init(&argc,&argv);
107#endif
108
109 // Setup command line options
110 Teuchos::CommandLineProcessor CLP;
111 CLP.setDocString(
112 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
113 int d = 3;
114 CLP.setOption("dimension", &d, "Stochastic dimension");
115 int p = 5;
116 CLP.setOption("order", &p, "Polynomial order");
117 double drop = 1.0e-12;
118 CLP.setOption("drop", &drop, "Drop tolerance");
119 std::string file = "A.mm";
120 CLP.setOption("filename", &file, "Matrix Market filename");
122 CLP.setOption("basis", &basis_type,
124 "Basis type");
126 CLP.setOption("growth", &growth_type,
128 "Growth type");
129 ProductBasisType prod_basis_type = COMPLETE;
130 CLP.setOption("product_basis", &prod_basis_type,
133 "Product basis type");
134 OrderingType ordering_type = TOTAL_ORDERING;
135 CLP.setOption("ordering", &ordering_type,
138 "Product basis ordering");
139 double alpha = 1.0;
140 CLP.setOption("alpha", &alpha, "Jacobi alpha index");
141 double beta = 1.0;
142 CLP.setOption("beta", &beta, "Jacobi beta index");
143 bool full = true;
144 CLP.setOption("full", "linear", &full, "Use full or linear expansion");
145 bool use_old = false;
146 CLP.setOption("old", "new", &use_old, "Use old or new Cijk algorithm");
147 bool print = false;
148 CLP.setOption("print", "no-print", &print, "Print Cijk to screen");
149 bool save_3tensor = false;
150 CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
151 "Save full 3tensor to file");
152 std::string file_3tensor = "Cijk.dat";
153 CLP.setOption("filename_3tensor", &file_3tensor,
154 "Filename to store full 3-tensor");
155 bool unique = false;
156 CLP.setOption("unique", "no-unique", &unique,
157 "Only save the unique non-zeros");
158 bool rcm = false;
159 CLP.setOption("rcm", "no-rcm", &rcm, "Reorder using RCM");
160
161 // Parse arguments
162 CLP.parse( argc, argv );
163
164 // Basis
165 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
166 for (int i=0; i<d; i++) {
167 if (basis_type == HERMITE)
168 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(
169 p, true, growth_type));
170 else if (basis_type == LEGENDRE)
171 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(
172 p, true, growth_type));
173 else if (basis_type == CC_LEGENDRE)
174 bases[i] =
176 p, true));
177 else if (basis_type == GP_LEGENDRE)
178 bases[i] =
180 p, true));
181 else if (basis_type == RYS)
182 bases[i] = Teuchos::rcp(new Stokhos::RysBasis<int,double>(
183 p, 1.0, true, growth_type));
184 else if (basis_type == JACOBI)
185 bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
186 p, alpha, beta, true, growth_type));
187 }
188 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
192 if (prod_basis_type == COMPLETE)
193 basis =
195 bases, drop, use_old));
196 else if (prod_basis_type == TENSOR) {
197 if (ordering_type == TOTAL_ORDERING)
198 basis =
200 bases, drop));
201 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
202 basis =
204 bases, drop));
205 else if (ordering_type == MORTON_Z_ORDERING)
206 basis =
208 bases, drop));
209 }
210 else if (prod_basis_type == TOTAL) {
211 if (ordering_type == TOTAL_ORDERING)
212 basis =
214 bases, drop));
215 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
216 basis =
218 bases, drop));
219 else if (ordering_type == MORTON_Z_ORDERING)
220 basis =
222 bases, drop));
223 }
224 else if (prod_basis_type == SMOLYAK) {
225 Stokhos::TotalOrderIndexSet<int> index_set(d, p);
226 if (ordering_type == TOTAL_ORDERING)
227 basis =
229 bases, index_set, drop));
230 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
231 basis =
233 bases, index_set, drop));
234 else if (ordering_type == MORTON_Z_ORDERING)
235 basis =
237 bases, index_set, drop));
238 }
239
240 // Triple product tensor
241 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
242 Teuchos::RCP<Cijk_type> Cijk;
243 if (full)
244 Cijk = basis->computeTripleProductTensor();
245 else
246 Cijk = basis->computeLinearTripleProductTensor();
247
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 if (rcm) {
259 Teuchos::RCP<Cijk_type> Cijk3 = Teuchos::rcp(new Cijk_type);
260 {
261 Cijk_type::k_iterator k_begin = Cijk->k_begin();
262 Cijk_type::k_iterator k_end = Cijk->k_end();
263 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
264 int k = index(k_it);
265 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
266 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
267 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
268 int j = index(j_it);
269 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
270 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
271 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
272 int i = index(i_it);
273 double cijk = value(i_it);
274 if (i != 0 || (i == 0 && j == 0 && k == 0))
275 Cijk3->add_term(i, j, k, cijk);
276 }
277 }
278 }
279 }
280 Cijk3->fillComplete();
281
282 Teuchos::RCP<Epetra_CrsGraph> graph =
283 Stokhos::sparse3Tensor2CrsGraph(*basis, *Cijk3, comm);
284 Epetra_CrsMatrix mat(Copy, *graph);
285 mat.FillComplete();
286 mat.PutScalar(1.0);
287 Ifpack_RCMReordering ifpack_rcm;
288 ifpack_rcm.SetParameter("reorder: root node", basis->size()-1);
289 ifpack_rcm.Compute(mat);
290
291 Teuchos::RCP<Cijk_type> Cijk2 = Teuchos::rcp(new Cijk_type);
292 Cijk_type::k_iterator k_begin = Cijk->k_begin();
293 Cijk_type::k_iterator k_end = Cijk->k_end();
294 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
295 int k = ifpack_rcm.Reorder(index(k_it));
296 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
297 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
298 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
299 int j = ifpack_rcm.Reorder(index(j_it));
300 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
301 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
302 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
303 int i = ifpack_rcm.Reorder(index(i_it));
304 double cijk = value(i_it);
305 Cijk2->add_term(i, j, k, cijk);
306 }
307 }
308 }
309 Cijk2->fillComplete();
310 Cijk = Cijk2;
311 }
312
313 if (print) {
314 std::cout << *Cijk << std::endl;
315 }
316
317 // Print triple product sparsity to matrix market file
318 Stokhos::sparse3Tensor2MatrixMarket(*basis, *Cijk, comm, file);
319
320 // Print full 3-tensor to file
321 if (save_3tensor) {
322 std::ofstream cijk_file(file_3tensor.c_str());
323 cijk_file.precision(14);
324 cijk_file.setf(std::ios::scientific);
325 cijk_file << "i, j, k, cijk" << std::endl;
326 Cijk_type::k_iterator k_begin = Cijk->k_begin();
327 Cijk_type::k_iterator k_end = Cijk->k_end();
328 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
329 int k = index(k_it);
330 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
331 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
332 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
333 int j = index(j_it);
334 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
335 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
336 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
337 int i = index(i_it);
338 double cijk = value(i_it);
339 if (!unique || ( i >= j && j >= k ))
340 cijk_file << i << ", "
341 << j << ", "
342 << k << ", "
343 << cijk << std::endl;
344 }
345 }
346 }
347 cijk_file.close();
348 }
349
350 Teuchos::TimeMonitor::summarize(std::cout);
351
352 }
353 catch (std::exception& e) {
354 std::cout << e.what() << std::endl;
355 }
356
357 return 0;
358}
Copy
BasisType
ProductBasisType
OrderingType
int FillComplete(bool OptimizeDataStorage=true)
int PutScalar(double ScalarConstant)
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.
A comparison functor implementing a strict weak ordering based Morton Z-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.
ordinal_type num_entries() const
Return number of non-zero entries.
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< 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.
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)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
int main(int argc, char **argv)
@ MORTON_Z_ORDERING
@ 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_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[]
@ CC_LEGENDRE
@ LEGENDRE
@ HERMITE
@ GP_LEGENDRE
@ COMPLETE
@ SMOLYAK