Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
TestEpetraMatrixFreeApply.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 <string>
43#include <iostream>
44#include <cstdlib>
45
46#include "Epetra_config.h"
47#ifdef HAVE_MPI
48# include "Epetra_MpiComm.h"
49#else
50# include "Epetra_SerialComm.h"
51#endif
52
53// Stokhos Stochastic Galerkin
54#include "Stokhos_Epetra.hpp"
55#include "EpetraExt_BlockUtility.h"
56
57#include "Teuchos_VerboseObject.hpp"
58#include "Teuchos_GlobalMPISession.hpp"
59#include "Teuchos_CommandLineProcessor.hpp"
60#include "Teuchos_StandardCatchMacros.hpp"
61
62template< typename IntType >
63inline
64IntType map_fem_graph_coord( const IntType & N ,
65 const IntType & i ,
66 const IntType & j ,
67 const IntType & k )
68{
69 return k + N * ( j + N * i );
70}
71
72template < typename ordinal >
73inline
74ordinal generate_fem_graph( ordinal N ,
75 Teuchos::Array< Teuchos::Array<ordinal> > & graph )
76{
77 graph.resize( N * N * N , Teuchos::Array<ordinal>() );
78
79 ordinal total = 0 ;
80
81 for ( int i = 0 ; i < (int) N ; ++i ) {
82 for ( int j = 0 ; j < (int) N ; ++j ) {
83 for ( int k = 0 ; k < (int) N ; ++k ) {
84
85 const ordinal row = map_fem_graph_coord((int)N,i,j,k);
86
87 graph[row].reserve(27);
88
89 for ( int ii = -1 ; ii < 2 ; ++ii ) {
90 for ( int jj = -1 ; jj < 2 ; ++jj ) {
91 for ( int kk = -1 ; kk < 2 ; ++kk ) {
92 if ( 0 <= i + ii && i + ii < (int) N &&
93 0 <= j + jj && j + jj < (int) N &&
94 0 <= k + kk && k + kk < (int) N ) {
95 ordinal col = map_fem_graph_coord((int)N,i+ii,j+jj,k+kk);
96
97 graph[row].push_back(col);
98 }
99 }}}
100 total += graph[row].size();
101 }}}
102
103 return total ;
104}
105
106Teuchos::Array<double>
108 const Teuchos::Array<int> & var_degree ,
109 const int nGrid ,
110 const int iterCount ,
111 const bool print_flag ,
112 const bool test_block ,
113 const bool check )
114{
115 typedef double value_type ;
116 typedef Stokhos::OneDOrthogPolyBasis<int,value_type> abstract_basis_type;
118 typedef Stokhos::CompletePolynomialBasis<int,value_type> product_basis_type;
120
121 using Teuchos::rcp;
122 using Teuchos::RCP;
123 using Teuchos::Array;
124 using Teuchos::ParameterList;
125
126 // Create a communicator for Epetra objects
127 RCP<const Epetra_Comm> globalComm;
128#ifdef HAVE_MPI
129 RCP<const Epetra_MpiComm> globalMpiComm =
130 rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
131 //globalMpiComm->Print(std::cout);
132 globalComm = globalMpiComm;
133#else
134 RCP<const Epetra_SerialComm> globalSerialComm =
135 rcp(new Epetra_SerialComm);
136 //globalSerialComm->Print(std::cout);
137 globalComm = globalSerialComm;
138#endif
139
140 //int MyPID = globalComm->MyPID();
141
142 // Create Stochastic Galerkin basis and expansion
143 const size_t num_KL = var_degree.size();
144 Array< RCP<const abstract_basis_type> > bases(num_KL);
145 for (size_t i=0; i<num_KL; i++)
146 bases[i] = rcp(new basis_type(var_degree[i],true));
147 RCP<const product_basis_type> basis =
148 rcp(new product_basis_type(bases, 1e-12));
149 const size_t stoch_length = basis->size();
150 RCP<Cijk_type> Cijk = basis->computeTripleProductTensor();
151
152 // Create stochastic parallel distribution
153 ParameterList parallelParams;
154 RCP<Stokhos::ParallelData> sg_parallel_data =
155 rcp(new Stokhos::ParallelData(basis, Cijk, globalComm, parallelParams));
156 RCP<const EpetraExt::MultiComm> sg_comm =
157 sg_parallel_data->getMultiComm();
158 RCP<const Epetra_Comm> app_comm =
159 sg_parallel_data->getSpatialComm();
160 RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
161 sg_parallel_data->getEpetraCijk();
162 RCP<const Epetra_BlockMap> stoch_row_map =
163 epetraCijk->getStochasticRowMap();
164
165 //------------------------------
166 // Generate FEM graph:
167
168 Teuchos::Array< Teuchos::Array<int> > fem_graph ;
169
170 const size_t fem_length = nGrid * nGrid * nGrid ;
171 generate_fem_graph( nGrid , fem_graph );
172
173 //------------------------------
174
175 // Generate Epetra objects
176 RCP<const Epetra_Map> x_map = rcp(new Epetra_Map(static_cast<int>(fem_length), 0, *app_comm));
177 RCP<const Epetra_Map> sg_x_map =
178 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
179 *x_map, *stoch_row_map, *sg_comm));
180
181 RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy, *x_map, 27));
182 int *my_GIDs = x_map->MyGlobalElements();
183 int num_my_GIDs = x_map->NumMyElements();
184 for (int i=0; i<num_my_GIDs; ++i) {
185 int row = my_GIDs[i];
186 int num_indices = fem_graph[row].size();
187 int *indices = &(fem_graph[row][0]);
188 graph->InsertGlobalIndices(row, num_indices, indices);
189 }
190 graph->FillComplete();
191 int nnz = graph->NumGlobalNonzeros();
192
193 RCP<ParameterList> sg_op_params = rcp(new ParameterList);
194 RCP<Stokhos::MatrixFreeOperator> sg_A =
195 rcp(new Stokhos::MatrixFreeOperator(sg_comm, basis, epetraCijk,
196 x_map, x_map, sg_x_map, sg_x_map,
197 sg_op_params));
198 RCP<Epetra_BlockMap> sg_A_overlap_map =
199 rcp(new Epetra_LocalMap(
200 static_cast<int>(stoch_length), 0, *(sg_parallel_data->getStochasticComm())));
201 RCP< Stokhos::EpetraOperatorOrthogPoly > A_sg_blocks =
203 basis, sg_A_overlap_map, x_map, x_map, sg_x_map, sg_comm));
204 for (unsigned int i=0; i<stoch_length; i++) {
205 RCP<Epetra_CrsMatrix> A = rcp(new Epetra_CrsMatrix(Copy, *graph));
206 A->FillComplete();
207 A->PutScalar(1.0);
208 A_sg_blocks->setCoeffPtr(i, A);
209 }
210 sg_A->setupOperator(A_sg_blocks);
211
212 RCP<Stokhos::EpetraVectorOrthogPoly> sg_x =
214 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
215 RCP<Stokhos::EpetraVectorOrthogPoly> sg_y =
217 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
218 sg_x->init(1.0);
219 sg_y->init(0.0);
220
221
222 // Time apply
223 Teuchos::Time clock("apply") ;
224 clock.start();
225 for ( int iter = 0 ; iter < iterCount ; ++iter ) {
226 sg_A->Apply( *(sg_x->getBlockVector()), *(sg_y->getBlockVector()) );
227 }
228 clock.stop();
229 double seconds_per_iter =
230 clock.totalElapsedTime() / ((double) iterCount );
231
232 // Averge time across proc's
233 double average_seconds_per_iter;
234 globalComm->SumAll(&seconds_per_iter, &average_seconds_per_iter, 1);
235 average_seconds_per_iter /= globalComm->NumProc();
236
237 // Compute number of fem mat-vec's
238 int n_apply = 0;
239 int n_add = 0;
240 for (Cijk_type::k_iterator k_it=Cijk->k_begin();
241 k_it!=Cijk->k_end(); ++k_it) {
242 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
243 j_it != Cijk->j_end(k_it); ++j_it) {
244 ++n_apply;
245 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
246 i_it != Cijk->i_end(j_it); ++i_it) {
247 ++n_add;
248 }
249 }
250 }
251
252 const double flops = 1.0e-9*(2.0*static_cast<double>(n_apply)*nnz +
253 static_cast<double>(n_add)*fem_length);
254
255 //------------------------------
256
257 Teuchos::Array<double> perf(8);
258 perf[0] = stoch_length;
259 perf[1] = fem_length;
260 perf[2] = stoch_length * fem_length;
261 perf[3] = Cijk->num_entries();
262 perf[4] = nnz;
263 perf[5] = average_seconds_per_iter ;
264 perf[6] = flops/average_seconds_per_iter;
265 perf[7] = flops;
266
267 return perf;
268}
269
270void performance_test_driver_epetra( const int pdeg ,
271 const int minvar ,
272 const int maxvar ,
273 const int nGrid ,
274 const int nIter ,
275 const bool print ,
276 const bool test_block ,
277 const bool check ,
278 Teuchos::FancyOStream& out)
279{
280 out.precision(8);
281
282 //------------------------------
283
284 out << std::endl
285 << "\"#nGrid\" , \"NNZ\" "
286 << "\"#Variable\" , \"PolyDegree\" , \"#Bases\" , "
287 << "\"#TensorEntry\" , "
288 << "\"VectorSize\" , "
289 << "\"Original-Matrix-Free-Block-MXV-Time\" , "
290 << "\"Original-Matrix-Free-Block-MXV-GFLOPS\" , "
291 << std::endl ;
292
293 for ( int nvar = minvar ; nvar <= maxvar ; ++nvar ) {
294 Teuchos::Array<int> var_degree( nvar , pdeg );
295
296 const Teuchos::Array<double> perf_original_mat_free_epetra =
297 test_original_matrix_free_epetra( var_degree , nGrid , nIter , print , test_block , check );
298
299 out << nGrid << " , "
300 << perf_original_mat_free_epetra[4] << " , "
301 << nvar << " , " << pdeg << " , "
302 << perf_original_mat_free_epetra[0] << " , "
303 << perf_original_mat_free_epetra[3] << " , "
304 << perf_original_mat_free_epetra[2] << " , "
305 << perf_original_mat_free_epetra[5] << " , "
306 << perf_original_mat_free_epetra[6] << " , "
307 << std::endl ;
308
309 }
310
311
312}
313
314int main(int argc, char *argv[])
315{
316 bool success = true;
317
318 try {
319 Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
320 Teuchos::RCP< Teuchos::FancyOStream > out =
321 Teuchos::VerboseObjectBase::getDefaultOStream();
322
323 // Setup command line options
324 Teuchos::CommandLineProcessor CLP;
325 int p = 3;
326 CLP.setOption("p", &p, "Polynomial order");
327 int d_min = 1;
328 CLP.setOption("dmin", &d_min, "Starting stochastic dimension");
329 int d_max = 12;
330 CLP.setOption("dmax", &d_max, "Ending stochastic dimension");
331 int nGrid = 64;
332 CLP.setOption("n", &nGrid, "Number of spatial grid points in each dimension");
333 int nIter = 1;
334 CLP.setOption("niter", &nIter, "Number of iterations");
335 bool test_block = true;
336 CLP.setOption("block", "no-block", &test_block, "Use block algorithm");
337 CLP.parse( argc, argv );
338
339 bool print = false ;
340 bool check = false;
342 p , d_min , d_max , nGrid , nIter , print , test_block , check , *out );
343 }
344 TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
345
346 if (!success)
347 return -1;
348 return 0;
349}
350
Copy
int main(int argc, char *argv[])
ordinal generate_fem_graph(ordinal N, Teuchos::Array< Teuchos::Array< ordinal > > &graph)
void performance_test_driver_epetra(const int pdeg, const int minvar, const int maxvar, const int nGrid, const int nIter, const bool print, const bool test_block, const bool check, Teuchos::FancyOStream &out)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Teuchos::Array< double > test_original_matrix_free_epetra(const Teuchos::Array< int > &var_degree, const int nGrid, const int iterCount, const bool print_flag, const bool test_block, const bool check)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Legendre polynomial basis.
An Epetra operator representing the block stochastic Galerkin operator.
Abstract base class for 1-D orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
Stokhos::LegendreBasis< int, double > basis_type