[petsc-dev] petsc4py "A[i,j] += 1"
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
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
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
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
On Apr 14, 2010, at 6:59 PM, Lisandro Dalcin wrote:
> On 14 April 2010 17:57, Victor Eijkhout <eijkhout at tacc.utexas.edu>
>> 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
>>> 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.
>>> Victor, would you enter a discussion about this issue via private
>>> 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
> 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
More information about the petsc-dev