Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cijk_nonzeros.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
44#include "Stokhos_Epetra.hpp"
45#include "Teuchos_CommandLineProcessor.hpp"
46
47#ifdef HAVE_MPI
48#include "Epetra_MpiComm.h"
49#else
50#include "Epetra_SerialComm.h"
51#endif
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
69const int num_basis_types = 6;
72const char *basis_type_names[] = {
73 "hermite", "legendre", "clenshaw-curtis", "gauss-patterson", "rys", "jacobi" };
74
75// Growth policies
76const int num_growth_types = 2;
79const char *growth_type_names[] = { "slow", "moderate" };
80
81// Product Basis types
86const char *prod_basis_type_names[] = {
87 "complete", "tensor", "total", "smolyak" };
88
89using Teuchos::Array;
90using Teuchos::RCP;
91using Teuchos::rcp;
92
94 int i, total_nz;
95 Array< Array<int> > nz_tiles;
96};
97
98struct NZCompare {
99 bool operator() (const CijkNonzeros& a,
100 const CijkNonzeros& b) const {
101 return (a.total_nz == b.total_nz) ? (a.i < b.i) : (a.total_nz > b.total_nz) ;
102 }
103};
104
106 bool operator() (const std::pair<int,int>& a,
107 const std::pair<int,int>& b) const {
108 return (a.second == b.second) ? (a.first < b.first) : (a.second > b.second);
109 }
110};
111
112int main(int argc, char **argv)
113{
114 try {
115
116 // Initialize MPI
117#ifdef HAVE_MPI
118 MPI_Init(&argc,&argv);
119#endif
120
121 // Setup command line options
122 Teuchos::CommandLineProcessor CLP;
123 CLP.setDocString(
124 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
125 int d = 3;
126 CLP.setOption("dimension", &d, "Stochastic dimension");
127 int p = 5;
128 CLP.setOption("order", &p, "Polynomial order");
129 double drop = 1.0e-12;
130 CLP.setOption("drop", &drop, "Drop tolerance");
131 std::string file = "A.mm";
132 CLP.setOption("filename", &file, "Matrix Market filename");
134 CLP.setOption("basis", &basis_type,
136 "Basis type");
138 CLP.setOption("growth", &growth_type,
140 "Growth type");
141 ProductBasisType prod_basis_type = COMPLETE;
142 CLP.setOption("product_basis", &prod_basis_type,
145 "Product basis type");
146 double alpha = 1.0;
147 CLP.setOption("alpha", &alpha, "Jacobi alpha index");
148 double beta = 1.0;
149 CLP.setOption("beta", &beta, "Jacobi beta index");
150 bool full = true;
151 CLP.setOption("full", "linear", &full, "Use full or linear expansion");
152 bool use_old = false;
153 CLP.setOption("old", "new", &use_old, "Use old or new Cijk algorithm");
154 int tile_size = 100;
155 CLP.setOption("tile_size", &tile_size, "Tile size");
156
157 // Parse arguments
158 CLP.parse( argc, argv );
159
160 // Basis
161 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
162 for (int i=0; i<d; i++) {
163 if (basis_type == HERMITE)
164 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(
165 p, true, growth_type));
166 else if (basis_type == LEGENDRE)
167 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(
168 p, true, growth_type));
169 else if (basis_type == CC_LEGENDRE)
170 bases[i] =
172 p, true));
173 else if (basis_type == GP_LEGENDRE)
174 bases[i] =
176 p, true));
177 else if (basis_type == RYS)
178 bases[i] = Teuchos::rcp(new Stokhos::RysBasis<int,double>(
179 p, 1.0, true, growth_type));
180 else if (basis_type == JACOBI)
181 bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
182 p, alpha, beta, true, growth_type));
183 }
184 RCP<const Stokhos::ProductBasis<int,double> > basis;
185 if (prod_basis_type == COMPLETE)
186 basis =
188 bases, drop, use_old));
189 else if (prod_basis_type == TENSOR)
190 basis =
192 bases, drop));
193 else if (prod_basis_type == TOTAL)
194 basis =
196 bases, drop));
197 else if (prod_basis_type == SMOLYAK) {
198 Stokhos::TotalOrderIndexSet<int> index_set(d, p);
199 basis =
200 Teuchos::rcp(new Stokhos::SmolyakBasis<int,double>(
201 bases, index_set, drop));
202 }
203
204 // Triple product tensor
205 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
206 RCP<Cijk_type> Cijk;
207 if (full)
208 Cijk = basis->computeTripleProductTensor();
209 else
210 Cijk = basis->computeLinearTripleProductTensor();
211
212 int sz = basis->size();
213 std::cout << "basis size = " << sz
214 << " num nonzero Cijk entries = " << Cijk->num_entries()
215 << std::endl;
216
217 // Setup tiles
218 if (tile_size > sz)
219 tile_size = sz;
220 int j_sz = sz;
221 int k_sz = sz;
222 if (!full)
223 k_sz = basis->dimension()+1;
224 int nj_tiles = j_sz / tile_size;
225 int nk_tiles = k_sz / tile_size;
226 if (j_sz - nj_tiles*tile_size > 0)
227 ++nj_tiles;
228 if (k_sz - nk_tiles*tile_size > 0)
229 ++nk_tiles;
230 Array<CijkNonzeros> nz(sz);
231 for (int i=0; i<sz; ++i) {
232 nz[i].i = i;
233 nz[i].nz_tiles.resize(nj_tiles);
234 for (int j=0; j<nj_tiles; ++j)
235 nz[i].nz_tiles[j].resize(nk_tiles);
236 }
237
238 // Get number of nonzeros in Cijk for each i
239 Cijk_type::k_iterator k_begin = Cijk->k_begin();
240 Cijk_type::k_iterator k_end = Cijk->k_end();
241 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
242 int k = index(k_it);
243 int k_tile = k / tile_size;
244 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
245 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
246 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
247 int j = index(j_it);
248 int j_tile = j / tile_size;
249 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
250 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
251 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
252 int i = index(i_it);
253 ++nz[i].total_nz;
254 ++nz[i].nz_tiles[j_tile][k_tile];
255 }
256 }
257 }
258
259 // Sort based on total number of nonzeros
260 std::sort(nz.begin(), nz.end(), NZCompare());
261
262 // Print nonzeros
263 int w_index = 3;
264 int w_nz = 5;
265 int w_tile = 4;
266 for (int i=0; i<nz.size(); ++i) {
267 int idx = nz[i].i;
268 std::cout << std::setw(w_index) << idx << " "
269 << basis->term(idx) << ": "
270 << std::setw(w_nz) << nz[i].total_nz
271 << ", ";
272 for (int j=0; j<nj_tiles; ++j)
273 for (int k=0; k<nk_tiles; ++k)
274 std::cout << std::setw(w_tile) << nz[i].nz_tiles[j][k] << " ";
275 std::cout << std::endl;
276 }
277
278 // Add up the nonzeros for each (j,k) tile
279 Array< Array<int> > total_nz_tiles(nj_tiles);
280 int total_nz = 0;
281 for (int j=0; j<nj_tiles; ++j)
282 total_nz_tiles[j].resize(nk_tiles);
283 for (int i=0; i<nz.size(); ++i) {
284 total_nz += nz[i].total_nz;
285 for (int j=0; j<nj_tiles; ++j)
286 for (int k=0; k<nk_tiles; ++k)
287 total_nz_tiles[j][k] += nz[i].nz_tiles[j][k];
288 }
289 int w_total = (w_index+1) + (2*basis->dimension()+5) + w_nz;
290 std::cout << std::endl << std::setw(w_total) << total_nz << ", ";
291 for (int j=0; j<nj_tiles; ++j)
292 for (int k=0; k<nk_tiles; ++k)
293 std::cout << std::setw(w_tile) << total_nz_tiles[j][k] << " ";
294 std::cout << std::endl;
295
296 // Now partition Cijk for each tile
297 Array< Array< RCP<Cijk_type> > > Cijk_tile(nj_tiles);
298 for (int j=0; j<nj_tiles; ++j) {
299 Cijk_tile[j].resize(nk_tiles);
300 for (int k=0; k<nk_tiles; ++k)
301 Cijk_tile[j][k] = rcp(new Cijk_type);
302 }
303 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
304 int k = index(k_it);
305 int k_tile = k / tile_size;
306 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
307 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
308 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
309 int j = index(j_it);
310 int j_tile = j / tile_size;
311 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
312 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
313 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
314 int i = index(i_it);
315 double c = value(i_it);
316 Cijk_tile[j_tile][k_tile]->add_term(i,j,k,c);
317 }
318 }
319 }
320 for (int j=0; j<nj_tiles; ++j)
321 for (int k=0; k<nk_tiles; ++k)
322 Cijk_tile[j][k]->fillComplete();
323
324
325 Array< Array< std::map<int,int> > > nz_tile(nj_tiles);
326 Array< Array< Array< std::pair<int,int> > > > sorted_nz_tile(nj_tiles);
327 for (int j_tile=0; j_tile<nj_tiles; ++j_tile) {
328 nz_tile[j_tile].resize(nk_tiles);
329 sorted_nz_tile[j_tile].resize(nk_tiles);
330 for (int k_tile=0; k_tile<nk_tiles; ++k_tile) {
331
332 // Count nonzeros for each i, for each tile
333 Cijk_type::k_iterator k_begin = Cijk_tile[j_tile][k_tile]->k_begin();
334 Cijk_type::k_iterator k_end = Cijk_tile[j_tile][k_tile]->k_end();
335 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
336 //int k = index(k_it);
337 Cijk_type::kj_iterator j_begin =
338 Cijk_tile[j_tile][k_tile]->j_begin(k_it);
339 Cijk_type::kj_iterator j_end =
340 Cijk_tile[j_tile][k_tile]->j_end(k_it);
341 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
342 //int j = index(j_it);
343 Cijk_type::kji_iterator i_begin =
344 Cijk_tile[j_tile][k_tile]->i_begin(j_it);
345 Cijk_type::kji_iterator i_end =
346 Cijk_tile[j_tile][k_tile]->i_end(j_it);
347 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it){
348 int i = index(i_it);
349 if (nz_tile[j_tile][k_tile].count(i) == 0)
350 nz_tile[j_tile][k_tile][i] = 1;
351 else
352 ++(nz_tile[j_tile][k_tile][i]);
353 }
354 }
355 }
356
357 // Sort based on non-zeros for each i, for each tile
358 sorted_nz_tile[j_tile][k_tile].resize(nz_tile[j_tile][k_tile].size());
359 int idx=0;
360 for (std::map<int,int>::iterator it = nz_tile[j_tile][k_tile].begin();
361 it != nz_tile[j_tile][k_tile].end(); ++it) {
362 sorted_nz_tile[j_tile][k_tile][idx] =
363 std::make_pair(it->first, it->second);
364 ++idx;
365 }
366 std::sort( sorted_nz_tile[j_tile][k_tile].begin(),
367 sorted_nz_tile[j_tile][k_tile].end(),
368 NZPairCompare() );
369
370 // Print number of non-zeros for each i, for each tile
371 std::cout << std::endl
372 << "Tile (" << j_tile << ", " << k_tile << "):" << std::endl;
373 for (int i=0; i<sorted_nz_tile[j_tile][k_tile].size(); ++i) {
374 int idx = sorted_nz_tile[j_tile][k_tile][i].first;
375 std::cout << std::setw(w_index) << idx << " "
376 << basis->term(idx) << ": "
377 << std::setw(w_nz) << sorted_nz_tile[j_tile][k_tile][i].second
378 << std::endl;
379 if (i % 32 == 31)
380 std::cout << std::endl;
381 }
382 }
383 }
384
385 }
386 catch (std::exception& e) {
387 std::cout << e.what() << std::endl;
388 }
389
390 return 0;
391}
int main(int argc, char **argv)
const BasisType basis_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[]
BasisType
@ JACOBI
@ CC_LEGENDRE
@ LEGENDRE
@ HERMITE
@ RYS
@ GP_LEGENDRE
ProductBasisType
@ COMPLETE
@ SMOLYAK
@ TOTAL
@ TENSOR
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.
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.
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
Array< Array< int > > nz_tiles
bool operator()(const CijkNonzeros &a, const CijkNonzeros &b) const
bool operator()(const std::pair< int, int > &a, const std::pair< int, int > &b) const