82 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
83 for (
int i=0; i<d; i++) {
86 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
90 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
91 basis->computeTripleProductTensor();
94 Teuchos::RCP<const Epetra_Comm> globalComm;
102 int num_spatial_procs = -1;
103 int num_procs = globalComm->NumProc();
105 num_spatial_procs = num_procs / 2;
106 Teuchos::ParameterList parallelParams;
107 parallelParams.set(
"Number of Spatial Processors", num_spatial_procs);
108 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
111 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
112 sg_parallel_data->getMultiComm();
113 Teuchos::RCP<const Epetra_Comm> app_comm =
114 sg_parallel_data->getSpatialComm();
115 Teuchos::RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
116 sg_parallel_data->getEpetraCijk();
120 Teuchos::RCP<Epetra_Map> x_map =
121 Teuchos::rcp(
new Epetra_Map(num_x, 0, *app_comm));
124 Teuchos::RCP<Epetra_Map> x_overlap_map =
129 Teuchos::RCP<Epetra_Map> f_map =
130 Teuchos::rcp(
new Epetra_Map(num_f, 0, *app_comm));
133 Teuchos::RCP<const Epetra_BlockMap> stoch_row_map =
134 epetraCijk->getStochasticRowMap();
135 sg_x_map = Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
136 *x_map, *stoch_row_map, *sg_comm));
137 sg_f_map = Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
138 *f_map, *stoch_row_map, *sg_comm));
141 const int num_indices = num_x;
142 Teuchos::RCP<Epetra_CrsGraph> graph =
144 int indices[num_indices];
145 for (
int j=0;
j<num_indices;
j++)
146 indices[
j] = x_overlap_map->GID(
j);
147 for (
int i=0; i<f_map->NumMyElements(); i++)
148 graph->InsertGlobalIndices(f_map->GID(i), num_indices, indices);
149 graph->FillComplete(*x_map, *f_map);
152 Teuchos::RCP<Epetra_BlockMap> sg_overlap_map =
155 *(sg_parallel_data->getStochasticComm())));
156 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > mat_sg =
158 basis, sg_overlap_map, x_map, f_map,
sg_f_map, sg_comm));
159 for (
int block=0; block<basis->size(); block++) {
160 Teuchos::RCP<Epetra_CrsMatrix> mat =
162 TEUCHOS_TEST_FOR_EXCEPTION(!mat->IndicesAreLocal(), std::logic_error,
163 "Indices are not local!");
164 double values[num_indices];
165 for (
int i=0; i<f_map->NumMyElements(); i++) {
166 for (
int j=0;
j<num_indices;
j++) {
167 indices[
j] = x_overlap_map->GID(
j);
168 values[
j] = 0.1*(i+1)*(
j+1)*(block+1);
170 mat->ReplaceMyValues(i, num_indices, values, indices);
172 mat->FillComplete(*x_map, *f_map);
173 mat_sg->setCoeffPtr(block, mat);
177 Teuchos::RCP<Teuchos::ParameterList> op_params =
178 Teuchos::rcp(
new Teuchos::ParameterList);
181 sg_comm, basis, epetraCijk, x_map, f_map,