PyTrilinos Python Wrappers for Trilinos

Sandia National Laboratories
PyTrilinos Logo PyTrilinos.Epetra

PyTrilinos.Epetra

The Trilinos Epetra package, for historical portability reasons, does not use namespaces. Instead, "Epetra_" is prepended to every class and function name. The python implementation, however, does utilize the Epetra namespace, and the "Epetra_" prefix has been stripped from all of the class and function names. Therefore, Epetra_Object in C++ is implemented as Epetra.Object in python.

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

Fundamental Classes

  • Epetra.Object: this is the base class for the majority of Epetra classes. In C++, it supports a Print() method that takes an output stream as an argument. In the python implementation for this and all derived classes, this method takes an optional file object argument whose default value is standard out.

    A __str__() method has also been added to the Epetra.Object class that returns the results of Print() in a string, so that the print command will work on Epetra objects. The Print() methods are designed to run correctly in parallel, so do not execute print on an Epetra object conditionally on the processor number. For example, do not do

    if comm.MyPID() == 0: print epetra_obj

    or it will hang your code.

  • Epetra.Util: The Sort() method is not supported.

Communicators

If PyTrilinos was compiled with MPI support, then MPI_Init() will be called internally when the Epetra module is imported (the Epetra module will also arrange for MPI_Finalize() to be called upon termination of the python interpreter).

Parallelism in Trilinos (and PyTrilinos) is largely encapsulated in the Epetra communicator classes. In addition to the base class Epetra.Comm and the derived Epetra.SerialComm and Epetra.MpiComm classes, the Epetra.PyComm() function has been added to the python implementation that returns the appropriate communicator type, depending on whether or not PyTrilinos was compiled with MPI support.

The python interfaces for various communicator methods are slightly different than the C++ interfaces. If comm is an Epetra communicator, then the following methods have the given interfaces:

  • comm.Broadcast(numpy.ndarray obj, int root)
  • comm.GatherAll(PyObject obj) -> numpy.ndarray
  • comm.SumAll(PyObject obj) -> numpy.ndarray
  • comm.MaxAll(PyObject obj) -> numpy.ndarray
  • comm.MinAll(PyObject obj) -> numpy.ndarray
  • comm.ScanSum(PyObject obj) -> numpy.ndarray

In the Broadcast method, the numpy.ndarray data from the root processor is broadcast in-place to all of the other processors. In the remaining methods, the PyObject obj input argument must be a python object that can be converted to an integer, long, or double numpy.ndarray, and the return value is a numpy.ndarray of the same type and size. In C++, these routines have integer error return codes. In python, a non-zero return code is converted to an exception.

Maps

Epetra maps describe the distribution of global vector indexes across processors. The python interface for the following classes is slightly modified from the C++ implementations:

  • Epetra.BlockMap is the (concrete) base class for other Epetra maps. It supports a given number of global elements, where each element can support a (possibly variable) number of data points. But it is the elements that are distributed among the processors. The following constructors are supported:

    • BlockMap(numGE, elSize, iBase, comm)
    • BlockMap(numGE, numME, elSize, iBase, comm)
    • BlockMap(numGE, myGEs, elSize, iBase, comm)
    • BlockMap(numGE, myGEs, elSizes, iBase, comm)
    • BlockMap(map)

    where comm is an Epetra communicator; numGE is the integer number of global elements; numME is the integer number of elements on this processor; elSize is the integer element size; iBase is the integer index base (typically 0 or 1); myGEs is a sequence of integers representing the global indexes on this processor; elSizes is a sequence of integers representing the number of data points for each element on this processor; and map is a BlockMap.

    Instead of two C++ RemoteIDList methods, there is only one with the following interface:

    • BlockMap.RemoteIDLiSt(GIDList) -> (PIDList, LIDList, sizeList)

    where GIDList is a sequence of integer global IDs, and PIDList, LIDList and sizeList are numpy.ndarray objects of integers representing the processor IDs, local IDs and element sizes, respectively. Other BlockMap methods are altered to have the following python interfaces:

    • BlockMap.FindLocalElementID(pointID) -> (elementID, elementOffset)
    • BlockMap.MyGlobalElements() -> numpy.ndarray
    • BlockMap.FirstPointInElementList() -> numpy.ndarray
    • BlockMap.ElementSizeList() -> numpy.ndarray
    • BlockMap.PointToElementList() -> numpy.ndarray

  • Epetra.Map is a simpler form of BlockMap, in which the size of each element is restricted to 1. The python implementation supports the following constructors:

    • Map(numGE, iBase, comm)
    • Map(numGE, numME, iBase, comm)
    • Map(numGE, myGEs, iBase, comm)
    • Map(map)

    where comm is an Epetra communicator; numGE is the integer number of global elements; numME is the integer number of local elements on this processor; iBase is the integer index base (typically 0 or 1); myGEs is a sequence of integer global indexes on this processor; and map is an Epetra.Map

  • Epetra.Export The Export class has the following altered python method interfaces:

    • Export.PermuteFromLIDs() -> numpy.ndarray
    • Export.PermuteToLIDs() -> numpy.ndarray
    • Export.RemoteLIDs() -> numpy.ndarray
    • Export.ExportLIDs() -> numpy.ndarray
    • Export.ExportPIDs() -> numpy.ndarray

  • Epetra.Import The Import class has the following altered python method interfaces:

    • Import.PermuteFromLIDs() -> numpy.ndarray
    • Import.PermuteToLIDs() -> numpy.ndarray
    • Import.RemoteLIDs() -> numpy.ndarray
    • Import.ExportLIDs() -> numpy.ndarray
    • Import.ExportPIDs() -> numpy.ndarray

Vectors

One of the most fundamental data structures for scientific computing is the contiguous array of homogeneous scalars. For scientific python, this data structure has historically been provided by the Numeric module. This was followed by numarray, which added capabilities Numeric lacked, but never replicated all of the Numeric module's functionality. These two modules have now been replaced by the numpy module, whose purpose is to provide a single multidimensional array module to python, acceptable by all factions of the scientific python community.

In Epetra, the contiguous array data structure is provided by several vector classes. In general, these classes lack many of the bells and whistles of numpy.ndarray objects, but they provide at least one unique feature: extensive parallel support. Epetra vectors are assumed to be distributed over one or more processors, as specified by a map object.

To give PyTrilinos users maximum flexibility, the python implementations of the Epetra vector classes Epetra.Vector, Epetra.MultiVector, Epetra.FEVector and Epetra.IntVector inherit from the numpy.lib.user_array.container class. These classes are thus instances of both Epetra vectors and numpy.ndarray. The Epetra Print() method provides Epetra output, and the __str__() method provides numpy output. The key to this approach is providing internal constructors that create an Epetra vector and a numpy.ndarray that both point to the same data buffer. These constructors, as well as other methods with alternate python interfaces, are described below

  • Epetra.Vector provides an array of double precision data. PyTrilinos packages will always interpret an Epetra.Vector as a 1D array, but the python implementation allows any shape that is legal for the given length of the vector. Its constructors and methods with new interfaces are

    • Vector(BlockMap map, bool zeroOut=True)

      Create a new Vector with a size and distribution defined by map. Initialize to zero unless zeroOut is False.

    • Vector(BlockMap map, PyObject array)

      Create a Vector that uses a PyObject to initialize the data. The array object can be either a numpy.ndarray or any sequence that can be used to construct a numpy.ndarray, which then provides the data buffer. The length of the array object on each processor must match the number of elements on each processor specified by the map. If not, an exception is raised.

    • Vector(DataAccess CV, MultiVector source, int index)

      Create a Vector by extracting a single vector from a MultiVector, specified by index. CV should be either Epetra.Copy or Epetra.View.

    • Vector(PyObject array)

      Convert array to a numpy.ndarray if necessary and use this to provide the data buffer for a new Vector. The underlying communicator is Epetra.SerialComm, used as the basis for a simple map.

    • Vector(Vector source)

      Copy constructor.

    • ExtractCopy() -> numpy.ndarray

    • ExtractView() -> numpy.ndarray

    • Dot(Vector A) -> double

    • Norm1() -> double

    • Norm2() -> double

    • NormInf() -> double

    • NormWeighted(Vector weights) -> double

    • MinValue() -> double

    • MaxValue() -> double

    • MeanValue() -> double

    • ReplaceGlobalValues(PyObject values, PyObject indices) -> int

    • ReplaceGlobalValues(int blockOffset, PyObject values, PyObject indices) -> int

    • ReplaceMyValues(PyObject values, PyObject indices) -> int

    • ReplaceMyValues(int blockOffset, PyObject values, PyObject indices) -> int

    • SumIntoGlobalValues(PyObject values, PyObject indices) -> int

    • SumIntoGlobalValues(int blockOffset, PyObject values, PyObject indices) -> int

    • SumIntoMyValues(PyObject values, PyObject indices) -> int

    • SumIntoMyValues(int blockOffset, PyObject values, PyObject indices) -> int

  • Epetra.MultiVector is actually the base class for Vector. PyTrilinos packages will always interpret a MultiVector as a 2D array, but the python implementation allows them to have two or more dimensions. Thus the shape attribute can be changed to a tuple of integers of length at least two whose elements' product is the total size of the array.

    • MultiVector(BlockMap map, int n, bool zeroOut=True)

      Create a new MultiVector with n vectors, with length and distribution according to map. Initialize to zero unless zeroOut is False.

    • MultiVector(BlockMap map, PyObject array)

      Create a MultiVector that uses a PyObject to initialize the data. The array object is converted to a numpy.ndarray if necessary, which then provides the data buffer. If the numpy.ndarray is one-dimensional, it is upcast to 2D with a first dimension of one. The product of the second and any subsequent dimensions of the array object on each processor must match the number of elements on each processor specified by the map. If not, an exception is raised.

    • MultiVector(DataAccess CV, MultiVector source, PyObject range)

      Create a MultiVector by extracting a subset of vectors from a MultiVector. The range object should evaluate to a sequence of integers that specify the subset. CV should be either Epetra.Copy or Epetra.View.

    • MultiVector(PyObject array)

      Convert array to a numpy.ndarray if necessary and use this to provide the data buffer for a new MultiVector. The first dimension of the numpy.ndarray specifies the number of vectors. If the numpy.ndarray is one-dimensional, it will be upcast to a 2D array with a first dimension of one. The underlying communicator is Epetra.SerialComm, used as the basis for a simple map.

    • MultiVector(MultiVector source)

      Copy constructor.

    • ExtractCopy() -> numpy.ndarray

    • ExtractView() -> numpy.ndarray

    • Dot(MultiVector a) -> numpy.ndarray

    • Norm1() -> numpy.ndarray

    • Norm2() -> numpy.ndarray

    • NormInf() -> numpy.ndarray

    • NormWeighted(MultiVector weights) -> numpy.ndarray

    • MinValue() -> numpy.ndarray

    • MaxValue() -> numpy.ndarray

    • MeanValue() -> numpy.ndarray

  • Epetra.FEVector derives from Epetra.MultiVector (although it currently supports only a single vector), and provides some methods that are convenient for finite element vector assembly.

  • Epetra.IntVector provides an array of integers. PyTrilinos packages will always interpret an Epetra.IntVector as a 1D array, but the python implementation allows any shape that is legal for the given length of the vector.

    • IntVector(BlockMap map, bool zeroOut=True)

      Create a new IntVector with a size and distribution defined by map. Initialize to zero unless zeroOut is False.

    • IntVector(BlockMap map, PyObject array)

      Create an IntVector that uses a PyObject to initialize the data. The array object is converted to a numpy.ndarray if necessary, which then provides the data buffer. The length of the array object on each processor must match the number of elements on each processor specified by the map. If not, an exception is raised.

    • IntVector(PyObject array)

      Convert array to a numpy.ndarray and use this to provide the data buffer for a new IntVector. The underlying communicator is Epetra.SerialComm, used as the basis for a simple map.

    • IntVector(IntVector source)

      Copy constructor.

    • ExtractCopy() -> numpy.ndarray

    • ExtractView() -> numpy.ndarray

    • Values() -> numpy.ndarray

SerialDense Classes

As with Epetra vector objects, several of the Epetra SerialDense objects multiply inherit from the numpy.lib.user_array.container class. These classes are:

  • IntSerialDenseMatrix
  • IntSerialDenseVector
  • SerialDenseMatrix
  • SerialDenseVector

The following SerialDense class methods have an altered calling signature because the python version is smart enough to determine the dimensions automatically:

  • SerialDenseSolver.IPIV() -> numpy.ndarray
  • SerialDenseSolver.A() -> numpy.ndarray
  • SerialDenseSolver.B() -> numpy.ndarray
  • SerialDenseSolver.X() -> numpy.ndarray
  • SerialDenseSolver.AF() -> numpy.ndarray
  • SerialDenseSolver.FERR() -> numpy.ndarray
  • SerialDenseSolver.BERR() -> numpy.ndarray
  • SerialDenseSolver.R() -> numpy.ndarray
  • SerialDenseSolver.C() -> numpy.ndarray

The SerialDenseSolver.ReciprocalConditionEstimate() method takes no arguments, returns the double precision result and raises an exception if the underlying C++ int result is nonzero.

Graphs

The CrsGraph class has two new constructors:

  • CrsGraph(DataAccess CV, BlockMap rowMap, PyObject numIndicesList, bool staticProfile=False)
  • CrsGraph(DataAccess CV, BlockMap rowMap, BlockMap colMap, PyObject numIndicesList, bool staticProfile=False)

The argument CV is either Epetra.Copy or Epetra.View, rowMap and colMap are maps that describe the domain decomposition of global row indices and column indices respectively, numIndicesList is a python sequence of integers that lists the number of non-zeros for each row, and staticProfile is a boolean that flags whether the number of indices per row is exact or approximate.

The following CrsGraph methods have simplified argument lists:

  • CrsGraph.ExtractGlobalRowCopy(int globalRow) -> numpy.ndarray
  • CrsGraph.ExtractMyRowCopy(int myRow) -> numpy.ndarray
  • CrsGraph.InsertGlobalIndices(int globalRow, PyObject indices) -> int
  • CrsGraph.InsertMyIndices(int myRow, PyObject indices) -> int
  • CrsGraph.RemoveGlobalIndices(int globalRow, PyObject indices) -> int
  • CrsGraph.RemoveMyIndices(int myRow, PyObject indices) -> int

The following CrsGraph methods do not have python wrappers:

  • CrsGraph.ExtractGlobalRowView()
  • CrsGraph.ExtractMyRowView()

Operators

Perhaps the simplest Epetra operators to create and use are the CrsMatrix and VbrMatrix. You simply call their constructors in the normal way, populate them, and then call their FillComplete() methods. These classes have two new constructors:

  • CrsMatrix(DataAccess CV, Map rowMap, PyObject numIndicesList, bool staticProfile=False)
  • CrsMatrix(DataAccess CV, Map rowMap, Map colMap, PyObject numIndicesList, bool staticProfile=False)

and

  • VbrMatrix(DataAccess CV, Map rowMap, PyObject numBlockEntriesPerRow)
  • VbrMatrix(DataAccess CV, Map rowMap, Map colMap, PyObject numBlockEntriesPerRow)

The argument CV is either Epetra.Copy or Epetra.View, rowMap and colMap are maps that describe the domain decomposition of global row indices and column indices respectively, numIndicesList is a python sequence of integers that lists the number of non-zeros for each row, and staticProfile is a boolean that flags whether the number of indices per row is exact or approximate. For the VbrMatrix constructors, numBlockEntriesPerRow can be either an interger, constant for all rows, or a sequence with as many entries as the matrix has rows.

The following CrsMatrix methods have simplified argument lists:

  • CrsMatrix.ExtractGlobalRowCopy(int globalRow) -> (numpy.ndarray,numpy.ndarray)
  • CrsMatrix.ExtractMyRowCopy(int myRow) -> (numpy.ndarray,numpy.ndarray)
  • CrsMatrix.InsertGlobalValues(int globalRow, PyObject indices, PyObject values) -> int
  • CrsMatrix.InsertMyValues(int myRow, PyObject indices, PyObject values) -> int
  • CrsMatrix.RemoveGlobalValues(int globalRow, PyObject indices, PyObject values) -> int
  • CrsMatrix.RemoveMyValues(int myRow, PyObject indices, PyObject values) -> int
  • CrsMatrix.SumIntoGlobalValues(int globalRow, PyObject indices, PyObject values) -> int
  • CrsMatrix.SumIntoMyValues(int myRow, PyObject indices, PyObject values) -> int

The following CrsMatrix methods do not have python wrappers:

  • CrsMatrix.ExtractGlobalRowView()
  • CrsMatrix.ExtractMyRowView()

The CrsMatrix class also has indexing enabled. Thus, if A is an CrsMatrix, matrix elements can be assigned with the syntax

>>> A[i,j] = 2.7182818284590451

Under certain conditions, index notation can be used to retrieve single matrix elements (with a row and column index), or all of the elements in the row (with a single row index). For single matrix elements, a column map must exist, which can be provided in the constructor, or computed automatically when FillComplete() is called. To extract row data with a single index, FillComplete() is required to have been called.

>>> comm = Epetra.PyComm()
>>> map = Epetra.Map(9,0,comm)
>>> A = Epetra.CrsMatrix(Epetra.Copy,map,3)
>>> for gid in map.MyGlobalElements():
...    lid = map.LID(gid)
...    if gid in (0,8): A.InsertGlobalIndices(lid,[1],[gid])
...    else: A.InsertGlobalIndices(lid,[-1,2,-1],[gid-1,gid,gid+1])
...
>>> A.FillComplete()
>>> print A[8,8]
1.0
>>> print A[4]
[ 0.  0.  0. -1.  2. -1.  0.  0.  0.]

A more flexible way of defining Epetra operators is to define your own classes by inheriting from pure virtual Epetra operator classes. The following classes can be sub-classed successfully in python. That is, you can define a python class that inherits from one of these classes, define the appropriate methods, such as Apply(), and the infrastructure is in place such that C++ solver routines (such as the AztecOO.Iterate() method) can call back to your python method successfully.

  • Operator
  • InvOperator
  • RowMatrix
  • BasicRowMatrix

Python classes derived from these callback-supporting base classes ("directors" in the parlance of SWIG, which is the tool that generates the wrapper code) must call their base class __init__() method:

from PyTrilinos import Epetra

class MyOperator(Epetra.Operator):
    def __init__(self):
        Epetra.Operator.__init__(self)
        self.__label = "MyOperator"

At a bare minimum, you must define a Label() method and an Apply() method that support the argument prototypes found in the Epetra documentation:

def Label(self):
    return self.__label
def Apply(self, x, y):
    try:
        y[:,0   ] = x[:,0]
        y[:,1:-1] = 2*x[:,1:-1] - x[:,:-2] - x[:,2:]
        y[:,-1  ] = x[:,-1]
        return 0
    except Exception, e:
        print "A python exception was raised in MyOperator.Apply:"
        print e
        return -1

A few notes about the Apply() method:

  • Arguments x and y will be sent by a C++ solver (such as AztecOO) with C++ type Epetra_MultiVector. Before they are passed along to your python method, they will be converted to the numpy-hybrid type Epetra.MultiVector. Hence we are able, within the method, to access slices of x and y, which is a very efficient way to calculate with buffers of data.
  • Also, because x and y are MultiVector objects, they should be treated as 2-dimensional, with the first dimension representing the individual vectors. As in this example, this usually just means making sure the first index is a colon.
  • The calculations here are done within a try block. This is because if your Apply() accidentally raises an exception, the callback mechanism isn't currently able to pass the error message to you. Hence, we catch all exceptions (base class Exception), print the error message, and return a non-zero error code to tell the calling function that we failed. This technique is recommended for all callback functions of any complexity. It adds almost no overhead and is a great help for debugging.
  • In this particular example, we have taken advantage of python's negative indexing (which indexes from the end of a sequence). For this reason, this operator works on different size vectors. If the input vectors, x and y, have different shapes, then the assignment statements will raise an exception, which we will catch, print, and then inform the calling routine by returning -1.

The operative inheritance chain here is that SrcDistObject and Operator are both base classes for RowMatrix, which is a base class for BasicRowMap. Each of these classes has virtual methods you may want to implement with the given input prototypes and output argument:

  • SrcDistObject
    • Map() -> BlockMap
  • Operator
    • SetUseTranspose(bool useTranspose)
    • UseTranspose() -> bool
    • Apply(MultiVector x, MultiVector y) -> int
    • ApplyInverse(MultiVector x, MultiVector y) -> int
    • HasNormInf() -> bool
    • NormInf() -> double
    • Label() -> string
    • Comm() -> Comm
    • OperatorDomainMap() -> Map
    • OperatorRangeMap() -> Map
  • RowMatrix
    • NumMyRowEntries(int myRow, numpy.ndarray numEntries) -> int
    • MaxNumEntries() -> int
    • ExtractMyRowCopy(int myRow, int length, numpy.ndarray numEntries, numpy.ndarray values, numpy.ndarray indices) -> int
    • ExtractDiagonalCopy(Vector diagonal) -> int
    • Multiply(bool useTranspose, MultiVector x, MultiVector y) -> int
    • Solve((bool upper, bool trans, bool unitDiagonal, MultiVector x, MultiVector y) -> int
    • InvRowSum(Vector x) -> int
    • LeftScale(Vector x) -> int
    • InvColSums(Vector x) -> int
    • RightScale(Vector x) -> int
    • Filled() -> bool
    • NormOne() -> double
    • NumGlobalNonzeros() -> int
    • NumGlobalRows() -> int
    • NumGlobalCols() -> int
    • NumGlobalDiagonals() -> int
    • NumMyNonzeros() -> int
    • NumMyRows() -> int
    • NumMyCols() -> int
    • NumMyDiagonals() -> int
    • LowerTriangular() -> bool
    • UpperTriangular() -> bool
    • RowMatrixRowMap() -> Map
    • RowMatrixColMap() -> Map
    • RowMatrixImporter() -> Import

Of particular importance to the python interface of the RowMatrix and BasicRowMatrix classes are the NumMyRowEntries and ExtractMyRowCopy methods, because these method's C++ argument lists have output arguments: int &NumEntries, double *Values, and int *Indices. To convert these into mutable python objects that can properly act as output arguments, these arguments are converted to numpy.ndarray objects whose data buffers point to the passed-in C++ data. It is natural to access the values and indices as arrays with bracket notation for setting values, but the numEntries argument, which is an array of length 1, must also be accessed this way:

def ExtractMyRowCopy(self, myRow, length, numEntries, values, indices):
    globalRow = self.RowMatrixRowMap().GID(myRow)
    if globalRow == -1:
        return -1
    if globalRow == 0 or globalRow == self.NumGlobalRows()-1:
        if (length < 1):
            return -2
        numEntries[0] = 1
        values[0]     = 1.0
        indices[0]    = myRow
    else:
        if (length < 3):
            return -2
        numEntries[0] = 3
        values[:3]    = [   -1.0,   2.0,    -1.0]
        indices[:3]   = [myRow-1, myRow, myRow+1]
    return 0

This example assumes you have properly implemented the RowMatrixRowMap() and NumGlobalRows() methods.

Generated on: 2010-04-16 18:48 UTC. Generated by Docutils from reStructuredText source.

PyTrilinos Home
About | Downloads | Mail Lists | Documentation | Help

Trilinos Home | Trilinos Developer Home | Trilinos Packages
Sandia | Contact | Disclaimer
Last updated November 29, 2009 - 3:58 am MST