Iterative Sparse Solvers QuickStart Sample (VB.NET)

Illustrates the use of iterative sparse solvers for efficiently solving large, sparse systems of linear equations using the iterative sparse solver and preconditioner classes in Visual Basic .NET.

C# code F# code Back to QuickStart Samples

' 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

            '
            ' Non-symmetric systems
            '

            Console.WriteLine("Non-symmetric systems")

            ' We load a sparse matrix and right-hand side from a data file:
            Dim A As SparseCompressedColumnMatrix = CType(MatrixMarketReader.ReadMatrix("..\..\..\..\data\sherman3.mtx"),  _
                SparseCompressedColumnMatrix)
            Dim b As Vector = 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:
            solver = New BiConjugateGradientSolver(A)

            ' With default (Jacobi) preconditioner:
            solver.Solve(b)
            Console.WriteLine("Solved in {0} iterations.", solver.IterationsNeeded)
            Console.WriteLine("Estimated error: {0}", 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 = New IncompleteLUPreconditioner(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 As Vector = Vector.Create(ny, Function(i) (i / (nx + 1.0)) ^ 2)
            Dim right As Vector = Vector.Create(ny, Function(i) 1 - (i / (nx + 1.0)))
            Dim top As Vector = Vector.Create(nx, Function(i) Elementary.SinPi(5 * (i / (nx + 1.0))))
            Dim bottom As Vector = Vector.Create(nx, Function(i) Elementary.CosPi(5 * (i / (nx + 1.0))))

            ' We discretize the Laplace operator using the 5 point stencil.
            Dim laplacian As SparseCompressedColumnMatrix = Matrix.CreateSparse(nx * ny, nx * ny, 5 * nx * ny)
            Dim rhs As Vector = Vector.Create(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(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