# Efficiency of using petsc4py vs. c/c++/fortran

How much slower is petsc4py vs c/c++/fortran?

I realize it will depend significantly on the code being executed, but what about something simple like a matrix-vector product?

This is a widely-held concern in the scientific programming community, and I would consider the performance uncertainty to be one of the major "myths" in computational science.

As @fcruz discusses, petsc4py is a wrapper to the PETSc libraries, not a reimplementation of PETSc in Python. Therefore, you can expect any performance penalties to come from either copying arrays to and from PETSc, or from overhead in your driver code/function calls. petsc4py is very carefully implemented, and as long as you understand the numpy multi-dimensional array interfaces, you can avoid copy overhead. For the majority of use cases I work in, the performance penalty in working in Python is on the order of 10-40%, and I often gain substantially in other ways that more than make up for this performance hit. In fact, several more experienced HPC Python developers I have talked with hold the opinion that this performance difference can usually be reduced even further, and when Python is driving computationally expensive codes, this will certainly be the case.

The petsc4py repository itself features a number of useful examples to illustrate the performance/flexibility tradeoff. Look in the petsc4py source repository for the demo called perftest, which solves a nonlinear system of equations using both a Python driver and a C driver (over a Fortran kernel provided in App.f90 in that directory). The performance overhead here is on the order of 10%.

As a concrete example, I am part of a team of scientists working on PyClaw, a software package that interfaces into PETSc for parallel grid management and legacy Fortran kernels for solving Riemann problems on cell interfaces. We performed a fairly careful study of the performance degradation from switching over from a Fortran driver, and you can see the results on the bottom of page 5 in Table 1 in the conference paper. In our case, we traded a little bit of on-core performance for the capability to easily interface our code to PETSc and Fortran and run efficiently in parallel on tens of thousands of cores.

• I also have a similar concern w.r.t. a small unstructured code. PETSc only provides the data structures and solvers but I still have to read in the mesh (up to 4GB input file), partition, create mappings, loop over the elements, calculate local (element) stiffness matrices etc. before PETSc can assemble and solve. Wouldn't python be slower for this non PETSc related stuff specially I/O, mappings and element level calculations. Because the rest of the code is simple anyway. – stali May 7 '12 at 12:04
• The element-level calculations are usually passed in as a kernel (see the App.f90 source in perftest). There is no performance difference in I/O. Have you looked at FEniCS for a higher-level package? – Aron Ahmadia May 7 '12 at 13:57
• You are right. I do get the idea but in my particular case there are many such kernels (shape functions for different type of elements, element level calculations, mappings etc.) which is about 90% of the code. I did look at Fenics a while back and a lot of details such as dealing with external meshes and imposing BCs etc. were not quite clear at the first glance or seem more complicated (at least to me). Besides I use Fortran which is fairly straightforward to use (given the excellent PETSc documentation) for non CS folks like myself. I actually find it easier than python :) for my work. – stali May 7 '12 at 15:52

Petsc4py is just another way to acces PETSc but from python, or is the same to say that, petsc4py provides the bindings so that, from python, you can access PETSc data structures and routines that are meant to reduce the effort of developing parallel PDE solvers (that scale).

PETSc provides several levels of abstractions to their solvers, and you can even use PETSc to implement your own solver. At the lowest level of software abstraction, PETSc uses BLAS, LAPACK, and MPI, and at best it will be as fast as the implementation of these.

Now, pets4py uses cython to implement the bindings to PETSc. The overhead of using cython is relative to how much computations are going to be done from PETSc. If you use the high level PDE solvers from PETSc, the overheads should be small enough that you don't need to worry about them.

A maybe more important question than the performance comparison of PETSc vs GEMV is if PETSc is the right tool for your job. If you need to implement non-trivial parallel PDE solvers, then most probably, PETSc will really help you. However, if you need to do a bunch of GEMV, you want a BLAS library. Good luck!