Testing Program for Thyra/Epetra Adapters
[Assorted Epetra to Thyra Operator/Vector Adapter Example Code]

Collaboration diagram for Testing Program for Thyra/Epetra Adapters:

Here we show the testing program that is used to perform regression tests on the Epetra adapter subclasses and tests a number of other classes as well. More...
This testing program is contained in the source file


where ./ is the base source directory for epetra.

This testing program will be shown in its entirety below but first we will point out some of the more significant parts of the program and highlight some interesting functionality.

First, the testing program reads in a number of options as shown below:

    int     local_dim            = 1000;
    int     num_mv_cols          = 4;
    double  max_rel_err          = 1e-13;
    double  max_rel_warn         = 1e-15;
    double  scalar               = 1.5;
    double  max_flop_rate        = 2.0e8;
#ifdef RTOp_USE_MPI
    bool    useMPI               = true;
    CommandLineProcessor  clp;
    clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
    clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
    clp.setOption( "local-dim", &local_dim, "Number of vector elements per process." );
    clp.setOption( "num-mv-cols", &num_mv_cols, "Number columns in each multi-vector (>=4)." );
    clp.setOption( "max-rel-err-tol", &max_rel_err, "Maximum relative error tolerance for tests." );
    clp.setOption( "max-rel-warn-tol", &max_rel_warn, "Maximum relative warning tolerance for tests." );
    clp.setOption( "scalar", &scalar, "A scalar used in all computations." );
    clp.setOption( "max-flop-rate", &max_flop_rate, "Approx flop rate used for loop timing." );
#ifdef RTOp_USE_MPI
    clp.setOption( "use-mpi", "no-use-mpi", &useMPI, "Actually use MPI or just run independent serial programs." );
    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

Next, some MPI functions are called to get the context of the parallel program (if built with MPI support and RTOp_USE_MPI is defined):

#ifdef RTOp_USE_MPI
    MPI_Comm mpiComm;
    int numProc;
    if(useMPI) {
      mpiComm = MPI_COMM_WORLD;
      MPI_Comm_size( mpiComm, &numProc );
      MPI_Comm_rank( mpiComm, &procRank );
    else {
      mpiComm = MPI_COMM_NULL;
      numProc = 1;
      procRank = 0;

This program then creates two different (parallel or serial) vector spaces where one is of type Thyra::EpetraVectorSpace and the other is some non-Epetra vector space:

    RCP<const Epetra_Comm> epetra_comm;
    RCP<const Epetra_Map> epetra_map;
    RCP<const Thyra::VectorSpaceBase<Scalar> > epetra_vs;
    RCP<const Thyra::VectorSpaceBase<Scalar> > non_epetra_vs;
#ifdef RTOp_USE_MPI
    if(useMPI) {
      // Create parallel vector spaces with compatible maps
      // Epetra vector space
      if(verbose) *out << "\nCreating vector space using Epetra_MpiComm ...\n";
      epetra_comm = rcp(new Epetra_MpiComm(mpiComm));
      epetra_map = rcp(new Epetra_Map(-1,local_dim,0,*epetra_comm));
      epetra_vs = Thyra::create_VectorSpace(epetra_map);
      // Non-Epetra vector space
      if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
      non_epetra_vs = rcp(
        new Thyra::DefaultSpmdVectorSpace<Scalar>(
    else {
      // Create serial vector spaces (i.e. VectorSpaceBase::isInCore()==true)
      // Epetra vector space
      if(verbose) *out << "\nCreating vector space using Epetra_SerialComm ...\n";
      epetra_comm = rcp(new Epetra_SerialComm);
      epetra_map = rcp(new Epetra_LocalMap(local_dim,0,*epetra_comm));
      epetra_vs = Thyra::create_VectorSpace(epetra_map);
      // Non-Epetra vector space
      if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
      non_epetra_vs = rcp(new Thyra::DefaultSpmdVectorSpace<Scalar>(local_dim));
#ifdef RTOp_USE_MPI
#endif // end create vector spacdes [Doxygen looks for this!]

Given these vector spaces a couple Epetra and non-Epetra vectors and multi-vectors are created:

    RCP<Thyra::VectorBase<Scalar> >
      ev1 = createMember(epetra_vs),
      ev2 = createMember(epetra_vs);
    RCP<Thyra::VectorBase<Scalar> >
      nev1 = createMember(non_epetra_vs),
      nev2 = createMember(non_epetra_vs);

    RCP<Thyra::MultiVectorBase<Scalar> >
      eV1 = createMembers(epetra_vs,num_mv_cols),
      eV2 = createMembers(epetra_vs,num_mv_cols);
    RCP<Thyra::MultiVectorBase<Scalar> >
      neV1 = createMembers(non_epetra_vs,num_mv_cols),
      neV2 = createMembers(non_epetra_vs,num_mv_cols);

Most of the operations performed in this testing program are fairly mundane and we show a few of these simpler examples before going onto some more interesting ones.

These vectors and multi-vectors are then initialized using the Thyra::assign() wrapper functions:

    Thyra::assign( &*ev1, 0.0 );
    Thyra::assign( &*ev2, scalar );
    Thyra::assign( &*nev1, 0.0 );
    Thyra::assign( &*nev2, scalar );
    Thyra::assign( &*eV1, 0.0 );
    Thyra::assign( &*eV2, scalar );
    Thyra::assign( &*neV1, 0.0 );
    Thyra::assign( &*neV2, scalar );

The one norms of these vectors and multi-vectors are are then computed using various overloads of Thyra::norm_1() as:

      ev1_nrm = Thyra::norm_1(*ev1),
      ev2_nrm = Thyra::norm_1(*ev2),
      eV1_nrm = Thyra::norm_1(*eV1),
      eV2_nrm = Thyra::norm_1(*eV2),
      nev1_nrm = Thyra::norm_1(*nev1),
      nev2_nrm = Thyra::norm_1(*nev2),
      neV1_nrm = Thyra::norm_1(*neV1),
      neV2_nrm = Thyra::norm_1(*neV2);

The first set of more interesting operations are mixing objects of different concrete type in vector transformation operations. The following code shows assigning vectors and multi-vectors of various types

     Thyra::assign( &*ev1, *ev2 );
     Thyra::assign( &*eV1, *eV2 );
     Thyra::assign( &*ev1, *nev2 );
     Thyra::assign( &*nev1, *ev2 );
     Thyra::assign( &*nev1, *nev2 );
     Thyra::assign( &*eV1, *neV2 );
     Thyra::assign( &*neV1, *eV2 );

What is significant about some of the above function calls is that they mix and match objects with different concrete types. For example, a non-Epetra vector is assigned to an Epetra in the line:

Thyra::assign( &*ev1, *nev2 );

and an Epetra vector is assigned to a non-Epetra in the line:

Thyra::assign( &*nev1, *ev2 );

Similarly, non-Epetra and Epetra multi-vectors are assigned to and from in the calls:

Thyra::assign( &*eV1, *neV2 );
Thyra::assign( &*neV1, *eV2 );

Developer's that really want to understand how this type of interoperability between different concrete types works are cautiously encouraged to look at the implementations of the functions Thyra::SerialVectorBase::applyOp() and Thyra::SerialMultiVectorBase::applyOp() in a serial build and Thyra::MPIVectorBase::applyOp() and Thyra::MPIMultiVectorBase::applyOp().

The Thyra::LinearOpBase interface that all vectors and multi-vectors support is tested by the testing class Thyra::LinearOpTester which is created in the line:

    Thyra::LinearOpTester<Scalar> linearOpTester;

This testing class is then called for all of the vector and multi-vector objects is the following lines:

    result = linearOpTester.check(*ev1,verbose?&*out:NULL);

Now we show how an Epetra_Operator object is wrapped by a Thyra::EpetraLinearOp object to give the Epetra object a Thyra skin.

First we create a simple diagonal Epetra operator in the following lines of code:

    RCP<Epetra_Operator>  epetra_op;

      // Create a diagonal matrix with scalar on the diagonal
        epetra_mat = rcp(new Epetra_CrsMatrix(::Copy,*epetra_map,1));
      Scalar values[1] = { scalar };
      int indices[1];
      const int IB = epetra_map->IndexBase(), offset = procRank*local_dim;
      for( int k = 0; k < local_dim; ++k ) {
        indices[0] = offset + k + IB;  // global column
          offset + k + IB     // GlobalRow
          ,1                  // NumEntries
          ,values             // Values
          ,indices            // Indices
      epetra_op = epetra_mat;
    } // end epetra_op

We then wrap this Epetra operator object as:

    RCP<const Thyra::LinearOpBase<Scalar> >
      Op = Thyra::epetraLinearOp(epetra_op);

    if(verbose && dumpAll) *out << "\nOp=\n" << *Op;

    if(verbose) *out << "\n*** (B.6b) Going through partial then full initialization of EpetraLinearOp ...\n";

This object is then tested in the line:

    result = linearOpTester.check(*Op,verbose?&*out:NULL);

Before the next set of tests we create some Epetra and non-Epetra working vectors and mulit-vectors in the lines:

    RCP<Thyra::VectorBase<Scalar> >
      ey  = createMember(epetra_vs);
    RCP<Thyra::MultiVectorBase<Scalar> >
      eY  = createMembers(epetra_vs,num_mv_cols);
    RCP<Thyra::VectorBase<Scalar> >
      ney = createMember(non_epetra_vs);
    RCP<Thyra::MultiVectorBase<Scalar> >
      neY = createMembers(non_epetra_vs,num_mv_cols);

In the next set of lines of code shows how different Epetra and non-Epetra vectors and multi-vectors can be mixed and matched with with this Thyra::EpetraLinearOp object:

    apply( *Op, NOTRANS, *ev1, &*ey, 2.0 );
    apply( *Op, NOTRANS, *eV1, &*eY, 2.0 );
    apply( *Op, NOTRANS, *ev1, &*ney, 2.0 );
    apply( *Op, NOTRANS, *eV1, &*neY, 2.0 );
    apply( *Op, NOTRANS, *nev1, &*ey, 2.0 );
    apply( *Op, NOTRANS, *neV1, &*eY, 2.0 );
    apply( *Op, NOTRANS, *nev1, &*ney, 2.0 );
    apply( *Op, NOTRANS, static_cast<const Thyra::MultiVectorBase<Scalar>&>(*nev1), static_cast<Thyra::MultiVectorBase<Scalar>*>(&*ney), 2.0 );
    apply( *Op, NOTRANS, *neV1, &*neY, 2.0 );

Note that many of the above function calls involve using non-Epetra vector and multi-vector objects with the underlying Epetra_Operator object. This magic happens through the use of the utility function Thyra::get_Epetra_MultiVector() by the overridden function Thyra::EpetraLinearOp::euclideanApply().

ToDo: Finish this discussion.

To see the full listing of this click: test_epetra_adpaters.cpp .

Generated on Wed May 12 21:42:56 2010 for Epetra to Thyra Adapter Software by  doxygen 1.4.7