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);
207 success = std::fabs(nrm) <
setup.
tol;
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);
225 success = std::fabs(nrm) <
setup.
tol;
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);
int main(int argc, char *argv[])
int NormInf(double *Result) const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
static void SetTracebackMode(int TracebackModeValue)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Legendre polynomial basis.
An Epetra operator representing the block stochastic Galerkin operator.
TEUCHOS_UNIT_TEST(Stokhos_MatrixFreeOperator, ApplyUnitTest)
Teuchos::RCP< Stokhos::MatrixFreeOperator > mat_free_op
Teuchos::RCP< Stokhos::FullyAssembledOperator > assembled_op
Teuchos::RCP< const Epetra_Map > sg_x_map
Teuchos::RCP< const Epetra_Map > sg_f_map