Trilinos/Epetra: Linear Algebra Services Package.

TrilinosRelease4.0Branch

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, MPI 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.

  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.

  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. Both the Epetra_CrsMatrix and Epetra_VbrMatrix classes implement the Epetra_RowMatrix interface and therefore objects from either of these classes can be used with Epetra_LinearProblem and with AztecOO. (See the AztecOO home page.) Furthermore, any class that implements Epetra_RowMatrix can be used with Epetra_LinearProblem and AztecOO.

  7. Sparse row matrix class: Epetra_CrsMatrix - Real double precision sparse matrix class. Supports construction and use of row-wise sparse matrices.

  8. Sparse block row matrix class: Epetra_VbrMatrix - Real double precision block sparse matrix class. Supports construction and use of row-wise block sparse matrices.

  9. 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.

  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.

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 a 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.

Epetra C and Fortran Wrappers

Epetra supports a full set of C and Fortran wrappers. Currently these wrappers are not up to date.

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/examples/Epetra_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 <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <assert.h>
#include <string.h>
#include <math.h>
#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);
  int size, rank; // Number of MPI processes, My process ID

  MPI_Comm_size(MPI_COMM_WORLD, &size);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );

#else

  int size = 1; // Serial case (not using MPI)
  int rank = 0;
  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;
    exit(1);
   }
  int NumGlobalElements = atoi(argv[1]);

  if (NumGlobalElements < NumProc)
      {
     if (verbose)
       cout << "numGlobalBlocks = " << NumGlobalElements 
	    << " cannot be < number of processors = " << NumProc << endl;
     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();

  int * MyGlobalElements = new int[NumMyElements];
    Map.MyGlobalElements(MyGlobalElements);

  // 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

  int * NumNz = new int[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);
  
  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1


  double *Values = new double[2];
  Values[0] = -1.0; Values[1] = -1.0;
  int *Indices = new int[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;
      }
     assert(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
     // Put in the diagonal entry
     assert(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i)==0);
    }
  
  // Finish up
  assert(A.TransformToLocal()==0);

  // Create vectors for Power method


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

  // 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);
    double * Rowvals = new double [numvals];
    int    * Rowinds = new int    [numvals];
    A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
    for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;

    A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
    delete [] Rowvals;
    delete [] Rowinds;
  }
 
  // 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
  delete [] NumNz;
  delete [] Values;
  delete [] Indices;
  delete [] MyGlobalElements;
			
#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);
}

---