inner product (assemble) of set of states

Question

I need to compute the scalar product of each pair of a set of states. So far, what I've been doing is

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

The states Theta are here organized in a Python list. The whole thing could probably be sped up by fitting Theta in some sort of Multivector type, but if I recall correctly FEniCS doesn't support any of that.

Any other ideas of how to speed this up?

Answer 1

Use as_tensor or as_matrix to build tensor ufl epression out of Theta. Then check thread Integrate expression of rank 1 with free indices () on FEniCS mailing list. (The approach can be generalized to integrate expression of arbitrary rank.)

Comments

Hm, not quite there yet. I did

tv = as_vector(Theta)
B = assemble(outer(tv, tv) * dx)

and get

ufl.log.UFLException: Trying to integrate expression of rank 2 with free indices ().

Here's where the trial and test functions of "R" are supposed to come in, but clearly the "mesh" that this function space is supposed to live in is purely artificial. I wouldn't know how to chose it sensibly.

---

Obviously you don't make a mistake by taking a mesh of tv. Note that by

R = VectorFunctionSpace(mesh, "R", 0, dim=n)
i, j = TrialFunction(R), TestFunction(R)
assemble(dot(tv, i)*dot(tv, j)*dx)

you get a Matrix with no relation to some mesh. This is what you want.