python matrix construction and multiplication

Question

I would like to construct a PETScMatrix from a numpy.array and thus try something like

A, N = PETScMatrix(), 10
A.resize(N,N)
values = numpy.ones(2,dtype=numpy.float_)
rows = numpy.array([0],dtype=numpy.intc)
cols = numpy.array([1,2],dtype=numpy.intc)
A.set(values,rows,cols)
A.apply()

First question: how is resize used correctly?
Next, for a specific optimisation algorithm, the following expression B has to be evaluated

A, alpha = assemble(u*v*dx), 0.3
B = A*A-alpha*I

where I is the identity matrix. How can this be implemented?
Thanks for advice, Martin

Answer 1

PETScMatrix is a sparse matrix, so requires more than just its dimensions initialise it. It needs a sparsity pattern for intialisation. Building a sparsity pattern can be complicated, so it's best to let DOLFIN take care of this, if possible. For example, assemble will initialise a matrix with the appropriate sparsity pattern.

Comments

My comment regarding the use of uBLAS backend is to avoid conflicts with the Initialization of the PETSc ambient. If you set the PETSc backend, PETSc is initiate by FEniCS and also when the first call to the pets4py interface is made, giving an error.

I think it is solved if FEniCS is compiled with petsc4py support.

---

The file "demo_extrension.py" is just the sample from the documentation of compile_extension_module, i.e.,

from numpy import arange
from dolfin import *
code = '''namespace dolfin {
  void PETSc_exp(boost::shared_ptr<dolfin::PETScVector> vec) {
  boost::shared_ptr<Vec> x = vec->vec();
  assert(x); VecExp(*x); } }'''
ext_module = compile_extension_module(code, 
         dolfin_module_import=["la"], 
         additional_system_headers=["petscvec.h"])
vec = PETScVector(10)
vec[:] = arange(10)
print vec[-1]
ext_module.PETSc_exp(vec)
print vec[-1]

Will try the other example but this should work nevertheless, shouldn't it?

I see your point DiegoM but I find it rather unsatisfying to work around the problem like that. However, your remarks are appreciated!

---

The file "demo_extrension.py" is just the sample from the documentation of compile_extension_module

You need to refer to development version of documentation for use with development version.

---

Good point. That works.
Well, I assume I will try to handle things [i.e. PETSc objects] in extension modules since no better options seem to be available.

---

no better options seem to be available

Well, using PETSc directly is definitely the best option as DOLFIN is not a complete interface to PETSc, it is rather a FEM assembler.