[Numpy-discussion] numpy and roundoff(?)

[Charles R Harris]

On Sat, Mar 1, 2008 at 12:43 PM, Lisandro Dalcin <dalcinl at gmail.com> wrote:

> Dear all,
>
> 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.
>
> Additionally, PETSc have in its source distribution a similar example,
> but implemented in C and using some PETSc utilities for managing
> structured grids. In short, this code is in C and completelly
> unrelated to the previously commented code. After running this
> example, I get almost the same results that the one for my petsc4py +
> F90 code.
>
> All this surprised me. It seems that for some reason numpy is
> accumulating some roundoff, and this is afecting the acuracy of the
> aproximated Jacobian, and then the linear solvers need more iteration
> to converge.
>
> Unfortunatelly, I cannot offer a self contained example, as this code
> depends on having PETSc and petsc4py. Of course, I could write myself
> the nonlinear loop, and a CG solver, but I am really busy.
>
> Can someone comment on this? Is all this expected? Have any of you
> experienced somethig similar?
>
> Could you attach the pure numpy solution along with a test case (alpha=?).

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080301/d763bb51/attachment.html>