Iterative Sparse Solvers QuickStart Sample (IronPython)

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 
#/ the Extreme Optimization Numerical Libraries for .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