VB.NET Hermitian Band Matrix Example

Imports System

Imports CenterSpace.NMath.Core
Imports CenterSpace.NMath.Matrix

Namespace CenterSpace.NMath.Matrix.Examples.VisualBasic

  ' A .NET example in VB.NET demonstrating the features of the banded Hermitian matrix classes.
  Module HermitianBandMatrixExample

    Sub Main()

      ' Set up the parameters that describe the shape of a Hermitian banded matrix.
      Dim HalfBandwidth As Integer = 1
      Dim Order As Integer = 7

      ' Set up an Hermitian banded matrix B by creating a banded complex matrix and setting
      ' B equal to the product of the conjugate transpose of that matrix with itself.
      Dim A As New FloatComplexBandMatrix(Order, Order, HalfBandwidth, HalfBandwidth)
      Dim I As Integer
      For I = -A.LowerBandwidth To A.UpperBandwidth
        A.Diagonal(I).Set(Slice.All, New FloatComplex(I + 2, I - 1))
      Next
      Dim B As New FloatHermitianBandMatrix(MatrixFunctions.Product(A.ConjTranspose(), A))
      Console.WriteLine()
      Console.Write("B = ")
      Console.WriteLine(B.ToString())
      Console.WriteLine()

      ' B = 2 7x7 [ (10,0)  (10,6)  (3,6)   (0,0)   (0,0)   (0,0)   (0,0)  
      '             (10,-6) (19,0)  (10,6)  (3,6)   (0,0)   (0,0)   (0,0)  
      '             (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)   (0,0)   (0,0)  
      '             (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)   (0,0)  
      '             (0,0)   (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)  
      '             (0,0)   (0,0)   (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)  
      '             (0,0)   (0,0)   (0,0)   (0,0)   (3,-6)  (10,-6) (14,0) ]

      ' Indexer accessor works just like it does for general matrices. 
      Console.Write("B[2,2] = ")
      Console.WriteLine(B(2, 2))
      Console.Write("B[5,0] = ")
      Console.WriteLine(B(5, 0))

      ' You can set the values of elements in the bandwidth
      ' of an Hermitian banded matrix using the indexer. Note that setting 
      ' the element in row i and column j to a value implicitly sets the 
      ' element in column j and row i to the complex conjugate of that value.
      Dim Z As New FloatComplex(99, -99)
      B(2, 1) = Z
      Console.Write("B[2,1] = ")
      Console.WriteLine(B(2, 1)) ' (99, -99)
      Console.Write("B[1,2] = ")
      Console.WriteLine(B(1, 2)) ' (99, 99 )

      ' But setting an element outside the bandwidth of the
      ' matrix raises a NonModifiableElementException exception.
      Try
        B(6, 0) = Z
      Catch E As NonModifiableElementException
        Console.WriteLine()
        Console.Write("NonModifiableElementException: " & E.Message)
      End Try

      ' Scalar multiplication and matrix addition/subtraction are supported.
      Z = New FloatComplex(0.0001, 0.0001)
      Dim C As FloatHermitianBandMatrix = Z * B
      Dim D As FloatHermitianBandMatrix = C - B
      Console.WriteLine()
      Console.Write("D = ")
      Console.WriteLine(D.ToString())

      ' Matrix/vector inner products too.
      Dim Rng As New RandGenUniform(-1, 1)
      Rng.Reset(&H124)
      Dim X As New FloatComplexVector(B.Cols, Rng)
      Dim Y As FloatComplexVector = MatrixFunctions.Product(B, X)
      Console.WriteLine()
      Console.Write("Bx = ")
      Console.WriteLine(Y.ToString())

      ' You can transform the non-zero elements of a banded matrix object by using
      ' the Transform() method on its data vector.
      C.DataVector.Transform(NMathFunctions.FloatComplexExpFunction)
      Console.WriteLine()
      Console.Write("exp(C) = ")
      Console.WriteLine(C.ToString())

      ' You can also solve linear systems.
      Dim X2 As FloatComplexVector = MatrixFunctions.Solve(B, Y)

      ' x and x2 should be the same. Let's look at the l2 norm of
      ' their difference.
      Dim Residual As FloatComplexVector = FloatComplexVector.Subtract(X, X2)
      Dim ResidualL2Norm As Double = Math.Sqrt(NMathFunctions.ConjDot(Residual, Residual).Real)
      Console.WriteLine()
      Console.Write("||x - x2|| = ")
      Console.WriteLine(ResidualL2Norm)

      ' You can calculate inverses too.
      Dim BInv As FloatComplexMatrix = MatrixFunctions.Inverse(B)
      Console.WriteLine()
      Console.Write("BInv = ")
      Console.WriteLine(BInv.ToString())

      ' If your Hermitian banded matrix is positive definite, you can invoke
      ' the Solve function with a third - the isPositiveDefinite 
      ' parameter - set to true. 
      ' First make sure B is positive definite.
      B = New FloatHermitianBandMatrix(MatrixFunctions.Product(A.ConjTranspose(), A))
      Y = MatrixFunctions.Product(B, X)
      X2 = MatrixFunctions.Solve(B, Y, True)   ' 3rd parameter isPositiveDefinite set to true.

      ' See how close Bx is to y by computing the l2 norm of their difference.
      Residual = X - X2
      ResidualL2Norm = Math.Sqrt(NMathFunctions.ConjDot(Residual, Residual).Real)
      Console.WriteLine()
      Console.Write("PD ||x - x2|| = ")
      Console.WriteLine(ResidualL2Norm)

      ' You can use the Resize() method to change the bandwidths.
      D.Resize(D.Order, 2)

      Console.WriteLine()
      Console.WriteLine("Press Enter Key")
      Console.Read()

    End Sub
  End Module
End Namespace