Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraOperatorWrapper_UnitTests.cpp
Go to the documentation of this file.
00001 #include "Thyra_EpetraOperatorWrapper.hpp"
00002 #include "Thyra_VectorStdOps.hpp"
00003 #include "Thyra_EpetraThyraWrappers.hpp"
00004 #include "Thyra_EpetraLinearOp.hpp"
00005 #include "Thyra_DefaultBlockedLinearOp.hpp"
00006 #include "Thyra_ProductVectorBase.hpp"
00007 #include "Thyra_SpmdVectorSpaceBase.hpp"
00008 #include "Thyra_DetachedSpmdVectorView.hpp"
00009 #include "Thyra_TestingTools.hpp"
00010 #include "Thyra_LinearOpTester.hpp"
00011 #include "Trilinos_Util_CrsMatrixGallery.h"
00012 #include "Teuchos_GlobalMPISession.hpp"
00013 #include "Teuchos_VerboseObject.hpp"
00014 #include "Teuchos_XMLParameterListHelpers.hpp"
00015 #include "Teuchos_CommandLineProcessor.hpp"
00016 #include "Teuchos_StandardCatchMacros.hpp"
00017 
00018 #ifdef HAVE_MPI
00019 #  include "Epetra_MpiComm.h"
00020 #else
00021 #  include "Epetra_SerialComm.h"
00022 #endif
00023 #include "Epetra_Vector.h"
00024 #include "Epetra_CrsMatrix.h"
00025 
00026 #include "Teuchos_UnitTestHarness.hpp"
00027 
00028 
00029 namespace Thyra {
00030 
00031 
00032 using Teuchos::null;
00033 using Teuchos::RCP;
00034 using Teuchos::rcp;
00035 using Teuchos::rcp_dynamic_cast;
00036 using Teuchos::inOutArg;
00037 using Teuchos::as;
00038 
00039 
00040 TEUCHOS_UNIT_TEST( EpetraOperatorWrapper, basic )
00041 {
00042 
00043 #ifdef HAVE_MPI
00044    Epetra_MpiComm comm(MPI_COMM_WORLD);
00045 #else
00046    Epetra_SerialComm comm;
00047 #endif
00048 
00049    out << "\nRunning on " << comm.NumProc() << " processors\n";
00050 
00051    int nx = 39; // essentially random values
00052    int ny = 53;
00053 
00054    out << "Using Trilinos_Util to create test matrices\n";
00055 
00056    // create some big blocks to play with
00057    Trilinos_Util::CrsMatrixGallery FGallery("recirc_2d",comm);
00058    FGallery.Set("nx",nx);
00059    FGallery.Set("ny",ny);
00060    RCP<Epetra_CrsMatrix> F = rcp(FGallery.GetMatrix(),false);
00061 
00062    Trilinos_Util::CrsMatrixGallery CGallery("laplace_2d",comm);
00063    CGallery.Set("nx",nx);
00064    CGallery.Set("ny",ny);
00065    RCP<Epetra_CrsMatrix> C = rcp(CGallery.GetMatrix(),false);
00066 
00067    Trilinos_Util::CrsMatrixGallery BGallery("diag",comm);
00068    BGallery.Set("nx",nx*ny);
00069    BGallery.Set("a",5.0);
00070    RCP<Epetra_CrsMatrix> B = rcp(BGallery.GetMatrix(),false);
00071 
00072    Trilinos_Util::CrsMatrixGallery BtGallery("diag",comm);
00073    BtGallery.Set("nx",nx*ny);
00074    BtGallery.Set("a",3.0);
00075    RCP<Epetra_CrsMatrix> Bt = rcp(BtGallery.GetMatrix(),false);
00076 
00077    // load'em up in a thyra operator
00078    out << "Building block2x2 Thyra matrix ... wrapping in EpetraOperatorWrapper\n";
00079    const RCP<const LinearOpBase<double> > A =
00080      Thyra::block2x2<double>(
00081        Thyra::epetraLinearOp(F),
00082        Thyra::epetraLinearOp(Bt),
00083        Thyra::epetraLinearOp(B),
00084        Thyra::epetraLinearOp(C),
00085        "A"
00086        );
00087 
00088    const RCP<Thyra::EpetraOperatorWrapper> epetra_A =
00089      rcp(new Thyra::EpetraOperatorWrapper(A));
00090 
00091    // begin the tests!
00092    const Epetra_Map & rangeMap  = epetra_A->OperatorRangeMap();
00093    const Epetra_Map & domainMap = epetra_A->OperatorDomainMap();
00094 
00095    // check to see that the number of global elements is correct
00096    TEST_EQUALITY(rangeMap.NumGlobalElements(), 2*nx*ny);
00097    TEST_EQUALITY(domainMap.NumGlobalElements(), 2*nx*ny);
00098 
00099    // largest global ID should be one less then the # of elements
00100    TEST_EQUALITY(rangeMap.NumGlobalElements()-1, rangeMap.MaxAllGID());
00101    TEST_EQUALITY(domainMap.NumGlobalElements()-1, domainMap.MaxAllGID());
00102 
00103    // create a vector to test: copyThyraIntoEpetra
00104    {
00105       const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
00106       Thyra::randomize(-100.0, 100.0, tv.ptr());
00107       const RCP<const VectorBase<double> > tv_0 =
00108         Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
00109       const RCP<const VectorBase<double> > tv_1 =
00110         Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
00111       const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
00112       const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
00113 
00114       int off_0 = vv_0.globalOffset();
00115       int off_1 = vv_1.globalOffset();
00116       
00117       // create its Epetra counter part
00118       Epetra_Vector ev(epetra_A->OperatorDomainMap());
00119       epetra_A->copyThyraIntoEpetra(*tv, ev);
00120 
00121       // compare handle_tv to ev!
00122       TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
00123       const int numMyElements = domainMap.NumMyElements();
00124       double tval = 0.0;
00125       for(int i=0; i < numMyElements; i++) {
00126          int gid = domainMap.GID(i);
00127          if(gid<nx*ny)
00128             tval = vv_0[gid-off_0];
00129          else
00130             tval = vv_1[gid-off_1-nx*ny];
00131          TEST_EQUALITY(ev[i], tval);
00132       }
00133    }
00134 
00135    // create a vector to test: copyEpetraIntoThyra
00136    {
00137       // create an Epetra vector
00138      Epetra_Vector ev(epetra_A->OperatorDomainMap());
00139      ev.Random();
00140 
00141       // create its thyra counterpart
00142       const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
00143       const RCP<const VectorBase<double> > tv_0 =
00144         Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
00145       const RCP<const VectorBase<double> > tv_1 =
00146         Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
00147       const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
00148       const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
00149 
00150       int off_0 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
00151         tv_0->space())->localOffset();
00152       int off_1 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
00153         tv_1->space())->localOffset();
00154 
00155       epetra_A->copyEpetraIntoThyra(ev, tv.ptr());
00156    
00157       // compare tv to ev!
00158       TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
00159       int numMyElements = domainMap.NumMyElements();
00160       double tval = 0.0;
00161       for(int i=0;i<numMyElements;i++) {
00162          int gid = domainMap.GID(i);
00163          if(gid<nx*ny)
00164             tval = vv_0[gid-off_0];
00165          else
00166             tval = vv_1[gid-off_1-nx*ny];
00167          TEST_EQUALITY(ev[i], tval);
00168       }
00169    }
00170 
00171    // Test using Thyra::LinearOpTester
00172    const RCP<const LinearOpBase<double> > thyraEpetraOp = epetraLinearOp(epetra_A);
00173    LinearOpTester<double> linearOpTester;
00174    linearOpTester.show_all_tests(true);
00175    const bool checkResult = linearOpTester.check(*thyraEpetraOp, inOutArg(out));
00176    TEST_ASSERT(checkResult);
00177 
00178 }
00179 
00180 
00181 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines