[petsc-dev] petsc4py "A[i,j] += 1"

[Barry Smith]

    Victor wanting a nice A[i,j] += 1 in Python reminds me of the days  
when "operator overloading" was going to turn C++ into THE great  
language for numerical computing.

    For example, if x,y,z,w are vectors one could write w = a*x + b*y  
+ z; in C++ to represent what is an ugly loop in Fortran 77 or C. The  
problem is that operator overloading is always ONE operation at a time  
which makes any kind of optimization difficult/impossible. (Partially  
because lots of temporary vectors would be created in evaluating the  
line of code.)

    The "solution" to this problem was "expression templates" where  
the operator overloading actually produced a list of operations that  
were then triggered by the = operator. For example, the operator a*x  
didn't define the product rather it saved the a and x and the product  
that would then be evaluated later.

    Either getting expression templates to work for this was either  
impossible or very hard so we still don't have efficient C++ vector  
math libraries.

    Now for sparse PETSc matrices one can overload [i,j] to call  
MatGetValue() and maybe A[i,j] = x to overload MatSetValues() (not  
sure if you can do this) But given how += is managed in Python  
operator overloading cannot get A[i,j] += x to be  
MatSetValues(  ADD_VALUES).

     I gave up on operator overloading in 1993 because the usual way  
is too inefficient and expression templates too damn hard. I viewed  
operator overloading as the sirens of greek mythology; always tempting  
the sailors to crash their ship's on the island even though a part of  
the men's minds know it is nuts. The beauty of operator overloading so  
overwhelms our senses.

     So operator overloading is useless; what I want instead is a  
language where a series of operations can be defined (but not failing  
like C++ expression templates for mathematics), thus allowing  
beautiful code that is efficient.

     Is there a way that the python interpreter (or something) could  
be mucked with so that A[i,j] += 1 gets translated into MatSetValues()  
ADD_VALUES? People can do a lot of weird %^^ with Python. In other  
words, can we resist the Sirens again or will they tempt us with our  
impossible dreams?

    Barry



On Apr 14, 2010, at 6:59 PM, Lisandro Dalcin wrote:

> On 14 April 2010 17:57, Victor Eijkhout <eijkhout at tacc.utexas.edu>  
> wrote:
>>
>> On 2010/04/13, at 7:56 PM, Lisandro Dalcin wrote:
>>
>>> Not sure if I was clear enough. In short, supporting A[i,j] += v  
>>> means
>>> tmp=A[i,j] should return in tmp a "proxy" object implemented in  
>>> such a
>>> way that tmp+=v do call MatSetValues(i,j,v,ADD_VALUES).
>>
>> Tricky. But if you can make that work....
>>
>
> I can make it work, but then val = A[i,j] will no longer return a
> scalar entry, but that proxy beast. Not a big deal, provided that
> accessing matrix values is not so common, right? And anyway, I have to
> implement a working MatGetArray() call and do the special case for AIJ
> or dense.
>
>
>>> What do you think? I really need some feedback to take an action on
>>> this issue...
>>
>> It's not terribly elegant, but I could live with having a  
>> matsetvalues call.
>> That doesn't seem to exist right now, correct?
>>
>
> Hey!! There are A LOT of Mat.setValues calls:
> http://code.google.com/p/petsc4py/source/browse/src/PETSc/Mat.pyx#603
>
> Unfortunately, all that is undocumented... Hope you understand that I
> do not have enough time/funding to spend on this... I do not even have
> feedback from users! Perhaps because they are happy, or perhaps
> because few people is using petsc4py for serious work out there.
>
>>> Victor, would you enter a discussion about this issue via private  
>>> mail
>>> or petsc-dev ?
>>
>> Barry seems to be interested in keeping up with petsc4py  
>> developments, so
>> let's talk on petsc-dev, and I'll make sure to prefix  my subjects  
>> with
>> "petsc4py".
>>
>
> OK, I'll try to write to petsc-dev commenting on this A[i,j]+= v
> bussines and I'll CC.
>
> -- 
> Lisandro Dalcin
> ---------------
> CIMEC (INTEC/CONICET-UNL)
> Predio CONICET-Santa Fe
> Colectora RN 168 Km 472, Paraje El Pozo
> Tel: +54-342-4511594 (ext 1011)
> Tel/Fax: +54-342-4511169