VB.NET Herm PD Tri Diag Fact 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 factorization classes for
  ' Hermitian positive definite tridiagonal matrices.
  Module HermPDTriDiagFactExample

    Sub Main()

      ' Construct a positive definite tridiagonal matrix.
      Dim Rows As Integer = 5
      Dim Cols As Integer = 5
      Dim Two As New DoubleComplex(2, 0)
      Dim NegOne As New DoubleComplex(-1, 0)
      Dim Data1 As New DoubleComplexVector(Cols, Two)
      Dim Data2 As New DoubleComplexVector(Cols - 1, NegOne)
      Dim A As New DoubleComplexTriDiagMatrix(Rows, Cols)
      A.Diagonal()(Slice.All) = Data1
      A.Diagonal(1)(Slice.All) = Data2
      A.Diagonal(-1)(Slice.All) = Data2

      Console.WriteLine()
      Console.Write("A = ")
      Console.WriteLine(A.ToString())

      ' A = 5x5 [ (2,0)  (-1,0) (0,0)  (0,0)  (0,0)
      '           (-1,0) (2,0)  (-1,0) (0,0)  (0,0)
      '           (0,0)  (-1,0) (2,0)  (-1,0) (0,0)
      '           (0,0)  (0,0)  (-1,0) (2,0)  (-1,0)
      '           (0,0)  (0,0)  (0,0)  (-1,0) (2,0) ]

      ' Construct a positive definite tridiagonal factorization class.
      Dim Fact As New DoubleHermPDTriDiagFact(A)

      ' Check to see if A is positive definite.
      Dim IsPDString As String
      If Fact.IsPositiveDefinite Then
        IsPDString = "A is positive definite"
      Else
        IsPDString = "A is NOT positive definite"
      End If
      Console.WriteLine()
      Console.WriteLine(IsPDString)

      ' Retrieve information about the matrix A.
      Dim Det As DoubleComplex = Fact.Determinant()
      Dim RCond As Double = Fact.ConditionNumber()
      Dim AInv As DoubleComplexMatrix = Fact.Inverse()

      Console.WriteLine()
      Console.Write("Determinant of A = ")
      Console.WriteLine(Det)

      Console.WriteLine()
      Console.Write("Reciprocal condition number = ")
      Console.WriteLine(RCond)

      Console.WriteLine()
      Console.Write("A inverse = ")
      Console.WriteLine(AInv.ToString())

      ' Use the factorization to solve some linear systems Ax = y.
      Dim Rng As New RandGenUniform(-1, 1)
      Rng.Reset(&H124)
      Dim Y0 As New DoubleComplexVector(Fact.Cols, Rng)
      Dim Y1 As New DoubleComplexVector(Fact.Cols, Rng)
      Dim X0 As DoubleComplexVector = Fact.Solve(Y0)

      Dim X1 As DoubleComplexVector = Fact.Solve(Y1)

      Console.WriteLine()
      Console.Write("Solution to Ax = y0 is ")
      Console.WriteLine(X0.ToString())

      Console.WriteLine()
      Console.Write("y0 - Ax0 = ")
      Console.WriteLine((Y0 - MatrixFunctions.Product(A, X0)).ToString())

      Console.WriteLine()
      Console.Write("Solution to Ax = y1 is ")
      Console.WriteLine(X1.ToString())

      Console.WriteLine()
      Console.Write("y1 - Ax1 = ")
      Console.WriteLine((Y1 - MatrixFunctions.Product(A, X1)).ToString())

      ' You can also solve for multiple right-hand sides.
      Dim Y As New DoubleComplexMatrix(Y1.Length, 2)
      Y.Col(0)(Slice.All) = Y0
      Y.Col(1)(Slice.All) = Y1
      Dim X As DoubleComplexMatrix = Fact.Solve(Y)

      ' The first column of X should be x0 the second column should be x1.
      Console.WriteLine()
      Console.Write("X = ")
      Console.WriteLine(X.ToString())

      ' Factor a different matrix.
      Dim B As DoubleComplexTriDiagMatrix = A + A
      Fact.Factor(B)
      X0 = Fact.Solve(Y0)
      Console.WriteLine()
      Console.Write("Solution to Bx = y0 is ")
      Console.WriteLine(X0.ToString())

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

    End Sub
  End Module
End Namespace