[Numpy-discussion] numpy and roundoff(?)

[Pauli Virtanen]

Hi,

la, 2008-03-01 kello 16:43 -0300, Lisandro Dalcin kirjoitti:
> I want to comment some extrange stuff I'm experiencing with numpy.
> Please, let me know if this is expected and known.
> 
> I'm trying to solve a model nonlinear PDE, 2D Bratu problem (-Lapacian
> u - alpha * exp(u), homogeneus bondary conditions), using the simple
> finite differences with a 5-point stencil.
> 
> I implemented the finite diference scheme in pure-numpy, and also in a
> F90 subroutine, next wrapped with f2py.
> 
> Next, I use PETSc (through petsc4py) to solve the problem with a
> Newton method, a Krylov solver, and a matrix-free technique for the
> Jacobian (that is, the matrix is never explicitelly assembled, its
> action on a vector is approximated again with a 1st. order finite
> direrence formula).
> 
> And the, surprise! The pure-numpy implementation accumulates many more
> inner linear iterations (about 25%) in the complete nonlinear solution
> loop than the one using the F90 code wrapped with f2py.

A silly question: did you check directly that the pure-numpy code and
the F90 code give the same results for the Jacobian-vector product 
J(z0) z for some randomly chosen vectors z0, z?

-- 
Pauli Virtanen