compute Gramian matrix

Question

Given a number of functions $\theta_i$, I would like to compute the Gramian matrix
$$
G_{i,j} = \int_\Omega \overline{\theta_i} \theta_j.
$$
Up until now, I just did something along the lines of

n = len(Theta)
G = np.empty((n, n))
for i in range(n):
    G[i, i] = assemble(Theta[i]*Theta[i]*dx)
    for j in range(i+1, n):
        alpha = assemble(Theta[i]*Theta[j]*dx)
        G[i, j] = alpha
        G[j, i] = alpha.conjugate()

This can take quite long though, and I thought about speeding this up a bit. Ideally, one would formulate all those assemble()s as
$$
\int_\Omega \overline{\theta_i} \theta_j = \alpha_i M \alpha_j,
$$
where $\alpha_i$ are the coefficients of the state $\theta_i$ in the finite-element space, and $M$ is the corresponding mass matrix. The Gramian can thus be expressed as
$$
G = A^H M A
$$
where the $i$th column of $A$ is $\alpha_i$. This representation encourages the use of a DGEMM-like routine. Not sure if it exists in this form (e.g., if the columns are PETSc vectors), let alone if anything is made available through the Dolfin/UFL layer.

When using uBLAS, I could just arrange the above in NumPy (and use its dot() which calls DGEMM), but I'm not sure how to do M*A (column by column, I suppose) in Dolfin. -- How does GenericLinearOperator::mult (http://fenicsproject.org/documentation/dolfin/dev/cpp/programmers-reference/la/GenericLinearOperator.html#GenericLinearOperator::mult__GenericVectorCR.GenericVectorRC) work in Python?

Any hints, or other ideas for speeding up the above code?

Answer 1

Okay, I found a good solution for when the uBLAS backend is used, and still using something slower in all other cases.

def _build_gramian_matrix(Theta):

    # Build mass matrix.
    Q = Theta[0].function_space()
    u = TrialFunction(Q)
    v = TestFunction(Q)
    M = assemble(u*v*dx)

    n = len(Theta)

    # Build list of M*Theta[i].
    MTheta = []
    for i in range(n):
        MTheta.append(Function(Q))
        M.mult(Theta[i].vector(), MTheta[i].vector())

    if parameters['linear_algebra_backend'] == 'uBLAS':
        # Cast Theta and MTheta into NumPy arrays.
        A = np.empty((Theta[0].vector().size(), len(Theta)))
        for i, theta in enumerate(Theta):
            A[:, i] = Theta[i].vector().array()
        MA = np.empty((MTheta[0].vector().size(), len(Theta)))
        for i, theta in enumerate(Theta):
            MA[:, i] = MTheta[i].vector().array()
        G = np.dot(A.transpose().conjugate(), MA)
    else:
        # Compute all the inner products manually.
        G = np.empty((n, n))
        for i in range(n):
            print(i)
            G[i, i] = Theta[i].vector().inner(MTheta[i].vector())
            for j in range(i+1, n):
                alpha = Theta[i].vector().inner(MTheta[j].vector())
                G[i, j] = alpha
                G[j, i] = alpha.conjugate()

    return G