How to add linear combinations of FEM Functions Efficently

Question

In the code I'm trying to run there is a point where I solve a linear system and get an array of doubles of size N stored in a numpy array X. I then want to construct a new FEM function by using these values and a linear combination of already defined finite element functions as seen below.

for i in range(0,N):
solHold.vector()[:] = solHold.vector() + X[i]*POD_funcs[i].vector()

The issue I'm having is that this loop runs ridiculously slow. I have to 2 questions.

1. Why does what I'm trying to do above run so slowly.

2. What is an efficient manner in which I can perform this procedure?

Answer 1

The easiest way is to project everything into the same functionspace and then add all the degrees of freedom. Here I assume your numpy vector X is associated to some dolfin FunctionSpace X_function_space.

X_func = Function(X_function_space)
X_func.vector()[:] = X[:]
X_proj = project(X_func, solHold.function_space())

solHold.vector()[:] += X_proj.vector()[:]

This way you don't have to do a slow python for-loop. If by some chance X_function_space and solHold.function_space() are the same then you can just add the DOF vectors directly.

Answer 2

Hi,

The reason your code is running very slow is that you are creating many temporary vectors to store X[i]*POD_funcs[i].vector() and then to store the sum.

Depending on how large is N, this solution may be better:

for i in range(N):
    solHold.vector().axpy(X[i], POD_funcs[i].vector())

note: u.axpy(alpha, v) computes u = u + alpha v without creating temporary vectors.

If this is still slow, you will most likely have to implement the loop using the C++ interface of FEniCS.