dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Sunday January 2 2011 14:51:37 Garth N. Wells wrote:
> On 02/01/11 22:45, Johan Hake wrote:
> > On Sunday January 2 2011 13:08:05 Garth N. Wells wrote:
> >> 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.
> > 
> > This sounds good. If the slice came as an Array<uint>  it would be easy
> > to map this to a NumPy like interface.
> > 
> >>> 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.
> > 
> > We then need a more elaborated getrow method in GenericMatrix.
> 
> The present getrow is a very troublesome function in parallel.
> 
> > Btw, is it possible to create (or initialize) distributed vectors with
> > arbitrary distribution patterns?
> 
> Almost. There is a new constructor that takes the local range for the
> vector. It must be linear, but it is possible to specify the number of
> entries on each process.

Ok

> >If so it would be cool to support this
> >
> > through the GenericVector interface, so we could implement a sub vector
> > extension in the GenericVector interface.
> > 
> >> 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.
> > 
> > Ok, do you know if Epetra supports this too?
> 
> I don't think so. Epetra has very limited matrix functionality compared
> to PETSc.

Ok

> >> 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.
> > 
> > Ok, what would this syntax be?
> > 
> >    GenericMatrix::getsubmat(const Array<uint>  rows,
> >    
> >                             const Array<uint>  columns,
> >                             GenericMatrix&  values);
> >    
> >    GenericMatrix::setsubmat(const Array<uint>  rows,
> >    
> >                             const Array<uint>  columns,
> >                             const GenericMatrix&  values);
> >    
> >    GenericMatrix::getrow(uint row, const Array<uint>  columns,
> >    
> >                          GenericVector&  values);
> >    
> >    GenericMatrix::setrow(uint row, const Array<uint>  columns,
> >    
> >                          const GenericVector&  values);
> >    
> >    GenericVector::getsubvec(const Array<uint>  indices,
> >    
> >                             GenericVector&  values);
> >    
> >    GenericVector::getsubvec(const Array<uint>  indices,
> >    
> >                             GenericVector&  values);
> 
> Something like that.
> 
> > For scalars I guess we just return the global value?
> > 
> >> 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.
> > 
> > The indice vectors need to be global, and can then be mapped to NumPy
> > arrays, which can then be mapped to slice objects. FYI, previously we
> > returned GenericVectors and GenericMatrices not NumPy arrays.
> 
> True. I was thinking about the 'array' function.

Ok, for now the 'array' method return a dense local structure, which is 
probably the correct behavour. This might also be usefull for debugging.

> >> 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.
> > 
> > Yes that would be easy to do. For now this is only possible for the uBLAS
> > and MTL4 backends, using the GenericMatrix::data method. Is this data
> > available for the distributed backends too?
> 
> Yes, but for simplicity it would be advisable not to mess with it.

Agree.

Johan