Vector of unkowns

Question

Dear all,

I am reading the documentation about how to solve a multivariable problem, defined as $ A u = L$, where $u$ and $L$ are not single vectors, but matrices. I'd like, in other words, to solve for different values. I can do this in Wolfram Alpha, for instance:

$$
A = (\, ( -1, \; 1 ), \, (\, 0,\, 1 )\, )
$$
by rows (See the note below),

$$
x = \left( x_1, \; x_2 \right)\,, u = \left( u_1, \; u_2 \right)\,,
$$

and

$$
L_x = \left( 3, \; -1 \right)\,, L_u = \left( 1, \; 0 \right)\,,
$$

Can be rewritten easily as $A X = L$, with $A$ as before (obviously), and $X$ and $L$ being matrices whose columns are $X = \left( x, \; u \right)$, and $L = \left( L_x, \; L_u \right)$. In Wolfram Alpha is here only for x, only for u, and finally the whole system at once.

I'm not sure how I can implement this in fenics, can you help me in figuring out the right way to approach the problem?

Thanks & Cheers!

PS. I cannot write a matrix here:
$$
A =
\left(
\begin{matrix}
-1 & 1 \
0 & 1
\end{matrix}
\right)\,,
$$
I've tried array, matrix, \matrix, but none of those produced a matrix, instead, they all put the rows one after the other, no new line added. weird...

Answer 1

If you are solving small dense linear systems, then you can use numpy directly:

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html

Concatenating arrays is very easy in numpy as well.
If A is a dense n by n matrix and Lx and Lu are 2 vectors of size n, then you can write:

import numpy as np
L = np.zeros(n,2)
L[:,0] = Lx
L[:,1] = Lu
X = np.linalg.solve(A, L)
x = X[:,0]
u = X[:,1]

NOTE: From a dolfin Vector you can extract a numpy array with the method array.

(For further help with numpy, use one the numpy forums :) )

Comments

Yes, except that the problem I have is sparse and big :)

It comes from a FEM formulation so, I'd like to use assembly and solve...

---

Krylov methods for multiple RHS are an active area of research.
I don't think that PETSc supports multiple RHS, but Belos (in Trilinos) does. I think it should not be too hard to integreate the multiple RHS solvers in Belos if you are willing to write a cpp module for the FEniCS instant compiler.

As alternative you could use some direct solver in FEniCS (PETScLUSolver, or MUMPSLUSolver so you can reuse the factorization for multiple RHSs.

---

Thanks for the pointers!

The Belos way is longer, but I think feasible (I used Trilinos several times), but the last two appeal to me.

However, PETScLUSolver won't work in MPI, and I'd like to use a cluster. The MUMPSLUSolver seems to be unavailable... I get the following error:

error: unknown type name 'MUMPSLUSolver'
        MUMPSLUSolver solver(A);

So for the moment I'm stuck with PETScLUSolver, and no MPI. That's not good :(