Re: [Pysparse-users] reading the source code of pysparse

[Dominique O. <dom...@gm...>]

On Sat, Feb 11, 2012 at 19:15, Oz Nahum Tiram <na...@gm...> wrote:
> On Sat, Feb 11, 2012 at 6:57 PM, Dominique Orban
> <dom...@gm...> wrote:
>
> Hi Dominique,
> Thanks for the answer. A few more questions.
>
>
>> This is a macro that adjusts the call to various types of Fortran
>> compilers (this particular call is a call to the BLAS library, which
>> is written in Fortran). Some Fortran compilers add a trailing
>> underscore to symbols, some add two, some add none, etc.
>>
>
> So this macro is found in BLAS? not in pysprase itself ?

The function dnrm2 is in the BLAS. The macro F77 is in
pysparse/include/fortran.h

>> I wouldn't recommend coding bicgstab all over again. The only costly
>> operations in Bi-CGSTAB (and other Krylov-type methods) are vector
>> operations (most often, addition of vectors) and operator-vector
>> products (e.g., A*x or A.T*x). I would say that to speed things up,
>> you'll want to speed up your operator-vector operations; they are the
>> dominant cost.
>>
>> You can take a look at PyKrylov (https://github.com/dpo/pykrylov)
>> which contains a pure Python implementation of Bi-CGSTAB and allows
>> you to input your operator in different ways (a Pysparse matrix being
>> one of them). For instance, you could implement your operator in C or
>> in Cython and that should speed things up. I believe that is the way
>> to go.
>
> Why is pure python Bi-CGSTAB should be faster than the one writen in C
> contained in pysparse?

It shouldn't. I'm saying that there isn't much to gain by recoding it
in C because the dominant cost is in the operator-vector products.
That's what should be fast.

>> Of course, vector operations in PyKrylov could also be speeded up with
>> Cython. That's been on my list for a while.
>
> You mean doing stuff on numpy vectores instead of for loops? I didn't see
> so many of these in the code ...

Or simply declaring types.

> And one last thing, in pysparse solvers, I pass a Matrix, initial x vector,
> and b vector (representing my boundary conditions.) I was quite
> confused to see the example here:
> https://github.com/dpo/pykrylov/blob/master/examples/bmark.py
>
>    for KSolver in [CGS, TFQMR, BiCGSTAB]:
>        ks = KSolver( lambda v: A*v,
>                      #precon = dp,
>                      #verbose=False,
>                      reltol = 1.0e-8
>                      )
>        ks.solve(rhs, guess = 1+np.arange(n, dtype=np.float), matvec_max=2*n)
>
> you are not assigning the lambda v to anything? What does this do?

lambda is standard Python. The lambda function is stored inside the
object upon instantiation. This syntax lets you solve several systems
with the same operator and different right-hand sides, or from
different initial guesses.

> solveing is done only with rhs and guess ? is this because A is
> already contained
> in ks ? I'd be happy to see an example which does not confuse me so much :-)

What's so confusing about it?

-- 
Dominique