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

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

 ./thyra/test/test_epetra_adapters.cpp 

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;
#endif
    CommandLineProcessor  clp;
    clp.throwExceptions(false);
    clp.addOutputSetupOptions(true);
    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." );
#endif
    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;
    }
#endif

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:

    RefCountPtr<const Epetra_Comm> epetra_comm;
    RefCountPtr<const Epetra_Map> epetra_map;
    RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > epetra_vs;
    RefCountPtr<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>(
          Thyra::create_Comm(epetra_comm)
          ,local_dim,-1
          )
        );
    }
    else {
#endif
      //
      // 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:

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

    RefCountPtr<Thyra::MultiVectorBase<Scalar> >
      eV1 = createMembers(epetra_vs,num_mv_cols),
      eV2 = createMembers(epetra_vs,num_mv_cols);
    RefCountPtr<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:

    Scalar
      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 step 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:

    linearOpTester.set_all_warning_tol(max_rel_warn);
    linearOpTester.set_all_error_tol(max_rel_err);
    linearOpTester.dump_all(dumpAll);
    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:

    RefCountPtr<Epetra_Operator>  epetra_op;

    if(1) {
      // Create a diagonal matrix with scalar on the diagonal
      RefCountPtr<Epetra_CrsMatrix>
        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
        epetra_mat->InsertGlobalValues(
          offset + k + IB     // GlobalRow
          ,1                  // NumEntries
          ,values             // Values
          ,indices            // Indices
          );
      }
      epetra_mat->FillComplete();
      epetra_op = epetra_mat;
    } // end epetra_op

We then wrap this Epetra operator object as:

    RefCountPtr<const Thyra::LinearOpBase<Scalar> >
      Op = rcp(new Thyra::EpetraLinearOp(epetra_op));

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:

    RefCountPtr<Thyra::VectorBase<Scalar> >
      ey  = createMember(epetra_vs);
    RefCountPtr<Thyra::MultiVectorBase<Scalar> >
      eY  = createMembers(epetra_vs,num_mv_cols);
    RefCountPtr<Thyra::VectorBase<Scalar> >
      ney = createMember(non_epetra_vs);
    RefCountPtr<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 Thu Sep 18 12:37:49 2008 for Epetra to Thyra Adapter Software by doxygen 1.3.9.1