Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraOperatorWrapper_UnitTests.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include "Thyra_EpetraOperatorWrapper.hpp"
00045 #include "Thyra_VectorStdOps.hpp"
00046 #include "Thyra_EpetraThyraWrappers.hpp"
00047 #include "Thyra_EpetraLinearOp.hpp"
00048 #include "Thyra_DefaultBlockedLinearOp.hpp"
00049 #include "Thyra_ProductVectorBase.hpp"
00050 #include "Thyra_SpmdVectorSpaceBase.hpp"
00051 #include "Thyra_DetachedSpmdVectorView.hpp"
00052 #include "Thyra_TestingTools.hpp"
00053 #include "Thyra_LinearOpTester.hpp"
00054 #include "Trilinos_Util_CrsMatrixGallery.h"
00055 #include "Teuchos_GlobalMPISession.hpp"
00056 #include "Teuchos_VerboseObject.hpp"
00057 #include "Teuchos_XMLParameterListHelpers.hpp"
00058 #include "Teuchos_CommandLineProcessor.hpp"
00059 #include "Teuchos_StandardCatchMacros.hpp"
00060 
00061 #ifdef HAVE_MPI
00062 #  include "Epetra_MpiComm.h"
00063 #else
00064 #  include "Epetra_SerialComm.h"
00065 #endif
00066 #include "Epetra_Vector.h"
00067 #include "Epetra_CrsMatrix.h"
00068 
00069 #include "Teuchos_UnitTestHarness.hpp"
00070 
00071 
00072 namespace Thyra {
00073 
00074 
00075 using Teuchos::null;
00076 using Teuchos::RCP;
00077 using Teuchos::rcp;
00078 using Teuchos::rcp_dynamic_cast;
00079 using Teuchos::inOutArg;
00080 using Teuchos::as;
00081 
00082 
00083 TEUCHOS_UNIT_TEST( EpetraOperatorWrapper, basic )
00084 {
00085 
00086 #ifdef HAVE_MPI
00087    Epetra_MpiComm comm(MPI_COMM_WORLD);
00088 #else
00089    Epetra_SerialComm comm;
00090 #endif
00091 
00092    out << "\nRunning on " << comm.NumProc() << " processors\n";
00093 
00094    int nx = 39; // essentially random values
00095    int ny = 53;
00096 
00097    out << "Using Trilinos_Util to create test matrices\n";
00098 
00099    // create some big blocks to play with
00100    Trilinos_Util::CrsMatrixGallery FGallery("recirc_2d",comm,false); // CJ TODO FIXME: change for Epetra64
00101    FGallery.Set("nx",nx);
00102    FGallery.Set("ny",ny);
00103    RCP<Epetra_CrsMatrix> F = rcp(FGallery.GetMatrix(),false);
00104 
00105    Trilinos_Util::CrsMatrixGallery CGallery("laplace_2d",comm,false); // CJ TODO FIXME: change for Epetra64
00106    CGallery.Set("nx",nx);
00107    CGallery.Set("ny",ny);
00108    RCP<Epetra_CrsMatrix> C = rcp(CGallery.GetMatrix(),false);
00109 
00110    Trilinos_Util::CrsMatrixGallery BGallery("diag",comm,false); // CJ TODO FIXME: change for Epetra64
00111    BGallery.Set("nx",nx*ny);
00112    BGallery.Set("a",5.0);
00113    RCP<Epetra_CrsMatrix> B = rcp(BGallery.GetMatrix(),false);
00114 
00115    Trilinos_Util::CrsMatrixGallery BtGallery("diag",comm,false); // CJ TODO FIXME: change for Epetra64
00116    BtGallery.Set("nx",nx*ny);
00117    BtGallery.Set("a",3.0);
00118    RCP<Epetra_CrsMatrix> Bt = rcp(BtGallery.GetMatrix(),false);
00119 
00120    // load'em up in a thyra operator
00121    out << "Building block2x2 Thyra matrix ... wrapping in EpetraOperatorWrapper\n";
00122    const RCP<const LinearOpBase<double> > A =
00123      Thyra::block2x2<double>(
00124        Thyra::epetraLinearOp(F),
00125        Thyra::epetraLinearOp(Bt),
00126        Thyra::epetraLinearOp(B),
00127        Thyra::epetraLinearOp(C),
00128        "A"
00129        );
00130 
00131    const RCP<Thyra::EpetraOperatorWrapper> epetra_A =
00132      rcp(new Thyra::EpetraOperatorWrapper(A));
00133 
00134    // begin the tests!
00135    const Epetra_Map & rangeMap  = epetra_A->OperatorRangeMap();
00136    const Epetra_Map & domainMap = epetra_A->OperatorDomainMap();
00137 
00138    // check to see that the number of global elements is correct
00139    TEST_EQUALITY(rangeMap.NumGlobalElements(), 2*nx*ny);
00140    TEST_EQUALITY(domainMap.NumGlobalElements(), 2*nx*ny);
00141 
00142    // largest global ID should be one less then the # of elements
00143    TEST_EQUALITY(rangeMap.NumGlobalElements()-1, rangeMap.MaxAllGID());
00144    TEST_EQUALITY(domainMap.NumGlobalElements()-1, domainMap.MaxAllGID());
00145 
00146    // create a vector to test: copyThyraIntoEpetra
00147    {
00148       const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
00149       Thyra::randomize(-100.0, 100.0, tv.ptr());
00150       const RCP<const VectorBase<double> > tv_0 =
00151         Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
00152       const RCP<const VectorBase<double> > tv_1 =
00153         Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
00154       const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
00155       const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
00156 
00157       int off_0 = vv_0.globalOffset();
00158       int off_1 = vv_1.globalOffset();
00159       
00160       // create its Epetra counter part
00161       Epetra_Vector ev(epetra_A->OperatorDomainMap());
00162       epetra_A->copyThyraIntoEpetra(*tv, ev);
00163 
00164       // compare handle_tv to ev!
00165       TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
00166       const int numMyElements = domainMap.NumMyElements();
00167       double tval = 0.0;
00168       for(int i=0; i < numMyElements; i++) {
00169          int gid = domainMap.GID(i);
00170          if(gid<nx*ny)
00171             tval = vv_0[gid-off_0];
00172          else
00173             tval = vv_1[gid-off_1-nx*ny];
00174          TEST_EQUALITY(ev[i], tval);
00175       }
00176    }
00177 
00178    // create a vector to test: copyEpetraIntoThyra
00179    {
00180       // create an Epetra vector
00181      Epetra_Vector ev(epetra_A->OperatorDomainMap());
00182      ev.Random();
00183 
00184       // create its thyra counterpart
00185       const RCP<VectorBase<double> > tv = Thyra::createMember(A->domain());
00186       const RCP<const VectorBase<double> > tv_0 =
00187         Thyra::productVectorBase<double>(tv)->getVectorBlock(0);
00188       const RCP<const VectorBase<double> > tv_1 =
00189         Thyra::productVectorBase<double>(tv)->getVectorBlock(1);
00190       const Thyra::ConstDetachedSpmdVectorView<double> vv_0(tv_0);
00191       const Thyra::ConstDetachedSpmdVectorView<double> vv_1(tv_1);
00192 
00193       int off_0 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
00194         tv_0->space())->localOffset();
00195       int off_1 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(
00196         tv_1->space())->localOffset();
00197 
00198       epetra_A->copyEpetraIntoThyra(ev, tv.ptr());
00199    
00200       // compare tv to ev!
00201       TEST_EQUALITY(tv->space()->dim(), as<Ordinal>(ev.GlobalLength()));
00202       int numMyElements = domainMap.NumMyElements();
00203       double tval = 0.0;
00204       for(int i=0;i<numMyElements;i++) {
00205          int gid = domainMap.GID(i);
00206          if(gid<nx*ny)
00207             tval = vv_0[gid-off_0];
00208          else
00209             tval = vv_1[gid-off_1-nx*ny];
00210          TEST_EQUALITY(ev[i], tval);
00211       }
00212    }
00213 
00214    // Test using Thyra::LinearOpTester
00215    const RCP<const LinearOpBase<double> > thyraEpetraOp = epetraLinearOp(epetra_A);
00216    LinearOpTester<double> linearOpTester;
00217    linearOpTester.show_all_tests(true);
00218    const bool checkResult = linearOpTester.check(*thyraEpetraOp, inOutArg(out));
00219    TEST_ASSERT(checkResult);
00220 
00221 }
00222 
00223 
00224 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines