Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
62 template< typename IntType >
63 inline
64 IntType 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 
72 template < typename ordinal >
73 inline
74 ordinal 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 
106 Teuchos::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 
270 void 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 
314 int 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 
An Epetra operator representing the block stochastic Galerkin operator.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Stokhos::LegendreBasis< int, double > basis_type
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
int main(int argc, char *argv[])
expr expr expr expr j
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
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)
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)
Abstract base class for 1-D orthogonal polynomials.
ordinal generate_fem_graph(ordinal N, Teuchos::Array< Teuchos::Array< ordinal > > &graph)