# How to compute the condition number efficiently?

How to compute the condition number of a problem efficiently? I know one can convert the matrix to dense and use numpy, but I am looking for something more efficient.

from fenics import *
import numpy as np

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
F = inner(v, u) * dx
A = assemble(lhs(F))
np.linalg.cond(A.array())


Consider using the SLEPcEigenSolver:

from dolfin import *
import numpy as np
import time

mesh = UnitSquareMesh(60,60)
V = FunctionSpace(mesh, "CG", 1)
u,v = TrialFunction(V), TestFunction(V)
a = inner(u,v)*dx

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

start = time.time()
eigenSolver = SLEPcEigenSolver(A)
eigenSolver.parameters["spectrum"]="smallest magnitude"
eigenSolver.solve(1)
eigen_min = eigenSolver.get_eigenvalue(0)[0]

eigenSolver.parameters["spectrum"]="largest magnitude"
eigenSolver.solve(1)
eigen_max = eigenSolver.get_eigenvalue(0)[0]
end = time.time()
print("Condition number {0:.2e} ({1:.2e} s)".format(eigen_max/eigen_min,
end-start))

start = time.time()
cond_numpy = np.linalg.cond(A.array())
end = time.time()
print("Numpy condition number {0:.2e} ({1:.2e} s)".format(cond_numpy,
end-start))


Note that you might have to increase the argument in

eigenSolver.solve(1)


when increasing the number of cells in your mesh.

2 Likes

Nice! I was not patient enough to wait for the numpy result in your example . Is there some docs on the eigen solver? E.g. the following questions pop up:

• What does the 1 in eigenSolver.solve(1) mean?
• What is eigenSolver.parameters["spectrum"] about?
• What do the zeros in eigenSolver.get_eigenvalue(0)[0] mean?

Also where can I read about PETScMatrix? e.g. what sparsity structure does it have and how do I linear algebra with it?

See for instance this demo for information about SLEPcEigen Solvers.

PETScMatrix is just a wrapper around PETScâ€™s sparse matrices and its docs can be found here.

1 Like

Also how to add a boundary condition before computing cond?

from fenics import *
import numpy as np

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
F = inner(v, u) * dx
A = assemble(lhs(F))
bc = DirichletBC(V, 1, "on_boundary")
bc.apply(A)
np.linalg.cond(A.array())


1 Like

The SLEPcEigenSolver computes a partial spectrum. The â€śspetrumâ€ť argument lets you choose, which part of the spectrum should be computed. The eigenSolver.solve(1)  means that roughly one eigenvalue should be computed. re, im = solver.get_eigenvalue(i) gives the real and imaginary part of the ith eigenvalue.


def compute_eig0(solver, neig_min=1, neig_max=1000):
neig = neig_min
ret = None
while True:
neig = min(neig, neig_max)
solver.solve(neig)
try:
ret = solver.get_eigenvalue(0)
break
except Exception:
pass
if neig >= neig_max:
raise f"Reached neig_max = {neig_max}"
neig = 2*neig
return ret

def cond(F, bcs):
A = PETScMatrix()
dummy = rhs(F)
assemble_system(lhs(F), dummy, A_tensor=A, bcs=bcs)
eigenSolver = SLEPcEigenSolver(A)
eigenSolver.parameters["spectrum"]="smallest magnitude"
eigen_min, _ = compute_eig0(eigenSolver)

eigenSolver.parameters["spectrum"]="largest magnitude"
eigen_max, _ = compute_eig0(eigenSolver)
return abs(eigen_max / eigen_min)