C# Herm PD Tri Diag Fact Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the factorization classes for
  /// Hermitian positive definite tridiagonal matrices.
  /// </summary>
  class HermPDTriDiagFactExample
  {

    static void Main(string[] args)
    {
      // Construct a positive definite tridiagonal matrix.
      int rows = 5;
      int cols = 5;
      DoubleComplexVector data1 = new DoubleComplexVector(cols, 2);
      DoubleComplexVector data2 = new DoubleComplexVector(cols - 1, -1);
      DoubleComplexTriDiagMatrix A = new DoubleComplexTriDiagMatrix(rows, cols);
      A.Diagonal()[Slice.All] = data1;
      A.Diagonal(1)[Slice.All] = data2;
      A.Diagonal(-1)[Slice.All] = data2;

      Console.WriteLine();
      Console.WriteLine("A = {0}", A.ToString());
      Console.WriteLine();

      // 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.
      DoubleHermPDTriDiagFact fact = new DoubleHermPDTriDiagFact(A);

      // Check to see if A is positive definite.
      string isPDString = fact.IsPositiveDefinite ? "A is positive definite" : "A is NOT positive definite";
      Console.WriteLine(isPDString);

      // Retrieve information about the matrix A.
      DoubleComplex det = fact.Determinant();
      double rcond = fact.ConditionNumber();
      DoubleComplexMatrix AInv = fact.Inverse();

      Console.WriteLine();
      Console.WriteLine("Determinant of A = {0}", det);

      Console.WriteLine();
      Console.WriteLine("Reciprocal condition number = {0}", rcond);

      Console.WriteLine();
      Console.WriteLine("A inverse = {0}", AInv.ToString());

      // Use the factorization to solve some linear systems Ax = y.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      DoubleComplexVector y0 = new DoubleComplexVector(fact.Cols, rng);
      DoubleComplexVector y1 = new DoubleComplexVector(fact.Cols, rng);
      DoubleComplexVector x0 = fact.Solve(y0);
      DoubleComplexVector x1 = fact.Solve(y1);

      Console.WriteLine();
      Console.WriteLine("Solution to Ax = y0 is {0}", x0.ToString());

      Console.WriteLine();
      Console.WriteLine("y0 - Ax0 = {0}", (y0 - MatrixFunctions.Product(A, x0)).ToString());

      Console.WriteLine();
      Console.WriteLine("Solution to Ax = y1 is {0}", x1.ToString());

      Console.WriteLine();
      Console.WriteLine("y1 - Ax1 = {0}", (y1 - MatrixFunctions.Product(A, x1)).ToString());

      // You can also solve for multiple right-hand sides.
      DoubleComplexMatrix Y = new DoubleComplexMatrix(y1.Length, 2);
      Y.Col(0)[Slice.All] = y0;
      Y.Col(1)[Slice.All] = y1;
      DoubleComplexMatrix X = fact.Solve(Y);

      // The first column of X should be x0; the second column should be x1.
      Console.WriteLine();
      Console.WriteLine("X = {0}", X.ToString());

      // Factor a different matrix.
      DoubleComplexTriDiagMatrix B = A + A;
      fact.Factor(B);
      x0 = fact.Solve(y0);

      Console.WriteLine();
      Console.WriteLine("Solution to Bx = y0 is {0}", x0.ToString());

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