44 #include <Teuchos_ConfigDefs.hpp> 45 #include <Teuchos_UnitTestHarness.hpp> 46 #include <Teuchos_TimeMonitor.hpp> 47 #include <Teuchos_RCP.hpp> 63 #include "EpetraExt_BlockUtility.h" 64 #include "EpetraExt_RowMatrixOut.h" 69 Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
75 int numProc = comm->NumProc();
79 bool full_expansion =
false;
81 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
buildBasis(num_KL,porder);
82 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
83 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data;
84 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion;
87 Cijk = basis->computeTripleProductTensor();
89 Cijk = basis->computeLinearTripleProductTensor();
91 Teuchos::ParameterList parallelParams;
92 parallelParams.set(
"Number of Spatial Processors", numProc);
99 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm = sg_parallel_data->getMultiComm();
102 Teuchos::RCP<Epetra_Map> determRowMap = Teuchos::rcp(
new Epetra_Map(-1,10,0,*comm));
103 Teuchos::RCP<Epetra_CrsGraph> determGraph = Teuchos::rcp(
new Epetra_CrsGraph(
Copy,*determRowMap,1));
104 for(
int row=0;row<determRowMap->NumMyElements();row++) {
105 int gid = determRowMap->GID(row);
106 determGraph->InsertGlobalIndices(gid,1,&gid);
108 for(
int row=1;row<determRowMap->NumMyElements()-1;row++) {
109 int gid = determRowMap->GID(row);
110 int indices[2] = {gid-1,gid+1};
111 determGraph->InsertGlobalIndices(gid,2,indices);
113 determGraph->FillComplete();
115 Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(
new Teuchos::ParameterList);
116 params->set(
"Scale Operator by Inverse Basis Norms",
false);
117 params->set(
"Include Mean",
true);
118 params->set(
"Only Use Linear Terms",
false);
120 Teuchos::RCP<Stokhos::EpetraSparse3Tensor> epetraCijk =
122 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> W_sg_blocks =
124 for(
int i=0; i<W_sg_blocks->size(); i++) {
126 crsMat->PutScalar(1.0 + i);
127 W_sg_blocks->setCoeffPtr(i,crsMat);
130 Teuchos::RCP<const Epetra_Map> sg_map =
131 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
132 *determRowMap, *(epetraCijk->getStochasticRowMap()),
133 *(epetraCijk->getMultiComm())));
141 op.setupOperator(W_sg_blocks);
144 full_op.PutScalar(0.0);
150 for(
int i=0;i<100;i++) {
152 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_vec_blocks =
154 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f_vec_blocks =
156 Teuchos::RCP<Epetra_Vector> x_vec_blocked = x_vec_blocks->getBlockVector();
157 Teuchos::RCP<Epetra_Vector> f_vec_blocked = f_vec_blocks->getBlockVector();
158 x_vec_blocked->Random();
159 f_vec_blocked->PutScalar(0.0);
162 Teuchos::RCP<Epetra_Vector> x_vec_inter = Teuchos::rcp(
new Epetra_Vector(op.OperatorDomainMap()));
163 Teuchos::RCP<Epetra_Vector> f_vec_inter = Teuchos::rcp(
new Epetra_Vector(op.OperatorRangeMap()));
164 Teuchos::RCP<Epetra_Vector> f_vec_blk_inter = Teuchos::rcp(
new Epetra_Vector(op.OperatorRangeMap()));
166 f_vec_inter->PutScalar(0.0);
168 full_op.
Apply(*x_vec_blocked,*f_vec_blocked);
169 op.Apply(*x_vec_inter,*f_vec_inter);
176 double true_norm = 0.0;
177 f_vec_blk_inter->NormInf(&true_norm);
178 f_vec_blk_inter->Update(-1.0,*f_vec_inter,1.0);
179 f_vec_blk_inter->NormInf(&error);
181 out <<
"rel error = " << error/true_norm <<
" ( " << true_norm <<
" ), ";
182 result &= (error/true_norm < 1e-14);
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
TEUCHOS_UNIT_TEST(interlaced_op, test)
Orthogonal polynomial expansions limited to algebraic operations.
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Wrap Apply() to add a timer.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< int, double > > buildBasis(int num_KL, int porder)
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
static void copyToInterlacedVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x)