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.

Communicators

If PyTrilinos was compiled with MPI support, then MPI_Init() will be called 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. For the following constructors, 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.

    • 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)

    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 SizeArray 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 implementations:

    • 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, 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:

    • Map(numGE, iBase, comm)
    • Map(numGE, numME, iBase, comm)
    • Map(numGE, myGEs, iBase, comm)
    • Map(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 are currently being 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 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 new 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.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

  • Note: Simple python wrappers are provided for the Epetra.FEVector class, but it has not yet been integrated with numpy arrays.

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:

  • Epetra.IntSerialDenseMatrix
  • Epetra.IntSerialDenseVector
  • Epetra.SerialDenseMatrix
  • Epetra.SerialDenseVector

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

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

The Epetra.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 Epetra.CrsGraph class has two new constructors:

  • Epetra.CrsGraph(Epetra.DataAccess CV, Epetra.BlockMap rowMap, PyObject numIndicesList, bool staticProfile=False)
  • Epetra.CrsGraph(Epetra.DataAccess CV, Epetra.BlockMap rowMap, Epetra.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 Epetra.CrsGraph methods have simplified argument lists:

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

The following Epetra.CrsGraph methods do not have python wrappers:

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

Operators

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++ routines can call back to your python method successfully.

  • Epetra.Operator
  • Epetra.InvOperator
  • Epetra.RowMatrix

The following classes provide the same functionality, but will be deprecated in the future:

  • Epetra.PyOperator
  • Epetra.PyRowMatrix

The Epetra.CrsMatrix class has two new constructors:

  • Epetra.CrsMatrix(Epetra.DataAccess CV, Epetra.Map rowMap, PyObject numIndicesList, bool staticProfile=False)
  • Epetra.CrsMatrix(Epetra.DataAccess CV, Epetra.Map rowMap, Epetra.Map 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 Epetra.CrsMatrix methods have simplified argument lists:

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

The following Epetra.CrsMatrix methods do not have python wrappers:

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

The Epetra.CrsMatrix class also has indexing enabled. Thus, if A is an Epetra.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.]