Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
sparsity_slice_example.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#include <sstream>
48
49#ifdef HAVE_MPI
50#include "Epetra_MpiComm.h"
51#else
52#include "Epetra_SerialComm.h"
53#endif
54#include "Epetra_LocalMap.h"
55#include "Epetra_CrsMatrix.h"
56#include "EpetraExt_RowMatrixOut.h"
57
58// sparsity_example
59//
60// usage:
61// sparsity_slice_example [options]
62//
63// output:
64// prints the sparsity of the sparse 3 tensor specified by the basis,
65// dimension, and order, printing a matrix-market file for each k-slice.
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 = 3;
74const char *basis_type_names[] = { "hermite", "legendre", "rys" };
75
76int main(int argc, char **argv)
77{
78 int num_k = 0;
79
80 try {
81
82 // Initialize MPI
83#ifdef HAVE_MPI
84 MPI_Init(&argc,&argv);
85#endif
86
87 // Setup command line options
88 Teuchos::CommandLineProcessor CLP;
89 CLP.setDocString(
90 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
91 int d = 3;
92 CLP.setOption("dimension", &d, "Stochastic dimension");
93 int p = 5;
94 CLP.setOption("order", &p, "Polynomial order");
95 double drop = 1.0e-15;
96 CLP.setOption("drop", &drop, "Drop tolerance");
97 std::string file_base = "A";
98 CLP.setOption("base filename", &file_base, "Base filename for matrix market files");
100 CLP.setOption("basis", &basis_type,
102 "Basis type");
103 bool full = true;
104 CLP.setOption("full", "linear", &full, "Use full or linear expansion");
105 bool use_old = false;
106 CLP.setOption("old", "new", &use_old, "Use old or new Cijk algorithm");
107
108 // Parse arguments
109 CLP.parse( argc, argv );
110
111 // Basis
112 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
113 for (int i=0; i<d; i++) {
114 if (basis_type == HERMITE)
115 bases[i] = Teuchos::rcp(new Stokhos::HermiteBasis<int,double>(p));
116 else if (basis_type == LEGENDRE)
117 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
118 else if (basis_type == RYS)
119 bases[i] = Teuchos::rcp(new Stokhos::RysBasis<int,double>(p, 1.0,
120 false));
121 }
122 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
124 drop,
125 use_old));
126
127 // Triple product tensor
128 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
129 if (full) {
130 num_k = basis->size();
131 Cijk = basis->computeTripleProductTensor();
132 }
133 else {
134 num_k = basis->dimension()+1;
135 Cijk = basis->computeLinearTripleProductTensor();
136 }
137
138 std::cout << "basis size = " << basis->size()
139 << " num nonzero Cijk entries = " << Cijk->num_entries()
140 << std::endl;
141
142#ifdef HAVE_MPI
143 Epetra_MpiComm comm(MPI_COMM_WORLD);
144#else
146#endif
147
148 // Number of stochastic rows
149 int num_rows = basis->size();
150
151 // Replicated local map
152 Epetra_LocalMap map(num_rows, 0, comm);
153
154 // Loop over Cijk entries including a non-zero in the graph at
155 // indices (i,j) if Cijk is non-zero for each k
156 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
157 double one = 1.0;
158 for (Cijk_type::k_iterator k_it=Cijk->k_begin();
159 k_it!=Cijk->k_end(); ++k_it) {
160 int k = index(k_it);
161 Epetra_CrsMatrix mat(Copy, map, 1);
162 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
163 j_it != Cijk->j_end(k_it); ++j_it) {
164 int j = index(j_it);
165 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
166 i_it != Cijk->i_end(j_it); ++i_it) {
167 int i = index(i_it);
168 mat.InsertGlobalValues(i, 1, &one, &j);
169 }
170 }
171 mat.FillComplete();
172
173 // Construct file name
174 std::stringstream ss;
175 ss << file_base << "_" << k << ".mm";
176 std::string file = ss.str();
177
178 // Save matrix to file
179 EpetraExt::RowMatrixToMatrixMarketFile(file.c_str(), mat);
180 }
181
182 }
183 catch (std::exception& e) {
184 std::cout << e.what() << std::endl;
185 }
186
187 return num_k;
188}
Copy
BasisType
int FillComplete(bool OptimizeDataStorage=true)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Hermite polynomial basis.
Legendre polynomial basis.
Rys polynomial basis.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
int main(int argc, char **argv)
const BasisType basis_type_values[]
const int num_basis_types
const char * basis_type_names[]