#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_DefaultSpmdVectorSpace.hpp"
#include "Thyra_VectorSpaceTester.hpp"
#include "Thyra_TestingTools.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_DefaultComm.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#include "Teuchos_VerboseObject.hpp"
template <class Scalar>
bool run_product_space_tests(
const int n
,const int numBlocks
,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType &tol
,const bool showAllTests
,const bool dumpAll
,Teuchos::FancyOStream *out_arg
)
{
using Thyra::relErr;
using Teuchos::OSTab;
typedef Teuchos::ScalarTraits<Scalar> ST;
typedef typename ST::magnitudeType ScalarMag;
Teuchos::RefCountPtr<Teuchos::FancyOStream>
out = rcp(new Teuchos::FancyOStream(rcp(out_arg,false)));
if(out.get()) *out << "\n*** Entering run_product_space_tests<"<<ST::name()<<">(...) ...\n";
bool success = true, result;
Scalar sresult1, sresult2;
Thyra::VectorSpaceTester<Scalar> vectorSpaceTester;
vectorSpaceTester.warning_tol(ScalarMag(0.1)*tol);
vectorSpaceTester.error_tol(tol);
vectorSpaceTester.show_all_tests(showAllTests);
vectorSpaceTester.dump_all(dumpAll);
Teuchos::Array<Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > >
vecSpaces(numBlocks);
const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >
comm = Teuchos::DefaultComm<Thyra::Index>::getComm();
const int numProcs = size(*comm);
Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> >
spaceBlock = Teuchos::rcp(new Thyra::DefaultSpmdVectorSpace<Scalar>(comm,n,-1));
for( int i = 0; i < numBlocks; ++i )
vecSpaces[i] = spaceBlock;
if(out.get()) *out << "\nA) Performing basic tests on product vectors with SPMD constituent vectors ...\n";
if(out.get()) *out << "\nCreating a product space ps with numBlocks="<<numBlocks<<" and n="<<n<<"vector elements per block ...\n";
Thyra::DefaultProductVectorSpace<Scalar> ps(numBlocks,&vecSpaces[0]);
if(out.get()) *out << "\nps.numBlocks()=";
result = ps.numBlocks() == numBlocks;
if(!result) success = false;
if(out.get()) *out
<< ps.numBlocks() << " == numBlocks=" << numBlocks
<< " : " << ( result ? "passed" : "failed" ) << std::endl;
if(out.get()) *out << "\nTesting the product space ps ...\n";
if(out.get()) *out << "\nps.dim()=";
result = ps.dim() == numProcs*n*numBlocks;
if(!result) success = false;
if(out.get()) *out
<< ps.dim() << " == numProcs*n*numBlocks=" << numProcs*n*numBlocks
<< " : " << ( result ? "passed" : "failed" ) << std::endl;
if(out.get()) *out << "\nTesting the VectorSpaceBase interface of ps ...\n";
result = vectorSpaceTester.check(ps,out.get());
if(!result) success = false;
if(out.get()) *out << "\nB) Testing a nested product space of product vector spaces called pps ...\n";
Teuchos::Array<Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > >
blockVecSpaces(numBlocks);
for( int i = 0; i < numBlocks; ++i )
blockVecSpaces[i] = Teuchos::rcp(&ps,false);
Thyra::DefaultProductVectorSpace<Scalar> pps(numBlocks,&blockVecSpaces[0]);
if(out.get()) *out << "\nTesting the VectorSpaceBase interface of pps ...\n";
result = vectorSpaceTester.check(pps,out.get());
if(!result) success = false;
if(numProcs==1) {
if(out.get()) *out
<< "\nC) Test the compatibility of serial vectors and product vectors with serial blocks."
<< "\n These tests demonstrate the principle of how all in-core vectors are compatible ...\n";
const Scalar
one = ST::one(),
two = Scalar(2)*one,
three = Scalar(3)*one;
if(out.get()) *out << "\nCreating a serial vector space ss with numBlocks*n=" << numBlocks*n << " vector elements ...\n";
Thyra::DefaultSpmdVectorSpace<Scalar> ss(numBlocks*n);
if(out.get()) *out << "\nTesting the serial space ss ...\n";
result = vectorSpaceTester.check(ss,out.get());
if(!result) success = false;
if(out.get()) *out << "\nCreating product vectors; pv1, pv2 ...\n";
Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
pv1 = createMember(ps),
pv2 = createMember(ps);
if(out.get()) *out << "\nassign(&pv1,2.0) ...\n";
Thyra::assign( &*pv1, two );
if(out.get()) *out << "\nassign(&pv1,3.0) ...\n";
Thyra::assign( &*pv2, three );
if(out.get()) *out << "\nCreating serial vectors; sv1, sv2 ...\n";
Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> >
sv1 = createMember(ss),
sv2 = createMember(ss);
if(out.get()) *out << "\nassign(&sv1,*pv1) ...\n";
Thyra::assign( &*sv1, *pv1 );
if(out.get()) *out << "\nsum(sv1)=";
sresult1 = Thyra::sum(*sv1);
sresult2 = two*Scalar(ps.dim());
result = ( ST::magnitude( Thyra::relErr( sresult1, sresult2 ) )
< ST::magnitude( tol ) );
if(!result) success = false;
if(out.get()) *out
<< sresult1 << " == 2*ps.dim()=" << sresult2
<< " : " << ( result ? "passed" : "failed" ) << std::endl;
if(out.get() && dumpAll) *out
<< "\nsv1 =\n" << *sv1;
if(out.get()) *out << "\nassign(&pv2,*sv1) ...\n";
Thyra::assign( &*pv2, *sv1 );
if(out.get()) *out << "\nsum(pv2)=";
sresult1 = Thyra::sum(*pv2);
sresult2 = two*Scalar(ps.dim());
result = ( ST::magnitude( Thyra::relErr( sresult1, sresult2 ) )
< ST::magnitude( tol ) );
if(!result) success = false;
if(out.get()) *out
<< sresult1 << " == 2*ps.dim()=" << sresult2
<< " : " << ( result ? "passed" : "failed" ) << std::endl;
if(out.get() && dumpAll) *out
<< "\npv2 =\n" << *pv2;
}
if(out.get()) *out
<< "\n*** Leaving run_product_space_tests<"<<ST::name()<<">(...) ...\n";
return success;
}
int main( int argc, char* argv[] ) {
using Teuchos::CommandLineProcessor;
bool success = true;
bool verbose = true;
Teuchos::GlobalMPISession mpiSession(&argc,&argv);
Teuchos::RefCountPtr<Teuchos::FancyOStream>
out = Teuchos::VerboseObjectBase::getDefaultOStream();
try {
int n = 4;
int numBlocks = 4;
bool showAllTests = true;
bool dumpAll = false;
CommandLineProcessor clp;
clp.throwExceptions(false);
clp.addOutputSetupOptions(true);
clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not." );
clp.setOption( "n", &n, "Number of elements in each constituent vector." );
clp.setOption( "num-blocks", &numBlocks, "blocks to create." );
clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );
clp.setOption( "show-all-tests", "no-show-all-tests", &showAllTests, "Determines if all tests are printed or not." );
CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
if( !run_product_space_tests<float>(n,numBlocks,float(1e-5),showAllTests,dumpAll,verbose?&*out:NULL) ) success = false;
if( !run_product_space_tests<double>(n,numBlocks,double(1e-13),showAllTests,dumpAll,verbose?&*out:NULL) ) success = false;
#if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
if( !run_product_space_tests<std::complex<float> >(n,numBlocks,float(1e-5),showAllTests,dumpAll,verbose?&*out:NULL) ) success = false;
if( !run_product_space_tests<std::complex<double> >(n,numBlocks,double(1e-13),showAllTests,dumpAll,verbose?&*out:NULL) ) success = false;
#endif
#ifdef HAVE_TEUCHOS_GNU_MP
#endif
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
if(verbose) {
if(success)
*out << "\nAll of the tests seem to have run successfully!\n";
else
*out << "\nOh no! at least one of the test failed!\n";
}
return success ? 0 : 1;
}
