How to compute the condition number efficiently?

Linear Algebra
2

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