[Numpy-discussion] numpy and roundoff(?)
Sat Mar 1 13:43:56 CST 2008
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
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 +
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
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?
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
More information about the Numpy-discussion