GenericMatrix.set

Question

Hello,

Even though there are already some posts on the GenericMatrix.set function, things are still unclear to me...

I'd like to construct a new GenericMatrix as the matrix product of massMatrix.transpose * stiffnessMatrix, and then solve.
The point is I really don't know how to proceed.
Any help would be appreciated :)

The code below is just a minimal example of what I'd like to do

---

from dolfin import *
import numpy

mesh = UnitInterval(10)
V = FunctionSpace(mesh, 'CG', 1)

# Define boundary conditions for initial guess
def left_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)

def right_boundary(x, on_boundary):
    return on_boundary and near(x[0], 1)

Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
bcs = [Gamma_0, Gamma_1]

# Define variational problem ans solve
u = TrialFunction(V)
v = TestFunction(V)
Mass = assemble(u*v*dx)
Stiff = assemble(u.dx(0)*v.dx(0)*dx)
L = assemble(Constant(0.0)*v*dx)

A=assemble(Constant(0.0)*u*v*dx) # define a new GenericMatrix
# I'd like to define
A.set(Mass.array()*Stiff.array())

for _bc in bcs:
    _bc.apply(A,L)

_U_ = Function(V)
solve(A, _U_.vector(), L)
plot(_U_, interactive = True)

Thanks in advance,
Nicolas

Edit: just to mention I use Fenics v. 1.0.0

Answer 1

Instead of solving
$$ M S x = L $$
you can first solve
$$ M y = L $$
and then
$$ S x = y .$$

Comments

Thanks a lot Jan,

I try this ASAP.

However, since I have to run the $MSx=L$ equation beneath a fixed point loop, is this way efficient in term of CPU cost?

Best regards,
Nicolas.

---

is this way efficient in term of CPU cost?

Sure. Moreover, if one of the matrices is constant, you could reuse its LU factorization or preconditioner.

---

All right.

Your approach seems to solve my problem. However, I let the subject open: I'm still curious on how to construct a dolfin matrix from a numpy array (obviously with consistent dimensions) ;)

---

I don't think it is a good idea to achieve a matrix multiplication through numpy interface. Consider that matrices are stored sparsely in some highly optimized, distributed data structures (for example using PETSc). So a numpy manipulation does not definitely scale well and also has some issues in parallel. I see three possibilities.

1. You achieve solving of $MSx=L$ by matrix-free method (using Krylov sequence method) so you don't have to calculate a product $MS$ but you rather define action of the operator $MS: x \mapsto MSx.$

2. You translate your product into the language of functions and forms and assemble $MS$ directly as a single matrix.

3. In the worst case you can perform a multiplication $MS$ using PETSc. You can do this from pyDOLFIN using compile_extension_mdoule or using petsc4py python module (if you compile development vesrsion of DOLFIN or upcoming 1.3.0 release with petsc4py support).

If you still insist on creating matrix from numpy array, this is not supported directly because of these reasons, I guess. But you possibly can (I don't how exactly) create a matrix with custom sparsity pattern (but how does it look like?) and fill it using numpy interface. Alternatively you can use uBLAS dense matrix to avoid a need for knowledge of a sparsity pattern.

I reccomend you strongly to install (or better compile) FEniCS 1.3.0 when it gets released if you want to use some of these tricks.

---

Thanks a lot, Jan.

That solved my question.