Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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 "Kokkos_Timer.hpp"
53
54#ifdef HAVE_MPI
55#include "Epetra_MpiComm.h"
56#else
57#include "Epetra_SerialComm.h"
58#endif
59
60using Teuchos::rcp;
61using Teuchos::RCP;
62using Teuchos::Array;
63using Teuchos::ParameterList;
64
65template< typename IntType >
66inline
67IntType 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
75template < typename ordinal >
76inline
77ordinal 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
109void
110run_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;
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::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
193int 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}
Copy
int main(int argc, char *argv[])
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
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)
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.
A comparison functor implementing a strict weak ordering based lexographic ordering.
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.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Stokhos::LegendreBasis< int, double > basis_type