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

Barry Smith bsmith at mcs.anl.gov
Sat Apr 17 21:03:10 CDT 2010

```    Victor wanting a nice A[i,j] += 1 in Python reminds me of the days
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
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

I gave up on operator overloading in 1993 because the usual way
is too inefficient and expression templates too damn hard. I viewed
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.

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

```