New Version 6.0!

Try it for free with our fully functional 60-day trial version.

Download now!

QuickStart Samples

Iterative Sparse Solvers QuickStart Sample (Visual Basic)

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

C# code F# code IronPython code Back to QuickStart Samples

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