1. Introduction to PyKrylov¶

PyKrylov aims to provide a flexible implementation of the most common Krylov method for solving systems of linear equations in pure Python. The sole requirement is Numpy for fast array operations.

PyKyrlov is in the very early stages of development. If you give it a whirl, please let me know what you think at the LightHouse PyKrylov page. Feel free to post bug reports, feature requests and praises.

1.1. Example¶

The following code snippet solves a linear system with a 1D Poisson coefficient matrix. Since this matrix is symmetric and positive definite, the conjugate gradient algorithm is used:

import numpy as np
from math import sqrt
from pykrylov.cg import CG

def Poisson1dMatvec(x):
    # Matrix-vector product with a 1D Poisson matrix
    y = 2*x
    y[:-1] -= x[1:]
    y[1:] -= x[:-1]
    return y

n = 100
e = np.ones(n)
rhs = Poisson1dMatvec(e)
cg = CG(Poisson1dMatvec, matvec_max=200)
cg.solve(rhs)

print 'Number of matrix-vector products: ', cg.nMatvec
print 'Residual: %8.2e' % cg.residNorm
print 'Error: %8.2e' % (np.linalg.norm(e - CG.bestSolution)/sqrt(n))

On my machine, the above script produces:

Number of matrix-vector products: 50
Residual: 7.39e-14
Error: 2.06e-15