Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
46 #include "Teuchos_StandardCatchMacros.hpp"
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;
63 using Teuchos::ParameterList;
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
123  Array< RCP<const abstract_basis_type> > bases(d);
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));
134  RCP<const EpetraExt::MultiComm> sg_comm =
135  sg_parallel_data->getMultiComm();
136  RCP<const Epetra_Comm> app_comm =
137  sg_parallel_data->getSpatialComm();
138  RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
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);
148  RCP<Stokhos::MatrixFreeOperator> sg_A =
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())));
155  RCP< Stokhos::EpetraOperatorOrthogPoly > A_sg_blocks =
157  basis, sg_A_overlap_map, map, map, sg_map, sg_comm));
158  for (int i=0; i<stoch_length; i++) {
159  RCP<Epetra_CrsMatrix> A = rcp(new Epetra_CrsMatrix(Copy, *graph));
160  A->PutScalar(1.0);
161  A->FillComplete();
162  A_sg_blocks->setCoeffPtr(i, A);
163  }
164  sg_A->setupOperator(A_sg_blocks);
165 
166  RCP<Stokhos::EpetraVectorOrthogPoly> sg_x =
168  basis, stoch_row_map, map, sg_map, sg_comm));
169  RCP<Stokhos::EpetraVectorOrthogPoly> sg_y =
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...
Copy
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
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.
expr expr expr expr j
int main(int argc, char *argv[])
Definition: TestEpetra.cpp:193
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
void run_test(const int p, const int d, const int nGrid, const int nIter, const RCP< const Epetra_Comm > &globalComm, const RCP< const Epetra_Map > &map, const RCP< Epetra_CrsGraph > &graph)
Definition: TestEpetra.cpp:110
Abstract base class for 1-D orthogonal polynomials.
A comparison functor implementing a strict weak ordering based lexographic ordering.