dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Jan 19, 2011 at 11:08:01AM +0100, Anders Logg wrote:
> On Wed, Jan 19, 2011 at 09:48:07AM +0000, Garth N. Wells wrote:
> >
> >
> > On 19/01/11 09:21, Anders Logg wrote:
> > > On Tue, Jan 18, 2011 at 07:55:30PM +0000, Garth N. Wells wrote:
> > >>
> > >>
> > >> On 17/01/11 18:16, Anders Logg wrote:
> > >>> On Mon, Jan 17, 2011 at 06:07:19PM +0000, Garth N. Wells wrote:
> > >>>>
> > >>>>
> > >>>> On 17/01/11 18:03, Anders Logg wrote:
> > >>>>> On Mon, Jan 17, 2011 at 04:30:23PM +0000, Garth N. Wells wrote:
> > >>>>>>
> > >>>>>>
> > >>>>>> On 17/01/11 16:28, Anders Logg wrote:
> > >>>>>>> The following doesn't seem to work in Python any longer:
> > >>>>>>>
> > >>>>>>> A = Matrix(10, 10)
> > >>>>>>>
> > >>>>>>> Is matrix initialization broken?
> > >>>>>>>
> > >>>>>>
> > >>>>>> Probably not broken - more like
> > >>>>>>
> > >>>>>>   A = Matrix(10, 10)
> > >>>>>>
> > >>>>>> is not supported. Since a Matrix is in general sparse, it doesn't make
> > >>>>>> sense to initialise it as above (some backends even insist on the
> > >>>>>> sparsity being defined at construction).
> > >>>>>>
> > >>>>>> Garth
> > >>>>>
> > >>>>> It would be good to allow simple initialization of a Matrix without
> > >>>>> requiring to go through all the hassle of creating a SparsityPattern.
> > >>>>>
> > >>>>> I generally don't think we should disallow certain operations just
> > >>>>> because they are potentially slow.
> > >>>>
> > >>>> It's not just that they're slow. They cannot be supported by all backends.
> > >>>
> > >>> Then we throw an error.
> > >>>
> > >>
> > >> Better still, use an appropriate data structure ;). If a user wants a
> > >> dense matrix, we should make it easy to use a dense matrix.
> > >>
> > >> The Matrix class in DOLFIN is sparse and the premise should be that it's
> > >> distributed. Being consistent in this makes life a lot easier in parallel.
> > >
> > > Yes, but being sparse does not mean it should be impossible to assign
> > > individual entries. This can be very helpful in debugging (small)
> > > problems.
> > >
> >
> > If it's small, use a dense matrix. If you want to add entries use
> > GenericFoo::set_local. Everything is there.

Even I don't know how to do that (use a dense matrix which supports
index assignment + DOLFIN assembly) and expect most users don't. Is it
even possible?

Most users (I think, at least I do) would expect to be able to do

A = Matrix(10, 10)
A[0, 0] = 1.0

just to try out the Matrix class, or to just to manually set a
boundary condition.

> > >>>>> Some operations (like Matrix index
> > >>>>> access) will only be performed for toy problems or while testing and
> > >>>>> then speed is not very important (since the problem is small anyway).
> > >>>>>
> > >>>>
> > >>>> If it's made available, it will be used. This type of operation is main
> > >>>> reason for things breaking in parallel. We have getitem and setitem
> > >>>> which are sufficiently unfriendly that hopefully users get the idea that
> > >>>> they're for testing only.
> > >>>
> > >>> Just because something can be misused, it shouldn't be disallowed.
> > >>
> > >> Doesn't sound like a strong argument. It should not be easy to misuse a
> > >> library. That's why I can live with get/setitem, but not indexing into a
> > >> sparse distributed matrix via A(i, j).
> > >
> > > Everything can be misused. We still allow overloading eval for
> > > Expressions in the Python interface. It's very slow compared to
> > > jit-compiling an Expression, but it's very convenient to have.
> > >
> >
> > That is not misuse.
>
I don't see the difference.

> > Speed is not significant for many applications. Performance is only
> > secondary in my argument. The main reason for not permitting operations
> > (or making them difficult) is that they are inherently unsuited to the
> > data type and will lead to code breaking. We should instead provide
> > suitable data types. There is a reason why established parallel
> > libraries like PETSc and Trilinos do not facilitate very easy element
> > access.
>
I don't think comparisons with the PETSc/Trilinos interfaces are very
relevant. We make a big effort to provide a simpler/cleaner interfaces
to linear algebra.

I also don't understand why it will lead to code breaking. If a
warning is issued when getitem/setitem are used, then we've done all
we can. It would also be easy to check inside those functions that all
conditions are fulfilled for using them, like not running in parallel.
Extra checks inside those functions will lead to some overhead but
that's not a problem since they should not be used for
performance-critical operations anyway.

> > >>> It would be easy (and more helpful) to add a message the first time
> > >>> getitem/setitem is used in a program:
> > >>>
> > >>>   Warning: Index access to matrix/vector values is potentially very
> > >>>   slow and it breaks in parallel. To disable this warning, set the
> > >>>   parameter "warning_index_access" to false.
> > >>>
> > >>
> > >> Yes, I would like to see a set/getitem warning.
> > >
> > > Yes + operator[] mapping to those. I think that would be a good
> > > compromise, don't you think? ;-)
> > >
> >
> > Having cleaned up the parallel linear algebra side, I don't think that
> > it is a compromise. We should encourage the use of suitable data types.
> > We need to decide if a Matrix is dense or sparse, and it it's local or
> > distributed. If we want to run in parallel, the premise must be that a
> > Matrix is distributed.
> >
> > To make things work in parallel, as developers we need to make more
> > effort to bear in mind that objects are distributed, and that this makes
> > certain operations more complicated.

I agree, but I don't see any conflict between that and mapping
operator[] to getitem/setitem:

1. The use would be discouraged through a warning.

2. An informative error message would be issued when their use is illegal.

--
Anders