[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