Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cijk_ltb_partition_level.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.hpp"
44
45#include "Teuchos_CommandLineProcessor.hpp"
46#include "Teuchos_Array.hpp"
47#include "Teuchos_RCP.hpp"
48
49#include <fstream>
50#include <iostream>
51
52using Teuchos::Array;
53using Teuchos::RCP;
54using Teuchos::rcp;
55
56int main(int argc, char **argv)
57{
58 try {
59
60 // Setup command line options
61 Teuchos::CommandLineProcessor CLP;
62 CLP.setDocString(
63 "This example generates partitions the Cijk tensor for a lexicographic tree basis.\n");
64 int d = 3;
65 CLP.setOption("dimension", &d, "Stochastic dimension");
66 int p = 5;
67 CLP.setOption("order", &p, "Polynomial order");
68 double drop = 1.0e-12;
69 CLP.setOption("drop", &drop, "Drop tolerance");
70 bool symmetric = true;
71 CLP.setOption("symmetric", "asymmetric", &symmetric, "Use basis polynomials with symmetric PDF");
72 int level = 1;
73 CLP.setOption("level", &level, "Level to partition");
74 bool save_3tensor = false;
75 CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
76 "Save full 3tensor to file");
77 std::string file_3tensor = "Cijk.dat";
78 CLP.setOption("filename_3tensor", &file_3tensor,
79 "Filename to store full 3-tensor");
80
81 // Parse arguments
82 CLP.parse( argc, argv );
83
84 // Basis
85 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
86 const double alpha = 1.0;
87 const double beta = symmetric ? 1.0 : 2.0 ;
88 for (int i=0; i<d; i++) {
89 bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
90 p, alpha, beta, true));
91 }
94 RCP<const basis_type> basis = Teuchos::rcp(new basis_type(bases, drop));
95
96 // Build LTB Cijk
97 typedef Stokhos::LTBSparse3Tensor<int,double> Cijk_LTB_type;
98 typedef Cijk_LTB_type::CijkNode node_type;
99 Teuchos::RCP<Cijk_LTB_type> Cijk =
100 computeTripleProductTensorLTB(*basis, symmetric);
101
102 int sz = basis->size();
103 std::cout << "basis size = " << sz
104 << " num nonzero Cijk entries = " << Cijk->num_entries()
105 << std::endl;
106
107 // Setup partitions
108 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
109 Teuchos::Array< int > index_stack;
110 node_stack.push_back(Cijk->getHeadNode());
111 index_stack.push_back(0);
112 Teuchos::RCP<const node_type> node;
113 int child_index;
114 Teuchos::Array< Teuchos::RCP<const node_type> > partition_stack;
115 int my_level = 0;
116 while (node_stack.size() > 0) {
117 node = node_stack.back();
118 child_index = index_stack.back();
119
120 // Leaf -- If we got here, just push this node into the partitions
121 if (node->is_leaf) {
122 partition_stack.push_back(node);
123 node_stack.pop_back();
124 index_stack.pop_back();
125 --my_level;
126 }
127
128 // Put nodes into partition if level matches
129 else if (my_level == level) {
130 partition_stack.push_back(node);
131 node_stack.pop_back();
132 index_stack.pop_back();
133 --my_level;
134 }
135
136 // More children to process -- process them first
137 else if (child_index < node->children.size()) {
138 ++index_stack.back();
139 node = node->children[child_index];
140 node_stack.push_back(node);
141 index_stack.push_back(0);
142 ++my_level;
143 }
144
145 // No more children
146 else {
147 node_stack.pop_back();
148 index_stack.pop_back();
149 --my_level;
150 }
151
152 }
153
154 // Print statistics
155 int max_i_size = 0, max_j_size = 0, max_k_size = 0;
156 for (int part=0; part<partition_stack.size(); ++part) {
157 node = partition_stack[part];
158 if (node->i_size > max_i_size) max_i_size = node->i_size;
159 if (node->j_size > max_j_size) max_j_size = node->j_size;
160 if (node->k_size > max_k_size) max_k_size = node->k_size;
161 }
162 std::cout << "num partitions = " << partition_stack.size() << std::endl
163 << "max i size = " << max_i_size << std::endl
164 << "max j size = " << max_j_size << std::endl
165 << "max k size = " << max_k_size << std::endl;
166
167 // Build flat list of (i,j,k,part) tuples
169 Teuchos::Array< Teuchos::Array<int> > tuples;
170 for (int part=0; part<partition_stack.size(); ++part) {
171 node = partition_stack[part];
172 node_stack.push_back(node);
173 index_stack.push_back(0);
174 while (node_stack.size() > 0) {
175 node = node_stack.back();
176 child_index = index_stack.back();
177
178 // Leaf -- store (i,j,k,part) tuples
179 if (node->is_leaf) {
180 Cijk_Iterator cijk_iterator(node->p_i,
181 node->p_j,
182 node->p_k,
183 symmetric);
184 bool more = true;
185 while (more) {
186 Teuchos::Array<int> t(4);
187 int I = node->i_begin + cijk_iterator.i;
188 int J = node->j_begin + cijk_iterator.j;
189 int K = node->k_begin + cijk_iterator.k;
190 t[0] = I;
191 t[1] = J;
192 t[2] = K;
193 t[3] = part;
194 tuples.push_back(t);
195 more = cijk_iterator.increment();
196 }
197 node_stack.pop_back();
198 index_stack.pop_back();
199 }
200
201 // More children to process -- process them first
202 else if (child_index < node->children.size()) {
203 ++index_stack.back();
204 node = node->children[child_index];
205 node_stack.push_back(node);
206 index_stack.push_back(0);
207 }
208
209 // No more children
210 else {
211 node_stack.pop_back();
212 index_stack.pop_back();
213 }
214
215 }
216 }
217
218 // Print full 3-tensor to file
219 if (save_3tensor) {
220 std::ofstream cijk_file(file_3tensor.c_str());
221 cijk_file.precision(14);
222 cijk_file.setf(std::ios::scientific);
223 cijk_file << "i, j, k, part" << std::endl;
224 for (int i=0; i<tuples.size(); ++i) {
225 cijk_file << tuples[i][0] << ", "
226 << tuples[i][1] << ", "
227 << tuples[i][2] << ", "
228 << tuples[i][3] << std::endl;
229 }
230 cijk_file.close();
231 }
232
233 }
234 catch (std::exception& e) {
235 std::cout << e.what() << std::endl;
236 }
237
238 return 0;
239}
Kokkos::Serial node_type
int main(int argc, char **argv)
Jacobi polynomial basis.
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
A comparison functor implementing a strict weak ordering based lexographic ordering.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Stokhos::LegendreBasis< int, double > basis_type