47 #include "Teuchos_UnitTestHarness.hpp" 48 #include "Teuchos_TestingHelpers.hpp" 49 #include "Teuchos_UnitTestRepository.hpp" 50 #include "Teuchos_GlobalMPISession.hpp" 53 #include "EpetraExt_BlockUtility.h" 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,
204 diff.
Update(1.0, result1, -1.0, result2, 0.0);
208 out <<
"Apply infinity norm of difference: " << nrm << std::endl;
209 out <<
"Matrix-free result = " << std::endl << result1 << std::endl
210 <<
"Assebled result = " << std::endl << result2 << std::endl;
222 diff.
Update(1.0, result1, -1.0, result2, 0.0);
226 out <<
"Apply-transpose infinity norm of difference: " << nrm << std::endl;
227 out <<
"Matrix-free result = " << std::endl << result1 << std::endl
228 <<
"Assebled result = " << std::endl << result2 << std::endl;
234 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
236 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
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.
TEUCHOS_UNIT_TEST(Stokhos_MatrixFreeOperator, ApplyUnitTest)
static void SetTracebackMode(int TracebackModeValue)
int NormInf(double *Result) const
Teuchos::RCP< Stokhos::MatrixFreeOperator > mat_free_op
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Teuchos::RCP< const Epetra_Map > sg_f_map
Teuchos::RCP< Stokhos::FullyAssembledOperator > assembled_op
int main(int argc, char *argv[])
Teuchos::RCP< const Epetra_Map > sg_x_map
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...