MATLAB Left Matrix Multiply Equivalent

Question

Hi all,

I have to compute the supremizer operator T defined as

(Tw,v)_X = a(w,v)

where (w,v)_X is the inner product and a(w,v) is the bilinear form of my problem.

Now, let A and B be the matrices associated do the bilinear form and the inner product in my FE problem.

T can be computed as follows

T = B^(-1) A

In MATLAB it can be done by typing

T = B \ A

Now the question: Is there a way to perform the left matrix multiply in Fenics? Or, is there an efficient way to get the inverse of B, maybe solving a system like
B X = I
where I is the identity matrix? What is the python code?

Thanks in advance

Answer 1

Hi, here's a brute force solution that uses scipy.sparse.

from dolfin import *
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import inv
import matplotlib.pyplot as plt

# Choose backend that allows .data()
parameters.linear_algebra_backend = 'uBLAS'

mesh = UnitSquareMesh(30, 30)
V = FunctionSpace(mesh, 'CG', 1)

u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
m = inner(u, v)*dx
L = Constant(0.)*v*dx

bc = DirichletBC(V, Constant(0), DomainBoundary())
A, _ = assemble_system(a, L, bc)
M, _ = assemble_system(m, L, bc)

# Prepare A and M for scipy
rows, cols, values = A.data()
A = csr_matrix((values, cols, rows))

rows, cols, values = M.data()
M = csr_matrix((values, cols, rows))

print 'Inverting %d by %d matrix' % (V.dim(), V.dim())
timer = Timer('inverse')
timer.start()
T = inv(M).dot(A)
timer.stop()
print '\t\t done in %g' % timing('inverse')

# Plot
plt.figure()
plt.subplot(131); plt.spy(A);
plt.subplot(132); plt.spy(M);
plt.subplot(133); plt.spy(T);
plt.show()