Add Matrix to Assembled Matrix

1

Hello,

I am having trouble modifying the assembled matrix.

Specifically:

  1. I have an assembled matrix A (A = assemble(…)).
  2. I have another matrix B (that I choose)
  3. I want to get A + B

I have seen the functions ‘add’ and ‘add_local’ of the GenericMatrix class in the reference documents.
> add(block, num_rows, row)
> add_local(block, num_rows, row)

However, I am unfamiliar with the parameters and descriptions. For instance, I don’t know what these are: “block”, “global indices”, “local indices”

Also, should the rows be associated with the ‘dof_to_vertex_map(V)’?

If someone can shed some light on this through an example, it will be much appreciated.

Thanks,
Trevor

1 Like
2

If I need to edit values of an assembled matrix, I usually go through the petsc4py API. See the following example:

from dolfin import *

mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh,"Lagrange",1)
A = assemble(TrialFunction(V)*TestFunction(V)*dx)

# If by "+ B" you mean adding numbers directly to specific values of A, I
# usually do this through petsc4py:
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc

# To add to existing values:
ADD_MODE = PETSc.InsertMode.ADD

# To replace values instead:
#ADD_MODE = PETSc.InsertMode.INSERT

print(A.array())  # Original A

# (Note: Printed output is somewhat unclear in parallel)

# What value to add/insert at matrix entry i,j (i.e., the definition of B):
i = 0
j = 1
value = 1.0

Am = as_backend_type(A).mat()

# If the value you want to modify is not allocated as a nonzero, you need to
# set this option (with some cost to performance).  Ideally, you would
# set up a matrix with the necessary space allocated, assemble into that one,
# and then edit values without this.
Am.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)

# In parallel, you want to make sure a given value is added only once, ideally
# on the process that owns that row of the matrix.  (Can skip the check for
# ownership range in serial.)
Istart, Iend = Am.getOwnershipRange()
if(i<Iend and i>=Istart):
    Am.setValue(i,j,value,addv=ADD_MODE)
Am.assemble()

print(A.array())  # A with 1.0 added to the i,j entry

A complete listing of petsc4py functions and arguments can be found here:

https://www.mcs.anl.gov/petsc/petsc4py-current/docs/apiref/index.html

(You can usually figure out what all the functions/arguments mean by looking up the regular PETSc documentation, since the function names are similar.)

1 Like
3

Hello, I tried your solution and it works. However, it is extremely slow. Have you experienced this slowness?

4

In the applications where I’ve used this, the limiting factor in performance was computing the values that I was adding in Python. One optimization might be to try setting many values at once with the setValues method instead of setValue. In setValues, you can pass numpy arrays of indices and values.

5

If you are adding an entire matrix instead of a few values it is probably better to use something like this (through the petsc4py Mat objects of the matrices) https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatAXPY.html

The PETSc documentation is decent, the petsc4py doc is almost non-existent (last I checked), so you have to translate the C code examples to Python yourself

Also, any reason why “C = A + B” does not work? I think that should create a new PETSc matrix for you, but I have not checked.

1 Like
6

Hello, I tried using the PETSc Mat objects but, as Kamensky mentioned, the setValues method is extremely slow.

Do you guys know how to use the add(…) functions of the GenericMatrix class? As I understand, these are built-in to FENICS and should be faster.

Thanks for the help!

Here’s what I’m currently doing (for reference):

**B_py is the matrix that I want to add (stored as a Numpy array)

# Convert to B_py PETSc matrix
B = PETSc.Mat().create()
B.setSizes([len(nodes),len(nodes)])
B.setType("aij")
B.setUp()
B.setValues(range(0,len(nodes)), range(0,len(nodes)), B_py)
B.assemble()

# Convert A (stiffness matrix) to PETSc matrix
S = PETSc.Mat().create()
S.setSizes([len(nodes),len(nodes)])
S.setType("aij")
S.setUp()
S.setValues(range(0,len(nodes)), range(0,len(nodes)), A.array())
S.assemble()

# Add the matrices
M_pet = S.__add__(B)
M = dolfin.PETScMatrix(M_pet)

solve(M == ...)
7

How large is your matrix? If you are creating a PETSc matrix directly from a numpy matrix then I assume you end up with a dense matrix stored on sparse Aij format (very non optimal)

If your matrix is large and sparse you should only insert the non-zero entries. I am not sure that setValues does this by default, I do not think so. If your matrix is large and dense then I do not know if PETSc is the correct choice. At least you should probably not use Aij format, try B.setType(“dense”)

Please surround your code with triple backticks ``` like this

# This is a comment and not a markdown heading
B.setType("dense")
1 Like
8

My matrix is 1000 by 1000 right now.

Changing B and S to dense matrices fixed the slowness of the assembly.
However, now when I call:

solve(M, u.vector(), b)

I get the following error:
Error: Unable to successfully call PETSc function ‘KSPSolve’
Reason: PETSc error code is: 92

Could this be because it doesn’t accept dense matrices?

9

PETSc solvers are all built for sparse matrices as far as I know. If you are working with dense matrices then LAPACK and similar libraries made for dense matrices should be used. Luckily 1000x1000 is very small, so you can just run numpy.solve which uses a dense solver (maybe LAPACK, depends on how numpy was built).