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();
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)
115 typedef double value_type ;
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));
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);
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)