Trilinos/Epetra: Linear Algebra Services Package.

Development

Outline

• Introduction
• Overview of Epetra.
  • Epetra Classes
  • Trilinos and Epetra
  • Fortran and C Support
  • Transition support for current Aztec users
  • Extending Epetra functionality
  • Example Codes
• Epetra adapters for Thyra and PyTrilinos
• Browse all of Epetra as a single doxygen collection

Introduction

Epetra provides the fundamental construction routines and services function that are required for serial and parallel linear algebra libraries. Epetra provides the underlying foundation for all Trilinos solvers.

Overview of Epetra.

Epetra Classes

Epetra contains a number of classes. They can be categorized as follows:

• Primary parallel user classes. These are typically the most important classes for most users.
  1. Communicator class: Epetra_Comm - Contains specific information about the parallel machine we are using. Currently supports serial (Epetra_SerialComm), MPI (Epetra_MpiComm) and prototype hybrid MPI/threaded parallel programming models.

  2. Map classes: Epetra_Map, Epetra_LocalMap, Epetra_BlockMap - Contain information used to distribute vectors, matrices and other objects on a parallel (or serial) machine.

  3. Vector class: Epetra_Vector - Real double precision vector class. Supports construction and use of vectors on a parallel machine. Epetra_FEVector is a specialization of Epetra_Vector to support construction of vectors from finite element applications. Epetra_Vector isa Epetra_MultiVector. Thus, any argument that requires an Epetra_MultiVector can also accept an Epetra_Vector.

  4. Multi-vector class: Epetra_MultiVector - Real double precision multi-vector class. Supports construction and use of multi-vectors on a parallel machine. A multi-vector is a collection vectors. It is a generalizaion of a 2D array.

  5. Sparse row graph class: Epetra_CrsGraph - Allows construction of a serial or parallel graph. The graph determines the communication pattern for subsequent matrix objects. All Epetra sparse matrix classes (Epetra_CrsMatrix, Epetra_VbrMatrix, Epetra_FECrsMatrix and Epetra_FEVbrMatrix) have an associated Epetra_CrsGraph. This graph is either constructed implicitly for the user as part of the construction of the matrix, or is passed in to the matrix constructor by the user. In the latter case the graph is either preconstructed by the user or extracted from an existing Epetra sparse matrix object using the Graph() accessor method.

  6. Pure virtual row matrix class: Epetra_RowMatrix - Pure virtual class that specifies interfaces needed to do most of the common operations required by a row matrix. The Epetra_LinearProblem expects the matrix as a Epetra_RowMatrix. The Epetra_CrsMatrix, Epetra_VbrMatrix, Epetra_FECrsMatrix and Epetra_FEVbrMatrix classes implement the Epetra_RowMatrix interface and therefore objects from these classes can be used with Epetra_LinearProblem and with AztecOO. AztecOO also provides the Epetra_MsrMatrix class as an Epetra_RowMatrix adapter that accepts an AZ_MATRIX struct (as defined in Aztec 2.1) and provides Epetra_RowMatrix functionality using the AZ_Matrix object, without copying matrix data. (See the AztecOO home page for details.) Furthermore, any class that implements Epetra_RowMatrix can be used with Epetra_LinearProblem and AztecOO.

  7. Easier-to-implement row matrix class: Epetra_BasicRowMatrix - Epetra_RowMatrix has many pure virtual functions that must be implemented by an adaptor. In many cases these functions can have a reasonable default implementation in terms of other functions, or using basic information that can be passed into a constructor. Epetra_BasicRowMatrix is intended to make implementation of Epetra_RowMatrix easier. All methods are virtual, but almost all methods have a default implementation that is reasonable. For someone who is implementing their first Epetra_RowMatrix adapter, Epetra_BasicRowMatrix is probably a good starting point.

  8. Sparse row matrix class: Epetra_CrsMatrix - Real double precision sparse matrix class. Supports construction and use of row-wise sparse matrices. Epetra_FECrsMatrix is a specialization that supports finite element applications more naturally.

  9. Sparse block row matrix class: Epetra_VbrMatrix - Real double precision block sparse matrix class. Supports construction and use of row-wise block sparse matrices. Epetra_FEVbrMatrix is a specialization that supports finite element applications more naturally.

  10. Jagged diagonal sparse matrix class: Epetra_JadMatrix - Real double precision sparse matrix class for vector processors. Constructs and updates a jagged-diagonal format sparse matrix from an existing Epetra_RowMatrix. The jagged-diagonal format is one of the better sparse matrix data structures for vector and vector-like processors. Matrix elements are organized so that long vector loops can be executed without memory conflicts. This approach is very effective for machines with vector hardware and an interleaved memory system that supports high bandwidth irregular-pattern memory access. Note that this class can take any existing Epetra sparse matrix as input, but does not require a constructed matrix as input. If the input Epetra_RowMatrix can construct elements on-demand and one-row-at-a-time, then no additional copy of the matrix is needed.

  11. Import/Export classes: Epetra_Import and Epetra_Export - Constructed from two Epetra_BlockMap (or Epetra_Map or Epetra_LocalMap). Allows efficient transfer of objects built using one map to a new object with a new map. Supports local and global permutations, overlapping Schwarz operations and many other data movement algorithms.

• Primary serial user classes. These classes provide object oriented interfaces to LAPACK capabilities, providing easy access to the most powerful numerical methods in LAPACK.
  1. General dense matrix/vector classes: Epetra_SerialDenseMatrix, Epetra_SerialDenseVector - Lightweight dense matrix and vector classes used to define matrices, left-hand-sides and right-hand-sides for the serial solver classes. Epetra_IntSerialDenseVector is also useful for construction of Epetra_Maps and advanced parallel data redistribution requirements. Also note that the C++ standard vector class can be a better alternative to these classes. However the STL vector class can be difficult to debug when memory errors are present.

  2. General dense solver class: Epetra_SerialDenseSolver - Provides dense matrix services such as factorizations, solves, QR, SVD, etc., with special attention focused on numerically robust solutions.

  3. Symmetric dense matrix class: Epetra_SerialSymDenseMatrix - Similar to Epetra_SerialDenseMatrix except focused specifically on symmetric matrices.

  4. Symmetric definite dense solver: Epetra_SerialSpdDenseSolver - Similar to Epetra_SerialDenseSolver except focused specifically on symmetric definite systems.

• Utility classes.
  1. Timing class: Epetra_Time - Provides timing functions for the purposes of performance analysis.

  2. Floating point operation class: Epetra_Flops - Provides floating point operations (FLOPS) counting and reporting functions for the purposes of performance analysis. All Epetra computational classes accumulate FLOP counts associated with the this object of the computations.

  3. Distributed directory class: Epetra_Directory - Allows construction of a distributed directory. Once constructed, a directory allows one to access randomly distributed objects in an efficient, scalable manner. This class is intended for support of general Epetra_BlockMap and Epetra_Map objects, but is useful in other settings as well.

  4. BLAS wrapper class: Epetra_BLAS - A ``thin'' layer of C++ code wrapping the basic linear algebra subprograms (BLAS). This class provides a single instance interface between Epetra and the BLAS. In this way we can easily switch BLAS interfaces and manage the C++/Fortran translation differences that exist between different computer systems. This class also provides a very convenient way to templatize the BLAS.

  5. LAPACK wrapper class: Epetra_LAPACK - A ``thin'' layer of C++ code wrapping LAPACK. Like Epetra_BLAS, it provides nice C++ access to LAPACK.

Fortran and C Support

Epetra provides a small set of user-extendable Fortran and C interfaces. These are found in Epetra_C_wrappers and Epetra_Fortran_wrappers. These wrappers are hand-generated, so adding to the list requires hand-coding. The Trilinos team is looking at automated techniques to support Fortran and C, but we have not made any signficant effort in this direction.

Trilinos and Epetra

Epetra can be used as a stand-alone package. However, it also provides the foundation for Trilinos. Trilinos is a collection of solver packages. The first available package is AztecOO, a C++ implementation of the preconditioned iterative solver package Aztec. This particular class can be used to solve a linear system as defined by a Epetra_LinearProblem object that is passed in at construction. Alternatively, one may pass in the matrix, initial guess and right-hand-side as Epetra objects and solves a linear system using preconditioner and solver options set by the user.

Note:

AztecOO supports all the functionality of Aztec except for explicit scaling options. Explicit scaling is done instead by the Epetra_LinearProblem class (which in turn uses the scaling methods in the Epetra_RowMatrix and Epetra_MultiVector classes). Documentation for Aztec can be found in your Trilinos archive in the file Trilinos/doc/aztec/manual.ps

Transition support for current Aztec users

In order to support existing Aztec users, we have developed a few modules that aid in the transition from Aztec to Trilinos/AztecOO. The first module is AZOO_iterate(). This function has an identical interface to AZ_iterate and, for most purpose, these two should be interchangeable. AZOO_iterate() differs from AZ_iterate in that it first converts the Aztec objects to Epetra objects using the module Aztec2Epetra(). It then creates an Epetra_LinearProblem and performs any requested scaling explicitly, and then calls AztecOO. The transition to AZOO_iterate() is meant to be trivial and temporary. We encourage users to customize their own version of AZOO_iterate(), and to eventually build Epetra objects directly.

AztecOO also provides the Epetra_MsrMatrix class as an Epetra_RowMatrix adapter that accepts an AZ_MATRIX struct (as defined in Aztec 2.1) and provides Epetra_RowMatrix functionality using the AZ_MATRIX object, without copying matrix data. If your application already stores matrix data in an AZ_MATRIX struct, then you can directly use AztecOO, ML and Ifpack without reimplementing the matrix construction process using native Epetra sparse matrix classes.

Extending Epetra functionality

Epetra, and therefore Trilinos, is easily extended to support new implementations of basic functionality. Although Epetra provides serial and MPI support, other parallel communication methods (such as Co-array Fortran, UPC, SHMEM, etc.) can be implemented because parallel machine functionality is accessed via the abstract Epetra_Comm interface. Similarly, solver algorithms that need access to the action of a linear operator on a vector get access to this functionality via the Epetra_Operator interface. Any class (including all Epetra matrix classes, Ifpack and ML preconditioners and AztecOO) that implements the Epetra_Operator interface can be used. The final major abstract interface in Epetra is Epetra_RowMatrix. Any solver algorithm that needs access to matrix coefficients can obtain this information via the Epetra_RowMatrix interface. All Epetra sparse matrix classes implement this interface. AztecOO, Ifpack, ML and other packages access matrix coefficients via the Epetra_RowMatrix interface.

By using abstract base classes to access external functionality and then providing a wide range of default implementations. Novel adapters to Epetra_Comm have been developed experimentally, but no production code is available beyond serial and MPI. Epetra_Operator is commonly extended by users who have complex operators that cannot be explicitly constructed but can be implicitly applied. Examples include ``matrix-free'' methods such as directional derivative approximations and implicit inverses via iterative methods. Most preconditioners are easily expressed as Epetra_Operator objects. Epetra_RowMatrix is often extended to provide support for alternative data structures, or to eliminate the need to explicitly store matrix coefficients but instead build them ``on-the-fly'' row-by-row as needed.

Note: Although Epetra_RowMatrix is the base class that is used to access matrix coefficient information, if you want to provide your own implementation of Epetra_RowMatrix, it is probably best to inherit from the Epetra_BasicRowMatrix. Epetra_BasicRowMatrix provides a reasonable implementation of most Epetra_RowMatrix methods and reduces the number of methods you must implement from 39 to between 4 and 6. To see an example of how to make your own Epetra_RowMatrix class using Epetra_BasicRowMatrix as a starting point, see the Epetra_JadMatrix class.

Why not abstract multivector and vector classes?

Although it may seem logical to have abstract multivector and vector classes, the practival value of such classes is not as obvious to the Epetra developers, at least at the Epetra level of functionality. Serious efforts in abstract multivectors/vectors is very important to Trilinos, however our efforts in this area have been focused on the Thyra package. In Thyra the user will find a broad set of interfaces and an incredibly useful capability to provide user-defined multivector/vector operations (see the RTOp package). Furthermore, Epetra multivectors/vectors are fully compatible with Thyra.

Example Codes

The following example code generates a simple tridiagonal matrix of dimension "n" where "n" is passed in as an argument. The matrix is then used by the power method to compute the largest eigenvalue and corresponding eigenvector.

The point of this example is to illustrate the flow of calls when using Epetra. This example program can be found in the file Trilinos/packages/epetra/examples/petra_power_method/cxx_main.cpp.

//@HEADER
// ************************************************************************
// 
//               Epetra: Linear Algebra Services Package 
//                 Copyright (2001) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ************************************************************************
//@HEADER

#include <cstdio>
#include <cstdlib>
#include <cassert>
#include <string>
#include <cmath>
#include <vector>
#include "Epetra_Map.h"
#include "Epetra_Time.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#ifdef EPETRA_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#endif
#ifndef __cplusplus
#define __cplusplus
#endif
#include "Epetra_Comm.h"
#include "Epetra_SerialComm.h"
#include "Epetra_Version.h"

// prototype
int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
		 bool verbose);

 
int main(int argc, char *argv[])
{
  int ierr = 0, i;

#ifdef EPETRA_MPI

  // Initialize MPI

  MPI_Init(&argc,&argv);

  Epetra_MpiComm Comm( MPI_COMM_WORLD );

#else

  Epetra_SerialComm Comm;

#endif

  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();
  bool verbose = (MyPID==0);

  if (verbose)
    cout << Epetra_Version() << endl << endl;

  cout << Comm << endl;

  // Get the number of local equations from the command line
  if (argc!=2)
   {
     if (verbose) 
       cout << "Usage: " << argv[0] << " number_of_equations" << endl;
    std::exit(1);
   }
  int NumGlobalElements = std::atoi(argv[1]);

  if (NumGlobalElements < NumProc)
      {
     if (verbose)
       cout << "numGlobalBlocks = " << NumGlobalElements 
	    << " cannot be < number of processors = " << NumProc << endl;
     std::exit(1);
      }

  // Construct a Map that puts approximately the same number of 
  // equations on each processor.

  Epetra_Map Map(NumGlobalElements, 0, Comm);
  
  // Get update list and number of local equations from newly created Map.

  int NumMyElements = Map.NumMyElements();

  std::vector<int> MyGlobalElements(NumMyElements);
    Map.MyGlobalElements(&MyGlobalElements[0]);

  // Create an integer vector NumNz that is used to build the Petra Matrix.
  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
  // on this processor

    std::vector<int> NumNz(NumMyElements);

  // We are building a tridiagonal matrix where each row has (-1 2 -1)
  // So we need 2 off-diagonal terms (except for the first and last equation)

  for (i=0; i<NumMyElements; i++)
    if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
      NumNz[i] = 2;
    else
      NumNz[i] = 3;

  // Create a Epetra_Matrix

  Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
  
  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1


  std::vector<double> Values(2);
  Values[0] = -1.0; Values[1] = -1.0;
  std::vector<int> Indices(2);
  double two = 2.0;
  int NumEntries;
  
  for (i=0; i<NumMyElements; i++)
    {
    if (MyGlobalElements[i]==0)
      {
	Indices[0] = 1;
	NumEntries = 1;
      }
    else if (MyGlobalElements[i] == NumGlobalElements-1)
      {
	Indices[0] = NumGlobalElements-2;
	NumEntries = 1;
      }
    else
      {
	Indices[0] = MyGlobalElements[i]-1;
	Indices[1] = MyGlobalElements[i]+1;
	NumEntries = 2;
      }
     ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
     assert(ierr==0);
     // Put in the diagonal entry
     ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
     assert(ierr==0);
    }
   
  // Finish up
  ierr = A.FillComplete();
  assert(ierr==0);

  // Create vectors for Power method


  // variable needed for iteration
  double lambda = 0.0;
  int niters = NumGlobalElements*10;
  double tolerance = 1.0e-2;

  // Iterate
  Epetra_Flops counter;
  A.SetFlopCounter(counter);
  Epetra_Time timer(Comm);
  ierr += power_method(A, lambda, niters, tolerance, verbose);
  double elapsed_time = timer.ElapsedTime();
  double total_flops =counter.Flops();
  double MFLOPs = total_flops/elapsed_time/1000000.0;

  if (verbose) 
    cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;

  // Increase diagonal dominance
  if (verbose) 
    cout << "\nIncreasing magnitude of first diagonal term, solving again\n\n"
		    << endl;

  if (A.MyGlobalRow(0)) {
    int numvals = A.NumGlobalEntries(0);
    std::vector<double> Rowvals(numvals);
    std::vector<int> Rowinds(numvals);
    A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]
    for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;

    A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
  }
 
  // Iterate (again)
  lambda = 0.0;
  timer.ResetStartTime();
  counter.ResetFlops();
  ierr += power_method(A, lambda, niters, tolerance, verbose);
  elapsed_time = timer.ElapsedTime();
  total_flops = counter.Flops();
  MFLOPs = total_flops/elapsed_time/1000000.0;

  if (verbose) 
    cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;


  // Release all objects
#ifdef EPETRA_MPI
  MPI_Finalize() ;
#endif

/* end main
*/
return ierr ;
}

int power_method(Epetra_CrsMatrix& A, double &lambda, int niters, 
		 double tolerance, bool verbose) {  

  Epetra_Vector q(A.RowMap());
  Epetra_Vector z(A.RowMap());
  Epetra_Vector resid(A.RowMap());

  Epetra_Flops * counter = A.GetFlopCounter();
  if (counter!=0) {
    q.SetFlopCounter(A);
    z.SetFlopCounter(A);
    resid.SetFlopCounter(A);
  }

  // Fill z with random Numbers
  z.Random();

  // variable needed for iteration
  double normz, residual;

  int ierr = 1;

  for (int iter = 0; iter < niters; iter++)
    {
      z.Norm2(&normz); // Compute 2-norm of z
      q.Scale(1.0/normz, z);
      A.Multiply(false, q, z); // Compute z = A*q
      q.Dot(z, &lambda); // Approximate maximum eigenvalue
      if (iter%100==0 || iter+1==niters)
	{
	  resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
	  resid.Norm2(&residual);
	  if (verbose) cout << "Iter = " << iter << "  Lambda = " << lambda 
			    << "  Residual of A*q - lambda*q = " 
			    << residual << endl;
	} 
      if (residual < tolerance) {
	ierr = 0;
	break;
      }
    }
  return(ierr);
}

Epetra adapters for Thyra and PyTrilinos

The Epetra pacakge can also optionally support adapters to various interfaces. Support for these adapters can be turned on at configure time.

• Epetra to Thyra Operator/Vector Adapters: These adapters take Epetra objects and turn them into Thyra objects.

• Epetra to PyTrinos Adapters: These are Python wrappers that allow Epetra to be used from Python.

Browse all of Epetra as a single doxygen collection

You can browse all of Epetra as a single doxygen collection. Warning, this is not the recommended way to learn about Epetra software. However, this is a good way to browse the directory structure of epetra, to locate files, etc.

---