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:
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.
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 inplace 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 nonzero return code is converted to an
exception.
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
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 onedimensional, 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 onedimensional, 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
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.
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 nonzeros 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()
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 nonzeros 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],[gid1,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 subclassed 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 callbacksupporting 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 numpyhybrid 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 2dimensional, 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 nonzero
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
 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
passedin 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] = [myRow1, myRow, myRow+1]
return 0
This example assumes you have properly implemented the
RowMatrixRowMap() and NumGlobalRows() methods.
