3. Generic Template for Krylov Methods¶

3.1. The generic Module¶

class generic.KrylovMethod(op, **kwargs)¶

Bases: object

A general template for implementing iterative Krylov methods. This module defines the KrylovMethod generic class. Other modules subclass KrylovMethod to implement specific algorithms.

Parameters:	op: an operator describing the coefficient matrix A. y = op * x must return the operator-vector product y = Ax for any given vector x. If required by the method, y = op.T * x must return the operator-vector product with the adjoint operator.
Keywords:	atol: absolute stopping tolerance. Default: 1.0e-8. rtol: relative stopping tolerance. Default: 1.0e-6. precon: optional preconditioner. If not None, y = precon(x) returns the vector y solution of the linear system M y = x. logger: a logging.logger instance. If none is supplied, a default null logger will be used.

For general references on Krylov methods, see [Demmel], [Greenbaum], [Kelley], [Saad] and [Templates].

References:

[Demmel]

J. W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.

[Greenbaum]

A. Greenbaum, Iterative Methods for Solving Linear Systems, number 17 in Frontiers in Applied Mathematics, SIAM, Philadelphia, 1997.

[Kelley]

C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, number 16 in Frontiers in Applied Mathematics, SIAM, Philadelphia, 1995.

[Saad]

Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, Philadelphia, 2003.

[Templates]

R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1993.

solve(rhs, **kwargs)¶

This is the solve() method of the abstract KrylovMethod class. The class must be specialized and this method overridden.

3.2. Writing a New Solver¶

Adding a new solver to PyKrylov should be done by subclassing KrylovMethod from the generic module and overriding the solve() method. A general template might look like the following:

import numpy as np
from pykrylov.generic import KrylovMethod

class NewSolver( KrylovMethod ):
    """
    Document your new class and give adequate references.
    """

    def __init__(self, matvec, **kwargs):
        KrylovMethod.__init__(self, matvec, **kwargs)

        self.name = 'New Krylov Solver'
        self.acronym = 'NKS'
        self.prefix = self.acronym + ': '  # Used when verbose=True


    def solve(self, rhs, **kwargs):
        """
        Solve a linear system with `rhs` as right-hand side by the NKS
        method. The vector `rhs` should be a Numpy array. An optional
        argument `guess` may be supplied, with an initial guess as a Numpy
        array. By default, the initial guess is the vector of zeros.
        """
        n = rhs.shape[0]
        nMatvec = 0

        # Initial guess is zero unless one is supplied
        guess_supplied = 'guess' in kwargs.keys()
        x = kwargs.get('guess', np.zeros(n))

        r0 = rhs  # Fixed vector throughout
        if guess_supplied:
            r0 = rhs - self.matvec(x)

        # Further initializations ...

        # Compute initial residual norm. For example:
        residNorm = np.linalg.norm(r0)
        self.residNorm0 = residNorm

        # Compute stopping threshold
        threshold = max( self.abstol, self.reltol * self.residNorm0 )

        if self.verbose:
            self._write('Initial residual = %8.2e\n' % self.residNorm0)
            self._write('Threshold = %8.2e\n' % threshold)

        # Main loop
        while residNorm > threshold and nMatvec < self.matvec_max:

            # Do something ...
            pass

        # End of main loop

        self.nMatvec = nMatvec
        self.bestSolution = x
        self.residNorm = residNorm