# A simple eigenvalue solver in DOLFIN-X

Hi,
I am trying to write a simple eigenvalue solver code in dolfin-x, but I have a problem in how to define a matrix as a PETScMatrix() and assembling the form into it. There is a proposed code in dolfin which I am trying to write it in dolfin-x:

``````from dolfin import *
mesh = …….
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A)

# Compute all eigenvalues
eigensolver.solve()

# Extract largest (first) eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(0)

print "Largest eigenvalue: ", r
``````

in the following, There is an ambiguity how to `# create eigensolver` and `# compute the eigenvalues` in DOLFIN-X? Thanks in advance.

1 Like

Hi Karo,

You can assemble the PETSc matrix by using the `dolfinx.fem.assemble_matrix` function in dolfinx. An example is given in this demo (https://github.com/FEniCS/dolfinx/blob/master/python/demo/stokes-taylor-hood/demo_stokes-taylor-hood.py). Refer to line 435-436.

Once you get the PETSc matrix. You can use slepc4py to compute the eigenvalues. You can refer to this slepc4py demo (https://gitlab.com/slepc/slepc4py/-/blob/master/demo/ex1.py) and adapt it to your need.

3 Likes

Disclaimer: This is something I tested quite sometime back with dolfinx-`2019.2.9.99` so no guarantees that will work as is, but should be a good starting point in case if it doesn’t too

``````import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh
from dolfinx.fem import assemble_matrix, Form
from dolfinx.plotting import plot
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
# form = Form(a)

# Assemble stiffness tensor
A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)

eigensolver.solve()
vr, vi = A.getVecs()

lmbda = eigensolver.getEigenpair(0, vr, vi)
u = Function(V)
u.vector.setArray(vr.array)

# Plot eigenfunction (0-th)
plt.figure(figsize=(8,8))
eig1 = plot(u, cmap=plt.cm.jet)
plot(mesh)
plt.colorbar(eig1)
``````

translated from

``````import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)
eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
u = Function(V)
u.vector()[:] = rx

# Plot eigenfunction
plt.figure(figsize=(8,8))
eig1 = plot(u, cmap=plt.cm.jet)
plot(mesh)
plt.colorbar(eig1)
``````

4 Likes

I read mentioned demos, thanks for your kind guidance.

Dear @bhaveshshrimali
Thanks for your help, it works for me properly too, and it would be really helpful.

As I found, `vr, vi = A.getVecs()` extracts the real and complex part of the corresponding eigenvector. Do you have an idea how is it possible to extract the real and complex part of the eigenvalue in dolfin-x?
In dolfin, `r, c, rx, cx = eigensolver.get_eigenpair(0)`, extracts both eigenvalue and eigenvector components.

`lmbda = eigensolver.getEigenpair(0, vr, vi)` and `lmbda.real`, `lmbda.imag` gives the real and imaginary part of eigen values. Check the example above for using slepc

Dear Karo,

`vr, vi = A.getVecs()` creates a placeholder vector to store the eigenvector of the matrix A. Following on the comment by bhaveshshrimali, once you called `lmbda = eigensolver.getEigenpair(0, vr, vi)`, it will gives you the first eigenvalue and the first eigenvector. If you would like the rest eigenvalue, you can use loop just like the slepc4py example above.

Here is the documentation on slepc4py (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS-class.html#getEigenpair) for the `getEigenpair` method.

It is exactly my new problem about extracting more eigenvalues.
In this code, dolfin returns selected numbers of eigenvalues:

``````import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

A = PETScMatrix()
assemble(a, tensor=A)

eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()

for i in range (10):
r, c, rx, cx = eigensolver.get_eigenpair(i)
print(r)
``````

The result is:

``````7.822379233170962
7.5635582373483015
7.5635582373482855
7.306318516404323
7.15614193137886
7.156030031719445
6.901377998228639
6.901377998228612
6.6339159619522885
6.6339159619522565
``````

But dolfin-x just returns the first eigenvalue:

``````import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)

eigensolver.solve()
self = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % self )
if self > 0:
for i in range (self):
l = eigensolver.getEigenpair(i ,vr, vi)
print(l.real)
``````

and The result is:

``````Number of converged eigenpairs 1
7.822379233170955
``````

and it doesn’t return more eigenvalues. Please let me know if there is a solution to solve this problem. Thanks in advance.

You can set the number of eigenvalues that you want using the `setDimensions()` method. Here is your code that I’ve modified:

``````import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.setDimensions(nev=10)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
for i in range (evs):
l = eigensolver.getEigenpair(i ,vr, vi)
print(l.real)
``````

You can refer to the documentation for more info.

2 Likes

Thanks a lot dear @adeebkor, @bhaveshshrimali,
The last topic about eigenvalue solver which I am interested in is :

`esolver.parameters["spectrum"] = "smallest real"`

this command helps to arrange eigenvalues from the smallest real part.
But I couldn’t find the command on the subject for dolfin-x, so please indicate if you have a source for this.

Hi,
A vague guess would be to just extract the eigen values using `esolver.getEigenvalue(k)` which should extract the `k` eigen values and then sort them yourself. However, looking up the corresponding function in old dolfin documentation for what it is doing under the hood should give you an idea. See this for the corresponding method, which should roughly translate to

``````eigensolver.setWhichEigenpairs(eigensolver.Which.LARGEST_REAL)
eigensolver.solve()
``````

Adeeb can maybe provide a more detailed reply

1 Like

Most of the commands that I’ve wrote are slepc4py command. You need to look in the slepc4py documentation which is a separate package from dolfin-x.

Nevertheless, the command that you are looking for is `eigensolver.setWhichEigenpairs()`. This is a modified code that may return the smallest real eigenvalue that you need.

``````import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.setDimensions(nev=10)
eigensolver.setWhichEigenpairs(4)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
for i in range (evs):
l = eigensolver.getEigenpair(i ,vr, vi)
print(l.real)
``````

This page list the desired piece of spectrum (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS.Which-class.html) and this page (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS-class.html) is the documentation for the `SLEPc.EPS` class.