45 #include "Thyra_VectorStdOps.hpp" 48 #include "Thyra_DefaultBlockedLinearOp.hpp" 49 #include "Thyra_ProductVectorBase.hpp" 50 #include "Thyra_SpmdVectorSpaceBase.hpp" 51 #include "Thyra_DetachedSpmdVectorView.hpp" 52 #include "Thyra_TestingTools.hpp" 53 #include "Thyra_LinearOpTester.hpp" 54 #include "Trilinos_Util_CrsMatrixGallery.h" 55 #include "Teuchos_GlobalMPISession.hpp" 56 #include "Teuchos_VerboseObject.hpp" 57 #include "Teuchos_XMLParameterListHelpers.hpp" 58 #include "Teuchos_CommandLineProcessor.hpp" 59 #include "Teuchos_StandardCatchMacros.hpp" 62 # include "Epetra_MpiComm.h" 64 # include "Epetra_SerialComm.h" 66 #include "Epetra_Vector.h" 67 #include "Epetra_CrsMatrix.h" 69 #include "Teuchos_UnitTestHarness.hpp" 78 using Teuchos::rcp_dynamic_cast;
79 using Teuchos::inOutArg;
87 Epetra_MpiComm comm(MPI_COMM_WORLD);
89 Epetra_SerialComm comm;
92 out <<
"\nRunning on " << comm.NumProc() <<
" processors\n";
97 out <<
"Using Trilinos_Util to create test matrices\n";
100 Trilinos_Util::CrsMatrixGallery FGallery(
"recirc_2d",comm,
false);
101 FGallery.Set(
"nx",nx);
102 FGallery.Set(
"ny",ny);
103 RCP<Epetra_CrsMatrix> F = rcp(FGallery.GetMatrix(),
false);
105 Trilinos_Util::CrsMatrixGallery CGallery(
"laplace_2d",comm,
false);
106 CGallery.Set(
"nx",nx);
107 CGallery.Set(
"ny",ny);
108 RCP<Epetra_CrsMatrix> C = rcp(CGallery.GetMatrix(),
false);
110 Trilinos_Util::CrsMatrixGallery BGallery(
"diag",comm,
false);
111 BGallery.Set(
"nx",nx*ny);
112 BGallery.Set(
"a",5.0);
113 RCP<Epetra_CrsMatrix> B = rcp(BGallery.GetMatrix(),
false);
115 Trilinos_Util::CrsMatrixGallery BtGallery(
"diag",comm,
false);
116 BtGallery.Set(
"nx",nx*ny);
117 BtGallery.Set(
"a",3.0);
118 RCP<Epetra_CrsMatrix> Bt = rcp(BtGallery.GetMatrix(),
false);
121 out <<
"Building block2x2 Thyra matrix ... wrapping in EpetraOperatorWrapper\n";
122 const RCP<const LinearOpBase<double> > A =
123 Thyra::block2x2<double>(
124 Thyra::epetraLinearOp(F),
125 Thyra::epetraLinearOp(Bt),
126 Thyra::epetraLinearOp(B),
127 Thyra::epetraLinearOp(C),
131 const RCP<Thyra::EpetraOperatorWrapper> epetra_A =
135 const Epetra_Map & rangeMap = epetra_A->OperatorRangeMap();
136 const Epetra_Map & domainMap = epetra_A->OperatorDomainMap();
139 TEST_EQUALITY(rangeMap.NumGlobalElements(), 2*nx*ny);
140 TEST_EQUALITY(domainMap.NumGlobalElements(), 2*nx*ny);
143 TEST_EQUALITY(rangeMap.NumGlobalElements()-1, rangeMap.MaxAllGID());
144 TEST_EQUALITY(domainMap.NumGlobalElements()-1, domainMap.MaxAllGID());
148 const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
149 Thyra::randomize(-100.0, 100.0, tv.ptr());
150 const RCP<const VectorBase<double> > tv_0 =
151 Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
152 const RCP<const VectorBase<double> > tv_1 =
153 Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
154 const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
155 const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
157 int off_0 = vv_0.globalOffset();
158 int off_1 = vv_1.globalOffset();
161 Epetra_Vector ev(epetra_A->OperatorDomainMap());
162 epetra_A->copyThyraIntoEpetra(*tv, ev);
165 TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
166 const int numMyElements = domainMap.NumMyElements();
168 for(
int i=0; i < numMyElements; i++) {
169 int gid = domainMap.GID(i);
171 tval = vv_0[gid-off_0];
173 tval = vv_1[gid-off_1-nx*ny];
174 TEST_EQUALITY(ev[i], tval);
181 Epetra_Vector ev(epetra_A->OperatorDomainMap());
185 const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
186 const RCP<const VectorBase<double> > tv_0 =
187 Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
188 const RCP<const VectorBase<double> > tv_1 =
189 Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
190 const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
191 const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
193 int off_0 = rcp_dynamic_cast<
const Thyra::SpmdVectorSpaceBase<double> >(
194 tv_0->space())->localOffset();
195 int off_1 = rcp_dynamic_cast<
const Thyra::SpmdVectorSpaceBase<double> >(
196 tv_1->space())->localOffset();
198 epetra_A->copyEpetraIntoThyra(ev, tv.ptr());
201 TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
202 int numMyElements = domainMap.NumMyElements();
204 for(
int i=0;i<numMyElements;i++) {
205 int gid = domainMap.GID(i);
207 tval = vv_0[gid-off_0];
209 tval = vv_1[gid-off_1-nx*ny];
210 TEST_EQUALITY(ev[i], tval);
215 const RCP<const LinearOpBase<double> > thyraEpetraOp = epetraLinearOp(epetra_A);
216 LinearOpTester<double> linearOpTester;
217 linearOpTester.show_all_tests(
true);
218 const bool checkResult = linearOpTester.check(*thyraEpetraOp, inOutArg(out));
219 TEST_ASSERT(checkResult);
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Implements the Epetra_Operator interface with a Thyra LinearOperator.