Introduction

slepc4py provides a python interface to the SLEPc eigensolver package.  It allows one to easily solve very large and matrix-free linear eigenvalue problems using iterative methods.  The philosophy behind slepc4py is to make powerful iterative methods easy to use even for the beginner.

Installation

slepc4py is developed and tested on Mac OS X and Linux, but might work well with other unix variants.  The prerequisites are:

• PETSc (2.3.3)
• numpy (1.0.0 or higher)
  
• petsc4py (0.7.5 or the latest from subversion)
  
• SLEPc (2.3.3)
• SWIG (1.3.29 or higher)
  

PETSc, and to a lesser extent SLEPc, has a wide variety of ways that it can be configured upon install.  To get slepc4py up and working initially, we recommend to do just the "Quick Install".

Currently, the only way to get slepc4py is to download it from the subversion interface.  First, check to see if the svn binary is available.  If not, and you are on a Redhat/Fedora linux box, you can usually do the following command:

yum install subversion

If you are on Mac OS X (or another linux variant), you might have to build subversion from source.
Once installed, check out the existing sources using this command:

svn co svn://t2.unl.edu/brian/slepc4py

A trac-based source code browser is available here:

http://t2.unl.edu:8092/browser/slepc4py

Once all the source has been downloaded, compile and install using the setup.py script:

python setup.py build install

This will generate the SWIG wrappers, compile them, and install them to your python's modules directory.  The final command - install - may require root permissions on your system.  In the future, when a release of slepc4py is made, we hope to include the generated SWIG wrappers in the tarball.  This will remove the dependency on SWIG for the end-user's system.

Once you believe you have done the install, see if you can import the SLEPc object; if so, you have been successful!

$ python
Python 2.5.1 (r251:54869, Apr 18 2007, 22:08:04) 
[GCC 4.0.1 (Apple Computer, Inc. build 5367)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from slepc4py import SLEPc
>>> 

Example Usage

There are simple, commented examples available in the examples/ directory of the source tree.  These will always be more up-to-date than this documentation.  However, we want to reproduce one or two here to show the programming style used in slepc4py.  This section walks through the example code found in ex1.py.  Please see that source file for the complete program.

The user should be able to access all the SLEPc functionality through the SLEPc module.  To import it, do this:

from slepc4py import SLEPc

(if this fails, your installation has not succeeded).  In addition, we want two more imports; PETSc and numpy:

from petsc4py import PETSc
import numpy

Again, if either of these imports fail, you'll want to try the installation again.  Finally, we want to use PETSc's options database.  This database is a flexible, powerful configuration method that PETSc uses; we recommend using it for your programs also:

opts = PETSc.Options()
n = opts.getInt('n', 30)

There are three general parts to solving an eigenproblem:

1. Set up the eigenproblem.  To do this, you will need to create a linear operator (i.e., PETSc.Mat class).  In this example, we create a matrix A representing the Laplacian operator.
   

   # Create the PETSc matrix, and set its size to `n`
   A = PETSc.Mat(); A.create()
   A.setSizes( (n, n) )
   A.setFromOptions( )
   
   # The ownership range represents the rows of the matrix we own.
   # In the case of sequential programs, this is all of them.
   o = A.getOwnershipRange()
   
   # If we own the first or last row, we must handle them separately
   # and not try to calculate them with the rest of the matrix.
   FirstBlock = 0
   LastBlock  = 0
   
   if (o[0] == 0): FirstBlock = 1
   if (o[1] == n): LastBlock  = -1
   
   # These are the values we will set; 2.0 is the diagonal.
   values = numpy.array( [-1.0,2.0,-1.0] )
   
   # Loop through the rows in our range, and set the values
   # where appropriate.  Look how easy this is!
   for i in range(o[0]+FirstBlock,o[1]+LastBlock):
     A[i,i-1:i+2] = values
   
   # Finally, our special cases for the first and last rows.
   if (LastBlock != 0):
     A[n-1,-2:] = numpy.array([-1.0,2.0])
   if (FirstBlock != 0):
     A[0,:2] = numpy.array([2.0,-1.0])
   
   # Tell PETSc we're done modifying values.
   A.assemble()

   All of the above work is exactly the same as what one might do when they use petsc4py; if any part seems especially puzzling, please refer to petsc4py's documentation.
   
2. Create the SLEPc.EPS (Eigenvalue Problem Solver) object, and set the appropriate options.  Run the solve method.
   

   # Create the python EPS object, then the underlying C one.
   # Both of these must be done for each EPS object used.
   E = SLEPc.EPS()
   E.create()
   
   # Set the operator for the eigenvalue problem to be the `A` object
   # created in the previous section.
   E.setOperators( A )
   
   # Set the problem type; `A` is hermitian, so we use problem type
   # `HEP`: Hermitian Eigenvalue Problem.  Refer to the SLEPc documentation
   # for other problem types.
   E.setProblemType( SLEPc.EPS.ProblemType.HEP )
   
   # Finally, set any additional options passed from the command line; this is
   # where the PETSc options database is extremely useful.
   # Examples of things that can be set from the command line: solution algorithm,
   # convergence tolerances, maximum iterations, number of eigenvalues requested, etc.
   E.setFromOptions()
   
   # Now the EPS object is set up, we call the `solve` routine.
   
   
   E.solve()
   

   This is the extent of the manipulation we must do to the EPS object in order to get the SLEPc machinery going.  Note that it is simple to use, yet there is lots of room for customization.  Also, if you compare the documentation for SLEPc, you will notice that we have tried to keep the Python bindings close to the original C code, except most methods need fewer arguments.  Hence, the bright user can use the SLEPc documentation for a significant portion of slepc4py also.
   

3. Use the output of the EPS object to do something useful!  In this case, we just print out a summary to the screen.
   

   PETSc.Print( "\n******************************\n" )
   PETSc.Print( "*** SLEPc Solution Results ***\n" )
   PETSc.Print( "******************************\n\n" )
   
   its = E.getIterationNumber()
   PETSc.Print( "Number of iterations of the method: %d\n" % its )
   
   eps_type = E.getType()
   PETSc.Print( "Solution method: %s\n\n" % eps_type )
   
   nev, ncv = E.getDimensions()
   PETSc.Print( "Number of requested eigenvalues: %d\n" % nev )
   
   tol, maxit = E.getTolerances()
   PETSc.Print( "Stopping condition: tol=%.4g, maxit=%d\n" % (tol, maxit) )
   
   nconv = E.getConverged()
   PETSc.Print( "Number of converged eigenpairs %d\n\n" % nconv )
   
   if nconv > 0:
       # Create the results vectors
       xr, tmp = A.getVecs()
       xi, tmp = A.getVecs()
   
       PETSc.Print( "           k          ||Ax-kx||/||kx||\n" \
                    "   ----------------- ------------------\n")
       for i in range(nconv):
           k = E.getEigenpair(i, xr, xi)
           error = E.computeRelativeError(i)
           if k.imag != 0.0:
               PETSc.Print( " %9f%+9f j %12g\n" % (k.real, k.imag, error) )
           else:
               PETSc.Print( " %12f       %12g\n" % (k.real, error) )
       PETSc.Print( "\n" )

   Note that we use PETSc.Print instead of python's print method.  This is because, in parallel runs, each process will try to print to stdout using print, while only the root process tries to print when using PETSc.Print.
   

slepc4py Binding Progress

Most of the functionality of the following objects has been wrapped:

• EPS: eigenvalue problem solver.
  
• ST: custom spectral transforms for the EPS object.
  

These all have been done by Brian Bockelman.  However, Lisandro Dalcin has helped me do a significant portion of the underlying framework for slepc4py, and indeed most of the slepc4py framework is taken from petsc4py.  We feel that the less code duplication between the two projects, the better.

The following have not been wrapped:

• IP: custom inner product approximations for the EPS object.  This is only probably used by advanced users; please contact us in case if you need it.
  
• SVD: singular vector decomposition for the matrix operator.  It is highly desired to get this wrapped.  The SLEPc maintainers have a summer student working for them that we hope we can sweet-talk into doing this for us.