[SciPy-User] gmres callback with restart unworking ?

2 Dec 2015

      Dear,

I would like to catch the full convergence history for all the iterations,
 - no restart: gmres(full)
 - with restart: gmres(restart)
and in the same time, be able to fix the maximum number of total
iterations, to avoid surprise.
I would like to track the residual of all the inner iterations.

What is the more efficient solution ? because I have remarked a
twisted behavior with 'callback'.
I guess that the test in 'ijob==3' is not fine enough.

'If callback' is None,
then the number of total iterations corresponds to 'maxiter*restart',
and 'maxiter' fixes the number of outer iterations,
'restart' fixes the number of inner iterations.
Note that if 'restart' is None, the default parameter is 20.

This behavior is the expected one, right ?

If 'callback' is not None,
then the number of total iterations corresponds to 'maxiter',
but the 'restart' parameter is never considered.

In both case, with or without 'callback',
if 'maxiter' is None, then GMRes is running more the size of the
matrix which does not make so much sense from a linear algebra point
of view.

Note that the number of matvec performed corresponds to the number of
total iterations plus one.
Which is expected due to the initial residual.
However, the 'callback' does not catch it. Why not add it ?

The preconditioner, is it applied at Left or at Right ?
What is the best strategy to apply a preconditioner depending on the
iteration number ?

Sould be a good idea or not to maintain a pythonic implementation of
gmres similar to lgmres ?

Last, but not least, where could I find the doc about the 'ijob'
numbering and the Fortran API ?

Maybe I am missing a point ?

All the best,
simon

~~python

#!/usr/env python
# coding: utf8

import numpy
import scipy.linalg
import scipy.sparse.linalg

def f(x):
    global countMatvec
    countMatvec += 1
    y = numpy.random.rand(x.shape[0])
    return y

def catch(res):
    global countRes, reslist
    countRes += 1
    reslist.append(res)

N = 50

A = scipy.sparse.linalg.LinearOperator((N, N), matvec=f, dtype=complex)
b = numpy.random.rand(A.shape[0])

m = 3
r = 7

if m is None: mm = 1
else: mm = m
if r is None: rr = 20
else: rr = r

print("\nno callback")
countMatvec = 0
countRes, reslist = 0, []
y, info = scipy.sparse.linalg.gmres(A, b,
                                    maxiter=m, restart=r)
print(countMatvec-1, countRes, mm*(rr+1)-1, info==mm, info)
print(numpy.linalg.norm(b - A.matvec(y)))

print("\ncallback")
countMatvec = 0
countRes, reslist = 0, []
y, info = scipy.sparse.linalg.gmres(A, b,
                                    maxiter=m, restart=r, callback=catch)

print(countMatvec-1, countRes, mm*(rr+1)-1, info==mm, info)
print(numpy.linalg.norm(b - A.matvec(y)))

[SciPy-User] gmres callback with restart unworking ?

zimon toutoune