46 #include "Teuchos_StandardCatchMacros.hpp" 50 #include "EpetraExt_BlockUtility.h" 52 #include "impl/Kokkos_Timer.hpp" 63 using Teuchos::ParameterList;
65 template<
typename IntType >
72 return k + N * (
j + N * i );
75 template <
typename ordinal >
78 std::vector< std::vector<ordinal> > & graph )
80 graph.resize( N * N * N , std::vector<ordinal>() );
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 ) {
90 graph[row].reserve(27);
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 ) {
100 graph[row].push_back(col);
103 total += graph[row].size();
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)
123 Array< RCP<const abstract_basis_type> > bases(d);
124 for (
int i=0; i<d; i++)
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();
131 ParameterList parallelParams;
132 RCP<Stokhos::ParallelData> sg_parallel_data =
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();
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 =
150 map, map, sg_map, sg_map,
152 RCP<Epetra_BlockMap> sg_A_overlap_map =
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++) {
162 A_sg_blocks->setCoeffPtr(i, A);
164 sg_A->setupOperator(A_sg_blocks);
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));
176 Kokkos::Impl::Timer clock;
177 for (
int iter=0; iter<nIter; ++iter)
178 sg_A->Apply( *(sg_x->getBlockVector()), *(sg_y->getBlockVector()) );
180 const double t = clock.seconds() / ((
double) nIter );
181 const double flops = sg_A->countApplyFlops();
182 const double gflops = 1.0e-9 * flops / t;
184 if (globalComm->MyPID() == 0)
185 std::cout << nGrid <<
" , " 199 MPI_Init(&argc,&
argv);
206 RCP<const Epetra_Comm> globalComm = rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
212 if (globalComm->MyPID() == 0)
213 std::cout << std::endl
215 <<
"\"#Variable\" , " 216 <<
"\"PolyDegree\" , " 218 <<
"\"MXV GFLOPS\" , " 222 const int nGrid = 32;
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;
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);
238 graph->FillComplete();
242 for (
int d=1; d<=12; ++d)
243 run_test(p, d, nGrid, nIter, globalComm, map, graph);
248 for (
int d=1; d<=6; ++d)
249 run_test(p, d, nGrid, nIter, globalComm, map, graph);
253 TEUCHOS_STANDARD_CATCH_STATEMENTS(
true, std::cerr, success);
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...
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)
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[])
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)
Abstract base class for 1-D orthogonal polynomials.
A comparison functor implementing a strict weak ordering based lexographic ordering.