Set values to an assembled matrix

Question

Hello,

I have a very basic question.

After the assembly of a linear system with matrix A and right hand side vector b, I want to insert changes in the matrix A.

To test the functions and methods I wrote this example, based on suggestions I found in this forum: https://answers.launchpad.net/dolfin/+question/190367

from dolfin import *
import numpy
set_log_level(ERROR)
parameters["linear_algebra_backend"] = 'PETSc'

# Create mesh and define function space
mesh = UnitSquareMesh(2,2)
V = FunctionSpace(mesh, 'Lagrange', 1)

# Define boundary conditions for initial guess
tol = 1E-14
def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol

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)
a = inner(nabla_grad(u), nabla_grad(v))*dx
f = Constant(0.0)
L = f*v*dx
A, b = assemble_system(a, L, bcs)
u_k = Function(V)
solve(A, u_k.vector(), b, 'lu')

print A
block = numpy.ones(2,dtype=numpy.float_)
print block
rows = numpy.array([0],dtype=numpy.uintc)
print rows
cols = numpy.array([1,2],dtype=numpy.uintc)
print cols
A.set(block,rows,cols)

When executed, I have the following error:

<Matrix wrapper of <PETScMatrix of size 9 x 9>>

[ 1.  1.]
[0]
[1 2]

Traceback (most recent call last):
  File "set_matrix_ex.py", line 39, in <module>
    A.set(block,rows,cols)
TypeError: contiguous numpy array of 'dolfin_index' expected. Make sure that the numpy array is contiguous, with 1 dimension, and uses dtype=intc.

Apparently there is a problem with the data type. I read all documentation I found and I cannot understand what I am doing wrong.

By the way, the matrix is of type "PETScMatrix" but I also try the "uBLAS" backend with the same error.

Thanks for your help!
Diego

Answer 1

First and most important: you should not manipulate the values of an assembled matrix. It will almost always be very slow.

Second, to get your example running, just change from uintc to intc as suggested in the error message.

You should also instantiate an Assembler object and set the parameter finalize_tensor to False (assembler.finalize_tensor = False), then assemble (assembler.assemble(...)) and finally call A.apply().

Comments

Thanks for your answer. That solved my question!
Regarding the computational time, I understand your point.

By assuming that I will not change the sparsity pattern of the matrix when setting the new values, it would be the process still very slow?

Another option I am considering is to copy the values of the assembled matrix to a new PETScMatrix, but considering the appropriate modifications. Do you think this could significate in some time improvement?

Thanks again!

---

If you don't change the sparsity pattern and you appropriately use finalize_tensor = False in combination with A.apply(), there shouldn't be any performance penalty (except for the overhead of Python function calls).

---

Thanks!
That solved my question....again.

---

If I want to change only some values of the petsc matrix?
For example

a21_lump = c*w*dx 
mass_matrix21 = assemble(a21_lump)
mass_matrix21.zero()
mass_form21= action(a21_lump, const)
vec=assemble(mass_form21)
print vec
N = len(vec.array())
for i in range(1, N-1): 
     mass_matrix21.array()[i][i-1]= vec.array()[i]

If I print, no changes occurred.
What should I do?