Fast computation of bilinear forms


I have the bilinear form q that can be represented as q(f, g) = int2d(Th)( b * ( dx(uf) * dx(ug) + dy(uf) * dy(ug) ) ) where b is a fixed function, uf and ug solve a given BVP (details of the specific problem are not important) with Neumann boundary conditions f and g, respectively.

I also have a collection of functions f1,...,fn and I want to compute the matrix B given by B_jk=q[fj,fk] for all pairs j,k.

The question is, what is the fastest way to do it?

Using the expression for q each time is extremely slow, so what I am doing now is to compute the matrix A = q(Vh, Vh) and computing the inner products <A * f[], g[]> whenever I want to compute q(f, g).

Thanks in advace!

How many js and ks do you have? Most often, you do not want to store this explicitly, and want to deal with this in a matrix-free fashion, or with another representation, like Kronecker product.

Give or take 40, I don’t think they are the problem. This algorithm fits into a larger process where this has be done a lot of times for different bs since this is an optimization algorithm for q.