# Fast computation of bilinear forms

Hi,

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)`.

How many `j`s and `k`s 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 `b`s since this is an optimization algorithm for `q`.