46 #include "Thyra_DefaultScaledAdjointLinearOp.hpp" 47 #include "Thyra_DefaultMultipliedLinearOp.hpp" 48 #include "Thyra_DefaultDiagonalLinearOp.hpp" 49 #include "Thyra_VectorStdOps.hpp" 50 #include "Thyra_TestingTools.hpp" 51 #include "Thyra_LinearOpTester.hpp" 54 #include "EpetraExt_readEpetraLinearSystem.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" 67 #include "Teuchos_UnitTestHarness.hpp" 76 using Thyra::epetraExtDiagScaledMatProdTransformer;
77 using Thyra::VectorBase;
78 using Thyra::LinearOpBase;
79 using Thyra::createMember;
80 using Thyra::LinearOpTester;
82 using Thyra::multiply;
83 using Thyra::diagonal;
86 std::string matrixFile =
"";
87 std::string matrixFile2 =
"";
92 Teuchos::UnitTestRepository::getCLP().setOption(
93 "matrix-file", &matrixFile,
94 "Defines the Epetra_CrsMatrix to read in." );
95 Teuchos::UnitTestRepository::getCLP().setOption(
96 "matrix-file-2", &matrixFile2,
97 "Defines the second Epetra_CrsMatrix to read in." );
101 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
102 buildBDGOperator(
int scenario,
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & B,
103 const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & G,
104 const double vecScale, std::ostream & out)
111 RCP<const LinearOpBase<double> > M;
112 RCP<VectorBase<double> > d;
114 d = createMember(B->domain(),
"d");
116 d = createMember(B->range(),
"d");
117 V_S( d.ptr(), vecScale );
123 M = multiply( scale(scalea,B), diagonal(d), scale(scaleb,G),
"M" );
124 out <<
" Testing B*D*G" << std::endl;
127 M = multiply( scale(scalea,B), diagonal(d), adjoint(G),
"M" );
128 out <<
" Testing B*D*adj(G)" << std::endl;
131 M = multiply( adjoint(B), diagonal(d), scale(scaleb,G),
"M" );
132 out <<
" Testing adj(B)*D*G" << std::endl;
135 M = multiply( adjoint(B), diagonal(d), adjoint(G),
"M" );
136 out <<
" Testing adj(B)*D*adj(G)" << std::endl;
139 M = multiply( B, scale(scaled,diagonal(d)), G,
"M" );
140 out <<
" Testing B*(52.0*D)*G" << std::endl;
143 TEUCHOS_ASSERT(
false);
148 out <<
"\nM = " << *M;
154 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
155 buildBGOperator(
int scenario,
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & B,
156 const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & G,
162 RCP<const LinearOpBase<double> > M;
168 M = multiply( scale(scalea,B), scale(scaleb,G),
"M" );
169 out <<
" Testing B*G" << std::endl;
172 M = multiply( scale(scalea,B), adjoint(G),
"M" );
173 out <<
" Testing B*adj(G)" << std::endl;
176 M = multiply( adjoint(B), scale(scaleb,G),
"M" );
177 out <<
" Testing adj(B)*G" << std::endl;
180 M = multiply( adjoint(B), adjoint(G),
"M" );
181 out <<
" Testing adj(B)*adj(G)" << std::endl;
184 TEUCHOS_ASSERT(
false);
189 out <<
"\nM = " << *M;
194 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
195 buildBDBOperator(
int scenario,
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & B,
196 const double vecScale, std::ostream & out)
198 return buildBDGOperator(scenario,B,B,vecScale,out);
210 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\' ...\n";
213 Epetra_MpiComm comm(MPI_COMM_WORLD);
215 Epetra_SerialComm comm;
217 RCP<Epetra_CrsMatrix> epetra_B;
218 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
224 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
233 for(
int scenario=1;scenario<=5;scenario++) {
234 const RCP<const Thyra::LinearOpBase<double> > M = buildBDBOperator(scenario,B,4.5,out);
240 const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
243 TEST_ASSERT(BtDB_transformer != null);
245 const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
247 BtDB_transformer->transform( *M, M_explicit.ptr() );
249 out <<
"\nM_explicit = " << *M_explicit;
255 LinearOpTester<double> M_explicit_tester;
256 M_explicit_tester.show_all_tests(
true);;
258 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
259 if (!result) success =
false;
271 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\'";
272 out <<
" and from the file \'"<<matrixFile2<<
"\' ...\n";
275 Epetra_MpiComm comm(MPI_COMM_WORLD);
277 Epetra_SerialComm comm;
279 RCP<Epetra_CrsMatrix> epetra_B;
280 RCP<Epetra_CrsMatrix> epetra_G;
281 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
282 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
288 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
289 const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
296 for(
int scenario=1;scenario<=3;scenario++) {
297 RCP<const Thyra::LinearOpBase<double> > M;
298 if(scenario==1 || scenario==3)
299 M = buildBDGOperator(scenario,B,G,4.5,out);
301 M = buildBDGOperator(scenario,G,B,4.5,out);
307 const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
310 TEST_ASSERT(BtDB_transformer != null);
312 const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
314 BtDB_transformer->transform( *M, M_explicit.ptr() );
316 out <<
"\nM_explicit = " << *M_explicit;
322 LinearOpTester<double> M_explicit_tester;
323 M_explicit_tester.show_all_tests(
true);;
325 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
326 if (!result) success =
false;
338 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile2<<
"\' ...\n";
341 Epetra_MpiComm comm(MPI_COMM_WORLD);
343 Epetra_SerialComm comm;
345 RCP<Epetra_CrsMatrix> epetra_G;
346 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
352 const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
359 int scenes[] = {1,4};
360 for(
int i=0;i<2;i++) {
361 int scenario = scenes[i];
362 const RCP<const Thyra::LinearOpBase<double> > M = buildBDGOperator(scenario,G,G,4.5,out);
368 const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
371 TEST_ASSERT(BtDB_transformer != null);
373 const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
375 BtDB_transformer->transform( *M, M_explicit.ptr() );
377 out <<
"\nM_explicit = " << *M_explicit;
383 LinearOpTester<double> M_explicit_tester;
384 M_explicit_tester.show_all_tests(
true);;
386 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
387 if (!result) success =
false;
399 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\'";
400 out <<
" and from the file \'"<<matrixFile2<<
"\' ...\n";
403 Epetra_MpiComm comm(MPI_COMM_WORLD);
405 Epetra_SerialComm comm;
407 RCP<Epetra_CrsMatrix> epetra_B;
408 RCP<Epetra_CrsMatrix> epetra_G;
409 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
410 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
416 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
417 const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
424 for(
int scenario=1;scenario<=4;scenario++) {
425 RCP<const Thyra::LinearOpBase<double> > M;
426 if(scenario==1 || scenario==3)
427 M = buildBGOperator(scenario,B,G,out);
429 M = buildBGOperator(scenario,G,B,out);
431 M = buildBGOperator(scenario,G,G,out);
437 const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
440 TEST_ASSERT(BtDB_transformer != null);
442 const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
444 BtDB_transformer->transform( *M, M_explicit.ptr() );
446 out <<
"\nM_explicit = " << *M_explicit;
452 LinearOpTester<double> M_explicit_tester;
453 M_explicit_tester.show_all_tests(
true);;
455 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
456 if (!result) success =
false;
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)