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())));
144 full_op.PutScalar(0.0);
145 full_op.setupOperator(W_sg_blocks);
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);
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);