Iterative Sparse Solvers QuickStart Sample (Visual Basic)

Option Infer On

' The DenseMatrix and DoubleVector classes resides in the 
' Extreme.Mathematics.LinearAlgebra namespace.
Imports Extreme.Mathematics
Imports Extreme.Mathematics.LinearAlgebra
Imports Extreme.Mathematics.LinearAlgebra.IterativeSolvers
Imports Extreme.Mathematics.LinearAlgebra.IterativeSolvers.Preconditioners
Imports Extreme.Mathematics.LinearAlgebra.IO

Namespace Extreme.Numerics.QuickStart.VB
    ' 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.
    Module IterativeSparseSolvers

        Sub Main()
            ' This QuickStart Sample illustrates how to solve sparse linear systems
            ' using iterative solvers.

            ' IterativeSparseSolver is the base class for all
            ' iterative solver classes:
            Dim solver As IterativeSparseSolver(Of Double)

            '
            ' Non-symmetric systems
            '

            Console.WriteLine("Non-symmetric systems")

            ' We load a sparse matrix and right-hand side from a data file:
            Dim A = CType(MatrixMarketReader.ReadMatrix("..\..\..\..\data\sherman3.mtx"),
                SparseCompressedColumnMatrix(Of Double))
            Dim b = MatrixMarketReader.ReadMatrix("..\..\..\..\data\sherman3_rhs1.mtx").GetColumn(0)

            Console.WriteLine("Solve Ax = b")
            Console.WriteLine("A is {0}x{1} with {2} nonzeros.", 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:

            ' #Region DOCIterativeSparseSolvers1
            solver = New BiConjugateGradientSolver(Of Double)(A)

            solver.Solve(b)
            Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded)
            Console.WriteLine("Estimated error: {0}", solver.SolutionReport.Error)
            ' #End Region

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

            ' With incomplete LU preconditioner
            solver.Preconditioner = New IncompleteLUPreconditioner(Of Double)(A)
            solver.Solve(b)
            Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded)
            Console.WriteLine("Estimated error: {0}", solver.EstimatedError)

            ' 
            ' Symmetrical systems
            ' 

            Console.WriteLine("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:
            Const nx = 99
            Const ny = 99

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

            ' We discretize the Laplace operator using the 5 point stencil.
            Dim laplacian = Matrix.CreateSparse(Of Double)(nx * ny, nx * ny, 5 * nx * ny)
            Dim rhs = Vector.Create(Of Double)(nx * ny)
            For j As Integer = 0 To ny - 1
                For i As Integer = 0 To nx - 1
                    Dim ix As Integer = j * nx + i
                    If (j > 0) Then laplacian(ix, ix - nx) = 0.25
                    If (i > 0) Then laplacian(ix, ix - 1) = 0.25
                    laplacian(ix, ix) = -1.0
                    If (i + 1 < nx) Then laplacian(ix, ix + 1) = 0.25
                    If (j + 1 < ny) Then laplacian(ix, ix + nx) = 0.25
                Next
            Next
            ' We build up the right-hand sides using the boundary conditions:
            For i As Integer = 0 To nx - 1
                rhs(i) = -0.25 * top(i)
                rhs(nx * (ny - 1) + i) = -0.25 * bottom(i)
            Next
            For j As Integer = 0 To ny - 1
                rhs(j * nx) -= 0.25 * left(j)
                rhs(j * nx + nx - 1) -= 0.25 * right(j)
            Next

            Console.WriteLine("A is {0}x{1} with {2} nonzeros.", laplacian.RowCount, laplacian.ColumnCount, laplacian.NonzeroCount)

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

            Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded)
            Console.WriteLine("Estimated error: {0}", solver.EstimatedError)


            Console.Write("Press Enter key to exit...")
            Console.ReadLine()
        End Sub

    End Module

End Namespace