dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 02/01/11 06:47, Johan Hake wrote:
> Hello!
>
> As pointed out by Garth in a couple of commit messages lately we need to
> consider the slicing interface for vectors and matrices in PyDOLFIN.
>
> Now that the parallel interface are starting to get more and more robust
> (great work Garth!) we need to consider what:
>
>    v[6], v[4:], v[index_array] # returning a scalar and a vector
>    A[5,6], A[5:,6], A[5:,6:] # returning a scalar, a vector and a matrix
>
> should return for serial and parallel structures.
>
> I guess Garth comments that the Matrix slicing "is just a toy" has some merits
> as it uses a dense representation of the Matrix (even though it is erasing
> zero entries). Regardless if we decide to discard the Matrix slicing or not we
> need to decide how off process indices should be treated.
>
> As it is now, a user need to check whether an index is owned by the local
> processor or not, by using GenericVector::owns_index, before he access the
> data.
>
> An other scenario would be to always allow a user to write:
>
>    v[5], v[4:] and v[index_array]
>
> and then return for example None (for single value access) and a reduced
> vector (only including local values) for the two other alternatives. Then the
> user need to deal with the result knowing how we treat this in parallel.

Another option is to not return a numpy vector, but to add a function to 
the C++ interface that creates a new reduced vector (which, in the 
general case, will be distributed across processes) based on a list of 
indices.

> Matrix slicing is more tricky especially now when GenericMatrix::resize is
> removed. There are probably ways of getting around this by building a
> SparsityPattern based on the slice. But I am not sure it makes any senses to
> do a matrix slice for a parallel structure, as it is only the rows that are
> local. Then we need to decide if the Matrix slicing should be:
>
>    1) removed

I think that for the fanciest usage cases this would be best.

>    2) reimplemented but disabled in parallel

Definitely no! A real selling point is that users can write simple 
Python programs that can run in parallel without knowing it. The code in 
the Python adaptivity module is broken now because it uses code that 
wouldn't work in parallel. We should make it hard for users to write 
code that won't work in parallel.

>    3) reimplemented for both serial and parallel cases using for example
>       a similar approach as sketched above for vectors, that is always
>       returning something.

We should be able to extract rows from matrices as distributed vectors 
for most backends. It looks like PETSc can also extract sub-matrices 
(and return them as distributed matrices). I guess all we need to do in 
the SWIG interface layer is convert the '[:,2]' part of A[:,2] into a 
list of indices it to pass to the C++ interface. Extracting sub-vectors 
from C++ should be easy.

I summary, I think that the best option is to add some more 
functionality to the C++ interface, and allow Python users to call these 
functions using NumPy-like syntax. This would provide most of the 
functionality that we had up to a few days ago. We should avoid using 
NumPy data structures since a DOLFIN vector/matrix should be presumed to 
be distributed, and NumPy objects can only be local.

For advanced users, would could make it easy to construct SciPy sparse 
matrices form the sparse data arrays. A user would then be on their own 
when it comes to parallel programs.

Garth

> Best regards,
>
> Johan
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp