PyTrilinos Logo PyTrilinos Tutorial

PyTrilinos Tutorial

PyTrilinos is a collection of python modules that allows a python programmer to access selected Trilinos packages, either dynamically or within a script. For example, once PyTrilinos is installed, a user running python interactively who wanted to gain access to Epetra classes would type

>>> from PyTrilinos import Epetra

To import a PyTrilinos module, use the from PyTrilinos import <module> syntax, where <module> generally corresponds to the Trilinos namespace you want to access. (Note the Epetra C++ package does not use namespaces, and a design decision was made to strip the Epetra_ prefix from Epetra classes and put them in a python namespace Epetra.

This tutorial will touch on the following packages:

No package is completely wrapped. See the python dir() or help() function or the PyTrilinos web documentation for a list of those classes within each package that are wrapped.

PyTrilinos supports MPI, allowing for parallel python scripts. The Epetra.PyComm() function will return an Epetra_MpiComm communicator if Trilinos was built with MPI support, and an Epetra_SerialComm communicator otherwise. In addition, if Trilinos was configured with MPI support, then the Epetra module will internally call MPI_Init() when imported and register MPI_Finalize() with the atexit python module. Thus, scripts can be written for either environment transparently.

Epetra

To see what is available in the Epetra namespace imported above, do a

>>> dir(Epetra)

This will print a list of all the attributes (classes, functions and objects) in the namespace. Currently, you will get a list of well over 600 strings. Some of these will be familiar Epetra names (with the Epetra_ prefix removed), such as 'SerialComm'. A SerialComm communicator can be created with

>>> comm = Epetra.SerialComm()

This object has (almost) all the methods of the C++ Epetra_SerialComm class and is recognized internally by python as an Epetra_SerialComm object. The dir() function works on objects, too:

>>> dir(comm)

gives a (much shorter) list of the attributes for comm, including all of the methods that can be called on the object. For example,

>>> comm.NumProc()
1
>>> comm.MyPID()
0
>>> comm.Label()
'Epetra::Comm'

For more detailed help, use

>>> help(comm)

or access the documentation pages of the appropriate Trilinos package. PyTrilinos provides a wrapper for the Epetra_MpiComm class and introduces the Epetra.PyComm() function, that returns an Epetra_MpiComm object if Trilinos was configured with MPI support and an Epetra_SerialComm otherwise.

We can use comm to create an Epetra_Map (which provides support for distributing global indexes across processors) by typing

>>> map = Epetra.Map(9,0,comm)

Typically, PyTrilinos classes support a python __str__() method, which internally calls the object's Print() method (or print method, as appropriate) and is in turn used by the python print command:

>>> print map

Number of Global Elements  = 9
Number of Global Points = 9
Maximum of all GIDs        = 8
Minimum of all GIDs        = 0
Index Base                 = 0
Constant Element Size      = 1

Number of Local Elements   = 9
Number of Local Points  = 9
Maximum of my GIDs         = 8
Minimum of my GIDs         = 0

         MyPID           Local Index        Global Index
             0                 0                 0
             0                 1                 1
             0                 2                 2
             0                 3                 3
             0                 4                 4
             0                 5                 5
             0                 6                 6
             0                 7                 7
             0                 8                 8

The Print() (or print) methods typically work by providing a python file object, default sys.stdout, where the C++ classes expect a stream.

We can now use map to create an Epetra_Vector:

>>> vect = Epetra.Vector(map)

The Epetra.Vector class actually inherits from both the Epetra_Vector C++ class and the UserArray python class (which exists in both the Numeric and numpy modules -- release 5.0 and 6.0 are Numeric-compatible, while the development branch and releases 7.0 and 8.0 are numpy-compatible). The constructors ensure that both base classes point to the same data buffer. Thus an Epetra.Vector has all the methods and capabilities of an Epetra_Vector, plus the capabilities of a python array:

>>> vect[4] = 3.14
>>> print vect
[ 0.    0.    0.    0.    3.14  0.    0.    0.    0.  ]
>>> vect.shape = (3,3)
>>> print vect
[[ 0.    0.    0.  ]
 [ 0.    3.14  0.  ]
 [ 0.    0.    0.  ]]
>>> vect[2,2] = 2.718
>>> print vect
[[ 0.     0.     0.   ]
 [ 0.     3.14   0.   ]
 [ 0.     0.     2.718]]

Epetra also allows the creation of (compressed row) sparse matrices. Let us consider here the definition of a matrix corresponding to the discretization of a 1D Laplace problem on a regular Cartesian grid. We store the matrix as an Epetra.CrsMatrix, whose constructor requires a map, and an estimation of the number of nonzeros per row (in this instance, 3):

>>> A = Epetra.CrsMatrix(Epetra.Copy,map,3)

The argument Epetra.Copy specifies CopyMode; see the Epetra documentation for more details.

We can now create the matrix, by adding one row at a time. For each row, indices contains global column indices, and values the corresponding values.

>>> numLocalRows = map.NumMyElements()
>>> for lid in range(numLocalRows):
...   gid = map.GID(lid)
...   if gid == 0:
...     indices = [gid, gid + 1]
...     values  = [2.0, -1.0   ]
...   elif gid == n - 1:
...     indices = [gid, gid - 1]
...     values  = [2.0, -1.0   ]
...   else:
...     indices = [gid, gid - 1, gid + 1]
...     values  = [2.0, -1.0   , -1.0   ]
...   A.InsertGlobalValues(lid, values, indices)
...

Finally, we transform the matrix representation into one based on local indices. The transformation is required in order to perform efficient parallel matrix-vector products and other matrix operations.

>>> ierr = A.FillComplete()

Once the matrix, the solution vector (x) and the right-hand side vector (b) have been created, it is convenient to store them in a Epetra.LinearProblem object, which can be created either as

>>> problem = Epetra.LinearProblem(A, x, b)

or as follows:

>>> problem = Epetra.LinearProblem()
>>> problem.SetOperator(A)
>>> problem.SetLHS(x)
>>> problem.SetRHS(b)

Methods GetMatrix(), GetLHS() and GetRHS() can be used to extract the linear system matrix, the solution vector, and the right-hand side vector, respectively.

Teuchos

The primary purpose for PyTrilinos.Teuchos is support for the ParameterList class, which is used by several Trilinos packages for setting solution parameters, flags and output behavior. Often, this is handled implicitly, as PyTrilinos has been designed to accept python dictionaries in place of ParameterList objects. However, you can build a ParameterList directly:

>>> from PyTrilinos import Teuchos
>>> pList = Teuchos.ParameterList()
>>> pList.set("maxiters", 100)
>>> pList.set("tol", 1.0e-6)
>>> pList.set("precond", "ILUT")

The python version of ParameterList is augmented to behave somewhat like dictionaries:

>>> for name in pList:
...   print name, ":", pList[name]
...
maxiters : 100
precond : ILUT
tol : 1.0e-6

You can convert a ParameterList to an XMLObject:

>>> writer = Teuchos.XMLParameterListWriter()
>>> xmlObj = writer.toXML(pList)
>>> print xmlObj
<ParameterList>
<Parameter name="maxiters" type="int" value="100"/>
<Parameter name="precond" type="string" value="ILUT"/>
<Parameter name="tol" type="double" value="1e-06"/>
</ParameterList>

>>> open("params.xml","w").write(xmlObj.toStr())

You can also read an XML ParameterList from disk:

>>> source = Teuchos.FileInputSource("params.xml")
>>> xmlObj = source.getObject()
>>> reader = Teuchos.XMLParameterListReader()
>>> pList  = reader.toParameterList(xmlObj)

EpetraExt

Module EpetraExt contains several utilities to read and write Epetra objects, in particular maps, vectors, and matrices.

An Epetra.Map object can be saved on file using the command

>>> from PyTrilinos import EpetraExt
>>> filename = "map.mm"
>>> EpetraExt.BlockMapToMatrixMarketFile(filename, map)

or can be read from a file as

>>> (ierr, map2) = EpetraExt.MatrixMarketFileToBlockMap(filename, comm)

where ierr is the return error code. Analogously, Epetra.MultiVector and Epetra.CrsMatrix objects can be saved in a file as follows:

>>> EpetraExt.MultiVectorToMatrixMarketFile("x.mm", x)
>>> EpetraExt.RowMatrixToMatrixMarketFile("A.mm", A)

then read as

>>> (ierr, x2) = EpetraExt.MatrixMarketFileToMultiVector("x.mm", map)
>>> (ierr, A2) = EpetraExt.MatrixMarketFileToCrsMatrix("A.mm", map)

These functions are a powerful tool to exchange data between codes written in C++ using Trilinos and codes written in python using PyTrilinos.

See the EpetraExt documentation for more details.

TriUtils

This module, imported with the command

>>> from PyTrilinos import TriUtils

allows the creation of several matrices, in a way that mimics the MATLAB gallery function, and it can be useful for examples and testing.

Several matrices can be easily generated using module TriUtils. For details, we refer to the Trilinos tutorial (Chapter 5). Here, we just show how to generate a matrix corresponding to a 3D Laplacian on a structured Cartesian grid. Let nx, ny, nz be the number of nodes along the x-, y- and z-axes, respectively, and comm be the communicator previously created (see the Epetra module). Then, we can simply write:

>>> nx = 100
>>> ny = 100
>>> nz = 100
>>> Gallery = TriUtils.CrsMatrixGallery("laplace_3d", comm)
>>> Gallery.Set("nx", nx)
>>> Gallery.Set("ny", ny)
>>> Gallery.Set("nz", ny)

The linear system matrix, solution and right-hand side are obtained as

>>> A = Gallery.GetMatrix()
>>> x = Gallery.GetStartingSolution()
>>> b = Gallery.GetRHS()

These objects are automatically destroyed when the Gallery object is deleted. See the TriUtils documentation for more details. Note also that there is also a more advanced gallery of test and example problems available in the PyTrilinos.Galeri module.

Amesos

All Amesos objects are constructed from the function class Amesos. The main goal of this class is to allow the user to select any supported and enabled direct solver, simply by changing an input parameter. Let us suppose that Amesos has been configured and compiled with support for SuperLU. To solve a linear system with SuperLU, we first need to create a Solver object,

>>> from PyTrilinos import Amesos, Epetra
>>> factory = Amesos.Factory()
>>> solver = factory.Create("Superlu", problem)

Then, we can perform the symbolic and numeric factorizations using methods

>>> solver.SymbolicFactorization()
>>> solver.NumericFactorization()

The numeric factorization phase will check whether a symbolic factorization exists or not. If not, method SymbolicFactorization() is invoked. Solution is computed using

>>> solver.Solve()

The solution phase will check whether a numeric factorization exists or not. If not, method NumericFactorization() is called. Users must provide the nonzero structure of the matrix for the symbolic phase, and the actual nonzero values for the numeric factorization. Right-hand side and solution vectors must be set before the solution phase.

Note that using the Amesos module the user can use within Python the following packages: KLU, LAPACK, UMFPACK, SuperLU, SuperLU_DIST, TAUCS, PARDISO, DSCPACK, MUMPS, DSCPACK. See the Amesos documentation for more details.

AztecOO

Often, large sparse and distributed linear systems are solved using iterative solvers of Krylov type, like for example conjugate gradient or GMRES. Within PyTrilinos, iterative solvers are accessed via the AztecOO module. As an example, let us consider the set of instructions required to adopt a non-preconditioned CG, with 1550 maximum iterations and a tolerance of 10^(-5) on the relative residual:

>>> from PyTrilinos import AztecOO
>>> solver = AztecOO.AztecOO(A, x, b)
>>> solver.SetAztecOption(AztecOO.AZ_solver, AztecOO.AZ_cg)
>>> solver.SetAztecOption(AztecOO.AZ_precond, AztecOO.AZ_none)
>>> solver.Iterate(1550, 1e-5)

Non-preconditioned methods rarely converge. AztecOO offers one-level overlapping domain decomposition preconditioner (with exact and inexact subdomain solvers) and multilevel preconditioners (see the ML module overview). The first can be specified as follows:

>>> solver.SetAztecOption(AztecOO.AZ_precond, AztecOO.AZ_dom_decomp)
>>> solver.SetAztecOption(AztecOO.AZ_subdomain_solve, AztecOO.AZ_ilu)
>>> solver.SetAztecOption(AztecOO.AZ_overalp, 1)
>>> solver.SetAztecOption(AztecOO.AZ_graph_fill, 1)

For more details on the available parameters, see the AztecOO documentation. Note also that Teuchos.ParameterList is now supported by AztecOO, so you can set the previous options equivalently with

>>> solver.SetParameters({"precond": "dom_decomp",
...                       "subdomain_solve": "ilu",
...                       "overlap": 1,
...                       "graph_fill": 1})

ML

To define a multilevel preconditioner as defined by the ML package, we first have to set up the required parameters in a python dictionary. A list of supported parameter can be found in the ML user's guide. Here, we specify 3 maximum levels, verbose output (10), and symmetric Gauss-Seidel smoother. Aggregates are computed using the Uncoupled scheme.

>>> from PyTrilinos import ML
>>> mlList = {
...  "max levels"        : 3,
...  "output"            : 10,
...  "smoother: type"    : "symmetric Gauss-Seidel",
...  "aggregation: type" : "Uncoupled"
... }

Then, we create the preconditioner and compute it,

>>> prec = ML.MultiLevelPreconditioner(A, False)
>>> prec.SetParameterList(mlList)
>>> prec.ComputePreconditioner()

Finally, we set up the solver, specifying prec as the preconditioner:

>>> solver = AztecOO.AztecOO(A, x, b)
>>> solver.SetPrecOperator(prec)
>>> solver.SetAztecOption(AztecOO.AZ_solver, AztecOO.AZ_cg)
>>> solver.SetAztecOption(AztecOO.AZ_output, 16)
>>> solver.Iterate(1550, 1e-5)

Please check the ML documentation for more details.