[Scipy-tickets] [SciPy] #360: Callback function is broken in linalg.gmres
SciPy
scipy-tickets@scipy....
Fri May 4 08:16:55 CDT 2007
#360: Callback function is broken in linalg.gmres
--------------------------+-------------------------------------------------
Reporter: nils | Owner: somebody
Type: defect | Status: new
Priority: normal | Milestone: 0.5.3 Release
Component: scipy.linalg | Version: devel
Severity: normal | Resolution:
Keywords: |
--------------------------+-------------------------------------------------
Comment (by bart):
The problem is that the wrapped GMRES routine does not need to build the
intermediate iterates x to compute the residual. See
http://www.netlib.org/linalg/html_templates/node29.html
For performance this is a good thing, but for tracking convergence not ...
What you can do, is to give the callback routine the residual (actually
the relative residual). The iter_ variable is not incremented during each
iteration but the code does return to the python wrapper. If you follow
the logic of a GMRES(m) algorithm, you can extract the new residuals by
checking the ijob variabele
matvec (ijob 1)
loop for j=1..maxit (this increases iter_)
psolve (ijob 2)
loop for i=1..m
matvec (ijob 3)
psolve (ijob 2)
new residual available
if res<tol: stop with ijob -1
update x (this is a new iterate x)
matvec (ijob 1)
So we have a new residual when
- ijob 3 after ijob 2, but not for the first time
- ijob -1
This is implemented in the attached new_iterative_gmres.py I'm pretty sure
that normal GMRES gives all the residuals, for restarted GMRES I'm not
sure...
--
Ticket URL: <http://projects.scipy.org/scipy/scipy/ticket/360#comment:1>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list