dolfin team mailing list archive

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

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.

>>>>> 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.

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.

>>> 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.

Garth

> --
> Anders