Vector and Matrix Operations in PETSc (via petsc4py)

• Post by: Philip Luke K
• September 30, 2021
• 1 Comment

For the past couple of days, I have been trying to form a basic understanding of PETSc and specifically, its implementation in Python, which is achieved through a library called petsc4py. A normal attempt on anything related to Python will bring you so many detailed web posts and documentation. But that was not the case for petsc4py. There are not many resources or detailed documentation available. Whatever is presently available has its own limitations and takes a coding pro to extract the nuances. So this post is an attempt to document my understanding and bridge the gap for a coding-wannabe like me.

Through this post, you will learn how to do basic vector and matrix operations through petsc4py and how to solve a small linear problem using petsc4py. In addition, FEniCS PETSc formats and conversions are also briefly mentioned.

At this point in time, I am an absolute beginner in Scientific Computing in general; PETSc and FEniCS in particular. If you find some mistakes in my understanding or know a better way out for doing certain things; kindly let me know in the comments.

What is PETSc

The Portable, Extensible Toolkit for Scientific Computation (PETSc, pronounced PET-see; the S is silent), is a suite of data structures and routines developed by Argonne National Laboratory for the scalable (parallel) solution of scientific applications modeled by partial differential equations. It employs the Message Passing Interface (MPI) standard for all message-passing communication. PETSc is the world’s most widely used parallel numerical software library for partial differential equations and sparse matrix computations.

Wikipedia

What is petsc4py

PETSc is primarily written on C++. petsc4py is simply a python binding for PETSc to make it accessible for python users.

Except that it does not! 

Unfortunately, petsc4py doesn’t have a well-established tutorial or documentation page like you would expect from any standard Python Library, despite its terrific applicability. An under-development documentation exists with only Overview and Installation Procedures listed. An unexplained list of syntax/API for petsc4py may be obtained here. The homepage for every official petsc4py development tools (documentation, git repo, API reference, source code etc… ) may be found here

Installing petsc4py

My particular application of PETSc is in FEniCS and FEniCS require Docker images to build. The Docker image, has already been made for our research group, which already contains the dependencies of petsc4py. Hence I need not install it again.

But if you are an absolute new user, you need to install petsc4py library in your system. The installation procedure is pretty straightforward as explained here.

Importing petsc4py

Once petsc4py has been successfully installed, the next task is to import it for use in a code. Go to your favorite text editor or IDE (jupyter, spyder) and type/execute the following :

import petsc4py
petsc4py.init()
from petsc4py import PETSc

Vector class in petsc4py

The vector class of petsc4py is concerned with the creation and operations of Vectors/1-D arrays. A list of Vector class functions, without any explanation whatsoever, may be found here. In this post, I will give out some code snippets for basic vector operations

Creating a vector – By default, the values of all elements of the vector of size n just after creation will be 0. As usual, indexing starts at 0.

n = 10 # Size of vector
x = PETSc.Vec().createSeq(n) # Faster way to create a sequential vector.

Assigning values to a vector – The major syntax is setValues(self, indices, values, addv=None) , where the parameters have the usual meaning. This is illustrated below, where the indices which ranges between 0-9 (which is n-1, n = 10) are given the values 0-9

x.setValues(range(n), range(n))

Accessing the vector – To display the entire vector, use getArray(self, readonly=False). To display only certain values, use getValue(self, index) or getValues(self, indices, values=None)

print(x.getArray())
print(x.getValues(3))
print(x.getValues([1, 2]))

Adding a scalar to the vector – shift(self, alpha) where alpha is the value to be added

x.setValues(range(n), range(n))
x.shift(1)
print(x.getArray())
x.shift(-1)
print(x.getArray())

Adding 2 vectors – Straightforward operation. Vec x + Vec y. Valid only if vector sizes are equal

x.setValues(range(n), range(n))
y.setValues(range(n), range(n))
z=x+y
print(z.getArray())

Sum/Minimum/Maximum – sum(self)/min(self)/max(self). Output of min/max will be a tuple with (position, value)

print(x.sum())
print(x.min())
print(x.max())

Dot product – dot(self, Vec vec)

print(x.dot(x)) # dot product with self
print(x.dot(y)) # dot product with another vector

L2 norm and Infinity norm – norm(self, norm_type=None)

print ('2-norm =', x.norm())
print ('Infinity-norm =', x.norm(PETSc.NormType.NORM_INFINITY))

Matrix class in petsc4py

The Matrix class of petsc4py is concerned with creation and operations on Matrices. A list of Matrix class functions, without any explanation whatsoever, may be found here. The creation and manipulation of matrices are slightly complex compared to vector operations.

Creating and Assembling a Matrix – PETSc provides the option to create many type of matrices such as a dense matrix (Dense), sparse matrix (AIJ. Default PETSc matrix representation)m sequential sparse matrix (SeqAIJ) etc.. A detailed information on the theory behind may be found at the PETSc documentation. Here I will limit to sparse matrices. Also, matrices in PETSc in general should be “assembled” before they can be used. An example syntax of making an m x n sparse matrix is given below. Note that, upon creation, all the elements are by default assigned the value 0

m, n = 4, 4 # size of the matrix
A = PETSc.Mat().createAIJ([m, n]) # AIJ represents sparse matrix
A.setUp()
A.assemble()

Accessing the Matrix – Similar to vectors, the getValues command is used to access and print the matrix

print(A.getValues(range(m), range(n)))
print(A.getValues(range(2), range(1)))

Assigning values to the matrix – setValue(self, row, col, value, addv=None) will assign the value to the mentioned [row][col] index of the matrix and setValues(self, rows, cols, values, addv=None) will assign values to the positions obtained from the tensor product of rows and cols.

A.setValue(1, 1, -1)
A.setValue(0, 0, -2)
A.setValue(2, 2, -5)
A.setValue(3, 3, 6)                        
A.setValues([0, 1], [2, 3], [1, 1, 1, 1])
A.assemble()
print(A.getValues(range(m), range(n)))

Size and Transpose of a vector – The command getSize(self) returns the shape of the matrix as a tuple and transpose(self, Mat out=None) returns the transpose of a matrix. Note that doing the transpose operation (even with as assignment operator) transposes the inherent size of the matrix. So it is advisable to copy the content of the matrix using copy(self, Mat result=None, structure=None) command before transposing.

print(A.getSize())
B = A.copy()
B.transpose()
print(A.getSize(), B.getSize())
print(B.getValues(range(4), range(4)))

Matrix Multiplication – If two matrices have the compatible size, they can be multiplied using the matMult(self, Mat mat, Mat result=None, fill=None) command.

C = A.matMult(B)
print(C.getValues(range(m), range(n)))

Multiplying matrix and vector – mult(self, Vec x, Vec y) multiply matrix with vector x where y is the resulting vector. The y vector also has to be defined

x = PETSc.Vec().createSeq(4) # making the x vector
x.set(1) # assigning value 1 to all the elements
y = PETSc.Vec().createSeq(4) # Put answer here.
A.mult(e, y) # A*e = y
print(y.getArray())

Solving a linear system in petsc4py

The KSP class denoting the krylov subspace solver is the simplest linear solver in PETSc. Here I try to solve a very simple linear system of the form $$Ax = b$$ using petsc4py with the help of the Vector, Matrix and KSP classes available in petsc4py. The code is self-sufficient.

print("Matrix A: ")
print(A.getValues(range(m), range(n))) # printing the matrix A defined above

b = PETSc.Vec().createSeq(4) # creating a vector
b.setValues(range(4), [10, 5, 3, 6]) # assigning values to the vector

print('\\n Vector b: ')
print(b.getArray()) # printing the vector 

x = PETSc.Vec().createSeq(4) # create the solution vector x

ksp = PETSc.KSP().create() # creating a KSP object named ksp
ksp.setOperators(A)

# Allow for solver choice to be set from command line with -ksp_type <solver>.
ksp.setFromOptions()
print ('\\n Solving with:', ksp.getType()) # prints the type of solver

# Solve!
ksp.solve(b, x) 

print('\\n Solution vector x: ')
print(x.getArray())

ksp.setOperators(A), ksp.setFromOptions(), ksp.getType() and ksp.solve(b, x) are functions in the PETSc KSP class. A detailed information may be found in the PETSc KSP documentation

PETSc and petsc4py in FEniCS

FEniCS/DOLFIN has several inherent PETSc based classes such as dolfin.cpp.la.PETScMatrix, dolfin.cpp.la.PETScVector, dolfin.cpp.la.PETScKrylovSolver etc… formed from wrapping PETSc. A rough documentation may be found here. But sometimes it may be needed to used petsc4py to have all the PETSc functionalities properly utilized in FEniCS. Here, I will show some basic code snippets to convert matrices from DOLFIN PETSc, and Numpy to petsc4py.

Consider an eigen value problem with the following variational form

k_form = inner(sigma(u),eps(v))*dx
m_form = rho*dot(u,v)*dx

If the system is large, we need to adopt to PETSc for further calculations and manipulations before solving. To convert the above system first into PETSc DOLFIN, we use the following code format

K = PETScMatrix()
assemble(k_form, tensor=K)
M = PETScMatrix()
assemble(m_form, tensor=M)
solve....
. 
.
.
print(type(K), type(M))

To convert a matrix of the above format into a petsc4py object, we can use the command as_backend_type(mat).mat()

Kp = as_backend_type(K).mat()
Mp = as_backend_type(M).mat()
print(type(Kp), type(Mp))

Suppose I have a numpy matrix X. To convert it into petsc4py matrix object, I can use the following code

Xp = PETSc.Mat().createAIJ(X.T.shape)
Xp.setUp()
Xp.setValues(range(0, X.shape[1]), range(0, X.shape[0]), X.T) # copy X matrix into petsc matrix
Xp.assemble()
Xp.transpose()
print(Xp.getValues(range(0, X.shape[0]), range(0, X.shape[10])))

Final Note

Unlike most other popular Python libraries, petsc4py doesn’t spoonfeed you. There aren’t many shortcuts out of this rabbit hole. Use the available forums and repositories to find your way around!

References

PETSc Homepage
PETSc Documentation
petsc4py Homepage
petsc4py under-development Documentation
petsc4py API Reference
DOLFIN PETSc wrappers