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" 55 #include "EpetraExt_readEpetraLinearSystem.h" 56 #include "Teuchos_GlobalMPISession.hpp" 57 #include "Teuchos_VerboseObject.hpp" 58 #include "Teuchos_XMLParameterListHelpers.hpp" 59 #include "Teuchos_CommandLineProcessor.hpp" 60 #include "Teuchos_StandardCatchMacros.hpp" 67 # include "Epetra_MpiComm.h" 69 # include "Epetra_SerialComm.h" 72 #include "Teuchos_UnitTestHarness.hpp" 81 using Thyra::epetraExtDiagScalingTransformer;
82 using Thyra::VectorBase;
83 using Thyra::LinearOpBase;
84 using Thyra::createMember;
85 using Thyra::LinearOpTester;
87 using Thyra::multiply;
88 using Thyra::diagonal;
91 std::string matrixFile =
"";
96 Teuchos::UnitTestRepository::getCLP().setOption(
97 "matrix-file", &matrixFile,
98 "Defines the Epetra_CrsMatrix to read in." );
102 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
103 buildADOperator(
int scenario,
const Teuchos::RCP<
const Thyra::LinearOpBase<double> > & A,
104 const double vecScale, std::ostream & out)
110 RCP<const LinearOpBase<double> > M;
111 RCP<VectorBase<double> > d;
113 d = createMember(A->domain(),
"d");
115 d = createMember(A->range(),
"d");
116 V_S( d.ptr(), vecScale );
122 M = multiply( scale(scalea,A), scale(scaled,diagonal(d)),
"M" );
123 out <<
" Testing A*D" << std::endl;
126 M = multiply( scale(scaled,diagonal(d)), scale(scalea,A),
"M" );
127 out <<
" Testing D*A" << std::endl;
130 M = multiply( A, scale(scaled,diagonal(d)),
"M" );
131 out <<
" Testing adj(A)*D" << std::endl;
134 M = multiply( scale(scaled,diagonal(d)), A,
"M" );
135 out <<
" Testing D*adj(A)" << std::endl;
138 TEUCHOS_ASSERT(
false);
143 out <<
"\nM = " << *M;
150 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\'\n";
153 Epetra_MpiComm comm(MPI_COMM_WORLD);
155 Epetra_SerialComm comm;
158 RCP<Epetra_CrsMatrix> epetra_A;
159 RCP<Epetra_CrsMatrix> epetra_B;
160 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
161 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
162 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
163 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
165 RCP<VectorBase<double> > d = createMember(B->domain(),
"d");
168 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
171 TEST_ASSERT(BD_transformer != null);
173 RCP<const LinearOpBase<double> > M;
176 TEST_ASSERT(not BD_transformer->isCompatible(*M));
178 M = multiply(A,scale(3.9,diagonal(d)));
179 TEST_ASSERT(BD_transformer->isCompatible(*M));
181 M = multiply(scale(3.0,diagonal(d)),scale(9.0,B));
182 TEST_ASSERT(BD_transformer->isCompatible(*M));
184 M = multiply(B,scale(3.0,diagonal(d)),scale(9.0,B));
185 TEST_ASSERT(not BD_transformer->isCompatible(*M));
194 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\'\n";
197 Epetra_MpiComm comm(MPI_COMM_WORLD);
199 Epetra_SerialComm comm;
201 RCP<Epetra_CrsMatrix> epetra_A;
202 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
208 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
210 RCP<const Thyra::LinearOpBase<double> > M;
213 for(
int scenario=1;scenario<=2;scenario++) {
214 out <<
"**********************************************************" << std::endl;
215 out <<
"**************** Starting Scenario = " << scenario << std::endl;
216 out <<
"**********************************************************" << std::endl;
221 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
224 TEST_ASSERT(BD_transformer != null);
230 const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
231 TEST_ASSERT(M_explicit != null);
234 M = buildADOperator(scenario,A,4.5,out);
235 BD_transformer->transform( *M, M_explicit.ptr() );
238 M = buildADOperator(scenario,A,7.5,out);
239 BD_transformer->transform( *M, M_explicit.ptr() );
243 TEST_EQUALITY(eOp_explicit0,eOp_explicit1);
245 out <<
"\nM_explicit = " << *M_explicit;
251 LinearOpTester<double> M_explicit_tester;
252 M_explicit_tester.show_all_tests(
true);;
254 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
255 if (!result) success =
false;
266 out <<
"\nReading linear system in Epetra format from the file \'"<<matrixFile<<
"\'\n";
269 Epetra_MpiComm comm(MPI_COMM_WORLD);
271 Epetra_SerialComm comm;
273 RCP<Epetra_CrsMatrix> epetra_A;
274 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
280 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
287 for(
int scenario=1;scenario<=4;scenario++) {
288 out <<
"**********************************************************" << std::endl;
289 out <<
"**************** Starting Scenario = " << scenario << std::endl;
290 out <<
"**********************************************************" << std::endl;
292 RCP<const Thyra::LinearOpBase<double> > M = buildADOperator(scenario,A,4.5,out);
297 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
300 TEST_ASSERT(BD_transformer != null);
302 BD_transformer->isCompatible(*M);
308 const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
309 BD_transformer->transform( *M, M_explicit.ptr() );
311 out <<
"\nM_explicit = " << *M_explicit;
317 LinearOpTester<double> M_explicit_tester;
318 M_explicit_tester.show_all_tests(
true);
320 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
321 if (!result) success =
false;
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.