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.

View this sample in: C# Visual Basic F#

```Python
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 
# 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

```