Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 ,
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 
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;
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 
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 ,
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);
321  Teuchos::VerboseObjectBase::getDefaultOStream();
322 
323  // Setup command line options
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.
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
Definition: TestEpetra.cpp:77
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Definition: TestEpetra.cpp:67
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
Stokhos::LegendreBasis< int, double > basis_type
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double stop()
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
int main(int argc, char **argv)
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.
size_type size() const
double totalElapsedTime(bool readCurrentTime=false) const