[Numpy-discussion] numpy and roundoff(?)

Lisandro Dalcin dalcinl@gmail....
Sat Mar 1 15:08:18 CST 2008

```Dear Charles,

As I said, I have no time to code the pure Python+numpy nonlinear and
linear loops, and the matrix-free stuff to mimic the PETSc
implementation. However, I post the F90 code and the numpy code, and a
small script for testing with random input. When I have some spare
time, I'll try to do the complete application in pure python.

Regards,

On 3/1/08, Charles R Harris <charlesr.harris@gmail.com> wrote:
>
>
>
> On Sat, Mar 1, 2008 at 12:43 PM, Lisandro Dalcin <dalcinl@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
>
>
>
> _______________________________________________
>  Numpy-discussion mailing list
>  Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
>

--
Lisandro Dalcín
---------------
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
Tel/Fax: +54-(0)342-451.1594
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bratu2dlib.f90
Type: application/octet-stream
Size: 861 bytes
Desc: not available
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bratu2dnpy.py
Type: text/x-python
Size: 454 bytes
Desc: not available