Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TestEpetra.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 
47 
48 #include "Stokhos_Epetra.hpp"
50 #include "EpetraExt_BlockUtility.h"
51 
52 #include "impl/Kokkos_Timer.hpp"
53 
54 #ifdef HAVE_MPI
55 #include "Epetra_MpiComm.h"
56 #else
57 #include "Epetra_SerialComm.h"
58 #endif
59 
60 using Teuchos::rcp;
61 using Teuchos::RCP;
62 using Teuchos::Array;
64 
65 template< typename IntType >
66 inline
67 IntType map_fem_graph_coord( const IntType & N ,
68  const IntType & i ,
69  const IntType & j ,
70  const IntType & k )
71 {
72  return k + N * ( j + N * i );
73 }
74 
75 template < typename ordinal >
76 inline
77 ordinal generate_fem_graph( ordinal N ,
78  std::vector< std::vector<ordinal> > & graph )
79 {
80  graph.resize( N * N * N , std::vector<ordinal>() );
81 
82  ordinal total = 0 ;
83 
84  for ( int i = 0 ; i < (int) N ; ++i ) {
85  for ( int j = 0 ; j < (int) N ; ++j ) {
86  for ( int k = 0 ; k < (int) N ; ++k ) {
87 
88  const ordinal row = map_fem_graph_coord((int)N,i,j,k);
89 
90  graph[row].reserve(27);
91 
92  for ( int ii = -1 ; ii < 2 ; ++ii ) {
93  for ( int jj = -1 ; jj < 2 ; ++jj ) {
94  for ( int kk = -1 ; kk < 2 ; ++kk ) {
95  if ( 0 <= i + ii && i + ii < (int) N &&
96  0 <= j + jj && j + jj < (int) N &&
97  0 <= k + kk && k + kk < (int) N ) {
98  ordinal col = map_fem_graph_coord((int)N,i+ii,j+jj,k+kk);
99 
100  graph[row].push_back(col);
101  }
102  }}}
103  total += graph[row].size();
104  }}}
105 
106  return total ;
107 }
108 
109 void
110 run_test(const int p, const int d, const int nGrid, const int nIter,
111  const RCP<const Epetra_Comm>& globalComm,
112  const RCP<const Epetra_Map>& map,
113  const RCP<Epetra_CrsGraph>& graph)
114 {
115  typedef double value_type ;
116  typedef Stokhos::OneDOrthogPolyBasis<int,value_type> abstract_basis_type;
119  typedef Stokhos::TotalOrderBasis<int,value_type,less_type> product_basis_type;
121 
122  // Create Stochastic Galerkin basis and expansion
124  for (int i=0; i<d; i++)
125  bases[i] = rcp(new basis_type(p,true));
126  RCP< product_basis_type> basis = rcp(new product_basis_type(bases, 1e-12));
127  int stoch_length = basis->size();
128  RCP<Cijk_type> Cijk = basis->computeTripleProductTensor();
129 
130  // Create stochastic parallel distribution
131  ParameterList parallelParams;
132  RCP<Stokhos::ParallelData> sg_parallel_data =
133  rcp(new Stokhos::ParallelData(basis, Cijk, globalComm, parallelParams));
135  sg_parallel_data->getMultiComm();
136  RCP<const Epetra_Comm> app_comm =
137  sg_parallel_data->getSpatialComm();
139  sg_parallel_data->getEpetraCijk();
140  RCP<const Epetra_BlockMap> stoch_row_map =
141  epetraCijk->getStochasticRowMap();
142 
143  // Generate Epetra objects
144  RCP<const Epetra_Map> sg_map =
145  rcp(EpetraExt::BlockUtility::GenerateBlockMap(
146  *map, *stoch_row_map, *sg_comm));
147  RCP<ParameterList> sg_op_params = rcp(new ParameterList);
149  rcp(new Stokhos::MatrixFreeOperator(sg_comm, basis, epetraCijk,
150  map, map, sg_map, sg_map,
151  sg_op_params));
152  RCP<Epetra_BlockMap> sg_A_overlap_map =
153  rcp(new Epetra_LocalMap(
154  stoch_length, 0, *(sg_parallel_data->getStochasticComm())));
157  basis, sg_A_overlap_map, map, map, sg_map, sg_comm));
158  for (int i=0; i<stoch_length; i++) {
160  A->PutScalar(1.0);
161  A->FillComplete();
162  A_sg_blocks->setCoeffPtr(i, A);
163  }
164  sg_A->setupOperator(A_sg_blocks);
165 
168  basis, stoch_row_map, map, sg_map, sg_comm));
171  basis, stoch_row_map, map, sg_map, sg_comm));
172  sg_x->init(1.0);
173  sg_y->init(0.0);
174 
175  // Apply operator
176  Kokkos::Impl::Timer clock;
177  for (int iter=0; iter<nIter; ++iter)
178  sg_A->Apply( *(sg_x->getBlockVector()), *(sg_y->getBlockVector()) );
179 
180  const double t = clock.seconds() / ((double) nIter );
181  const double flops = sg_A->countApplyFlops();
182  const double gflops = 1.0e-9 * flops / t;
183 
184  if (globalComm->MyPID() == 0)
185  std::cout << nGrid << " , "
186  << d << " , "
187  << p << " , "
188  << t << " , "
189  << gflops << " , "
190  << std::endl;
191 }
192 
193 int main(int argc, char *argv[])
194 {
195  bool success = true;
196 
197 // Initialize MPI
198 #ifdef HAVE_MPI
199  MPI_Init(&argc,&argv);
200 #endif
201 
202  try {
203 
204 // Create a communicator for Epetra objects
205 #ifdef HAVE_MPI
206  RCP<const Epetra_Comm> globalComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
207 #else
208  RCP<const Epetra_Comm> globalComm = rcp(new Epetra_SerialComm);
209 #endif
210 
211  // Print header
212  if (globalComm->MyPID() == 0)
213  std::cout << std::endl
214  << "\"#nGrid\" , "
215  << "\"#Variable\" , "
216  << "\"PolyDegree\" , "
217  << "\"MXV Time\" , "
218  << "\"MXV GFLOPS\" , "
219  << std::endl;
220 
221  const int nIter = 1;
222  const int nGrid = 32;
223 
224  // Generate FEM graph:
225  const int fem_length = nGrid * nGrid * nGrid;
226  RCP<const Epetra_Map> map = rcp(new Epetra_Map(fem_length, 0, *globalComm));
227  std::vector< std::vector<int> > fem_graph;
228  generate_fem_graph(nGrid, fem_graph);
229  RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy, *map, 27));
230  int *my_GIDs = map->MyGlobalElements();
231  int num_my_GIDs = map->NumMyElements();
232  for (int i=0; i<num_my_GIDs; ++i) {
233  int row = my_GIDs[i];
234  int num_indices = fem_graph[row].size();
235  int *indices = &(fem_graph[row][0]);
236  graph->InsertGlobalIndices(row, num_indices, indices);
237  }
238  graph->FillComplete();
239 
240  {
241  const int p = 3;
242  for (int d=1; d<=12; ++d)
243  run_test(p, d, nGrid, nIter, globalComm, map, graph);
244  }
245 
246  {
247  const int p = 5;
248  for (int d=1; d<=6; ++d)
249  run_test(p, d, nGrid, nIter, globalComm, map, graph);
250  }
251 
252  }
253  TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
254 
255 #ifdef HAVE_MPI
256  MPI_Finalize() ;
257 #endif
258 
259  if (!success)
260  return -1;
261  return 0 ;
262 }
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.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Teuchos::RCP< const EpetraExt::MultiComm > getMultiComm() const
Get global comm.
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
int MyGlobalElements(int *MyGlobalElementList) const
void init(const value_type &val)
Initialize coefficients.
Teuchos::RCP< EpetraExt::BlockVector > getBlockVector()
Get block vector.
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Definition: TestEpetra.cpp:67
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
virtual int MyPID() const =0
Teuchos::RCP< const Epetra_Comm > getStochasticComm() const
Get stochastic comm.
int FillComplete(bool OptimizeDataStorage=true)
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > getEpetraCijk() const
Get Epetra Cijk.
Stokhos::LegendreBasis< int, double > basis_type
int PutScalar(double ScalarConstant)
int NumMyElements() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const Epetra_BlockMap > getStochasticRowMap() const
Get stochastic row map.
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)
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector Input in Result as described ...
std::vector< double > run_test(const size_t num_cpu, const size_t num_core_per_cpu, const size_t num_threads_per_core, const size_t p, const size_t d, const size_t nGrid, const size_t nIter, const bool symmetric, SG_Alg sg_alg, const std::vector< double > &perf1=std::vector< double >())
Definition: HostScaling.cpp:62
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
int main(int argc, char **argv)
Abstract base class for 1-D orthogonal polynomials.
double countApplyFlops() const
Return number of FLOPS for each call to Apply()
A comparison functor implementing a strict weak ordering based lexographic ordering.
Teuchos::RCP< const Epetra_Comm > getSpatialComm() const
Get spatial comm.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.