Iterative Sparse Solvers in IronPython QuickStart Sample

Illustrates the use of iterative sparse solvers and preconditioners for efficiently solving large, sparse systems of linear equations in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to solve large sparse systems of linear equations using iterative solvers and preconditioners in Numerics.NET.

The sample shows two main examples:

  1. Solving a non-symmetric sparse system loaded from Matrix Market files using the BiConjugate Gradient solver. It demonstrates how to:
    • Load sparse matrices from files
    • Create and configure an iterative solver
    • Apply preconditioners to improve convergence
    • Monitor solution progress and error estimates
  2. Solving the Laplace equation on a rectangular grid using a symmetric solver. This example shows how to:
    • Construct a sparse matrix representing the discretized Laplace operator
    • Set up boundary conditions
    • Solve the resulting system using the Quasi-Minimal Residual method
    • Handle symmetric vs non-symmetric systems appropriately

The sample illustrates proper error handling, solver selection based on matrix properties, and the benefits of preconditioning for convergence acceleration. It provides a practical introduction to solving large-scale sparse linear systems efficiently.

The code

import numerics

from Extreme.Mathematics import *
# Sparse matrices are in the Extreme.Mathematics.LinearAlgebra
# namespace
from Extreme.Mathematics.LinearAlgebra import *
from Extreme.Mathematics.LinearAlgebra.IterativeSolvers import *
from Extreme.Mathematics.LinearAlgebra.IterativeSolvers.Preconditioners import *
from Extreme.Mathematics.LinearAlgebra.IO import *

# Illustrates the use of iterative sparse solvers for efficiently
# solving large, sparse systems of linear equations using the 
# iterative sparse solver and preconditioner classes from the
# Extreme.Mathematics.LinearAlgebra.IterativeSolvers namespace of 
# Extreme Numerics.NET.

# This QuickStart Sample illustrates how to solve sparse linear systems
# using iterative solvers.

# IterativeSparseSolver is the base class for all
# iterative solver classes:

#
# Non-symmetric systems
#

print "Non-symmetric systems"

# We load a sparse matrix and right-hand side from a data file:
A = MatrixMarketReader.ReadMatrix(r"..\data\sherman3.mtx")
b = MatrixMarketReader.ReadMatrix(r"..\data\sherman3_rhs1.mtx").GetColumn(0)

print "Solve Ax = b"
print "A is {0}x{1} with {2} nonzeros.".format(A.RowCount, A.ColumnCount, A.NonzeroCount)

# Some solvers are suitable for symmetric matrices only.
# Our matrix is not symmetric, so we need a solver that
# can handle this:
solver = BiConjugateGradientSolver(A)

solver.Solve(b)
print "Solved in {0} iterations.".format(solver.IterationsNeeded)
print "Estimated error:", solver.SolutionReport.Error

# Using a preconditioner can improve convergence. You can use
# one of the predefined preconditioners, or supply your own.

# With incomplete LU preconditioner
solver.Preconditioner = IncompleteLUPreconditioner(A)
solver.Solve(b)
print "Solved in {0} iterations.".format(solver.IterationsNeeded)
print "Estimated error:", solver.EstimatedError

# 
# Symmetrical systems
# 

print "Symmetric systems"

# In this example we solve the Laplace equation on a rectangular grid
# with Dirichlet boundary conditions.

# We create 100 divisions in each direction, giving us 99 interior points
# in each direction:
nx = 99
ny = 99

# The boundary conditions are just some arbitrary functions.
left = Vector.Create(ny, lambda i: (i / (nx - 1.0)) ** 2)
right = Vector.Create(ny, lambda i: 1 - (i / (nx - 1.0)))
top = Vector.Create(nx, lambda i: Elementary.SinPi(5 * (i / (nx - 1.0))) )
bottom = Vector.Create(nx, lambda i: Elementary.CosPi(5 * (i / (nx - 1.0))) )

# We discretize the Laplace operator using the 5 point stencil.
laplacian = Matrix.CreateSparse(nx * ny, nx * ny, 5 * nx * ny)
rhs = Vector.Create(nx * ny)
for j in range(0,ny):
    for i in range(0,nx):
        ix = j * nx + i
        if (j > 0):
            laplacian[ix, ix - nx] = 0.25
        if (i > 0):
            laplacian[ix, ix - 1] = 0.25
        laplacian[ix, ix] = -1.0
        if (i + 1 < nx):
            laplacian[ix, ix + 1] = 0.25
        if (j + 1 < ny):
            laplacian[ix, ix + nx] = 0.25

# We build up the right-hand sides using the boundary conditions:
for i in range(0,nx):
    rhs[i] = -0.25 * top[i]
    rhs[nx * (ny - 1) + i] = -0.25 * bottom[i]
for j in range(0,ny):
    rhs[j * nx] -= 0.25 * left[j]
    rhs[j * nx + nx - 1] -= 0.25 * right[j]

# Finally, we create an iterative solver suitable for
# symmetric systems...
solver = QuasiMinimalResidualSolver(laplacian)
# and solve using the right-hand side we just calculated:
solver.Solve(rhs)

print "Solve Ax = b"
print "A is {0}x{1} with {2} nonzeros.".format(A.RowCount, A.ColumnCount, A.NonzeroCount)
print "Solved in {0} iterations.".format(solver.IterationsNeeded)
print "Estimated error:", solver.EstimatedError