PyTrilinos Logo PyTrilinos Tutorial - Release 7.0

PyTrilinos Tutorial - Release 7.0

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

Release 7.0 of PyTrilinos supports the following packages:

No package is completely wrapped. See the python dir() or help() function or the PyTrilinos 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 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 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 pappropriate) 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 numpy.lib.user_array.container python class. 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 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 keyword Epetra.Copy is the 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 ii in range(numLocalRows):
...   i = map.GID(ii)
...   if i == 0:
...     indices = [i, i + 1]
...     values  = [2.0, -1.0]
...   elif i == n - 1:
...     indices = [i, i - 1]
...     values  = [2.0, -1.0]
...   else:
...     indices = [  i,  i - 1, i + 1]
...     values  = [2.0,   -1.0,  -1.0]
...   A.InsertGlobalValues(i, 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 (lhs) and the right-hand side vector (rhs) have been created, it is convenient to store them in a Epetra.LinearProblem object, which can be created either as

>>> Problem = Epetra.LinearProblem(A, lhs, rhs)

or as follows:

>>> Problem = Epetra.LinearProblem()
>>> Problem.SetOperator(A)
>>> Problem.SetLHS(lhs)
>>> Problem.SetRHS(rhs)

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

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

>>> 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. First, we need to import modules TriUtils and Epetra. 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

>>> Matrix = Gallery.GetMatrix()
>>> LHS = Gallery.GetStartingSolution()
>>> RHS = Gallery.GetRHS()

These objects are automatically destroyed then the Gallery object is deleted.

See the TriUtils documentation for more details.

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 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(Matrix, LHS, RHS)
>>> 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.

ML

To define a multilevel preconditioner as defined by the ML package, we first have to setup the required parameters in a Python's 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(Matrix, False)
>>> Prec.SetParameterList(MLList)
>>> Prec.ComputePreconditioner()

Finally, we set up the solver, and specifies to use Prec as preconditioner:

>>> Solver = AztecOO.AztecOO(Matrix, LHS, RHS)
>>> 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.