[SciPy-user] fsolve with sparse matrices

[Gergely Imreh]

2009/4/17 Jan Wicijowski <janislaw at o2.pl>:
> Hi all!
>
> I would like to ask you, if anyone was ever confronted with solving
> nonlinear system with scipy.optimize.fsolve, with huge no. of equations,
> say 10000.
>
> The usual approach fails, because the size of jacobian matrix is
> proportional to the square of equations number. The paging kills my OS.
>
> In my problem, which is solving a nonlinear electrical circuit, the matrix
> is mostly zeros, which has immediately led me to trying sparse matrices.
>
> The first obstacle was numpy.atleast_1d contained in
> scipy.optimize.minpack.check_func - instead of preserving sparse matrix as
> is, it has a suprising behaviour of boxing the matrix in a list. This may
> be a bug, but I may post it elsewhere.
>
> Injecting custom function instead of scipy.optimize.minpack.check_func
> didn't work as well with fsolve. This wasn't surprising, as I guessed,
> that FORTRAN hybrj won't be able to deal with interfacing with scipy
> matrices.
>
> So, am I left with rewriting the original FORTRAN hybrj source to python,
> or is there somebody, who dealt with such problem?
>
> Regards,
> Jan Wicijowski


Hi,

  I was working on a similar problem (large sparse matrix, solving
optical Bloch-equations), and been planning to check out interfacing
with Fortran. The idea was, that I know the NAG library [1] works very
well for these kinds of calculations, just have to get the interface
between the sparse matrices and Fortran working. I checked it now, and
they have some documents showing how to use the NAG with F2PY [2].
Though they only use numpy, so no mention of sparse matrices, but NAG
does sparse.
  This supposed to be just a heads up, need to check it out. Though in
my experience even the NAG people can be contacted, and they helped to
resolve issues (or extend functionality) before.

  Cheers,
     Greg

[1] http://www.nag.co.uk/
[2] http://www.nag.co.uk/doc/techrep/index.asp#np3665