IMPOSSIBLE operations on matrices objects : A*B, A.T !!!!!!!!!!!

Question

as a matter of fact my problem seems to be trivial
I spent the whole day searching on FENICS book and INTERNET to solve this issue

R = assemble(nabla_grad(u),nabla_grad(v)*dx)

I need to compute this : R*R and R.T

the only thing I found I can do is : R+k*R / k is a real number

It would be nice if someone could bring any help

Answer 1

Hi, the resulting vector from operations that you require can always be obtained without form the matrices explicitly. For this the linear algebra API of DOLFIN is sufficient. Should you require the matrices explicitely you can use PETSc API though petsc4py. Both usecases are illustrated below

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

e = Constant((1, 2))
A = assemble(inner(dot(grad(u), e), v)*dx)
B = A.copy()
A_mat, B_mat = as_backend_type(A).mat(), as_backend_type(B).mat()

x = interpolate(Expression(('sin(x[0])', 'cos(x[1])'), degree=2), V).vector()

#--- Transpose
# By action
y1 = x.copy()
A.transpmult(x, y1)  # y1 = A'*x

# With explicit matrix B = A'
A_mat.transpose(B_mat)

# We have non-sym matrix so sanity check
A -= B   # A = A - B
e = A.norm('linf')
info('|A-B|=%g' % e)

# Make sure they are the same
y2 = x.copy()
B.mult(x, y2)  # y2 = B*x

y2.axpy(-1, y1)   # y2 = y2-y1
e = y2.norm('linf')
info('|y1-y2|=%g' % e)

#--- Power
# By action, y1 = A*(A*x)
for i in range(2): A.mult(x, y1) 

# With explicit matrix
C = A.copy(); C_mat = as_backend_type(C).mat()
C_mat.matMult(C_mat)  # C = C**2
C.mult(x, y2)

# Make sure they are the same
y2.axpy(-1, y1)
e = y2.norm('linf')
info('|y1-y2|=%g' % e) 

Comments

Thank you MiroK for the details. This was worth accepting it as best answer.

Answer 2

Those operations are available in the linear algebra backend.
For example, with PETSc and petsc4py,

R = assemble(a, tensor = PETScMatrix())
RT = PETScMatrix(R.mat().transpose(R.mat().duplicate()))
R2 = PETScMatrix(R.mat().matMult(R.mat()))

Comments

Thank you Nordaas for the answer. Still just some things to fix
how to transpose a rectangular matrix
and one question last, how to form a block matrix from these element petsc Matrices ?

Let say for example I want this [[R,0],[0,RT]]

Best

---

For non-square matrices, you can construct the transpose matrix as follows.

from petsc4py.PETSc import Mat
RT = PETScMatrix(R.mat().transpose(Mat()))

Constructing block matrices is more involved. You can use PETSc MatNest or implement the block structure in Python. See also cbc.block.