C# Band 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
  /// banded matrices.
  /// </summary>
  class BandFactExample
  {

    static void Main(string[] args)
    {
      // Construct a banded matrix with random entries from a data vector
      // of the appropriate length,
      int upperBandwidth = 2;
      int lowerBandwidth = 1;
      int rows = 7;
      int cols = 7;
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      FloatVector data = new FloatVector((lowerBandwidth + upperBandwidth + 1) * cols, rng);
      FloatBandMatrix A = new FloatBandMatrix(data, rows, cols, lowerBandwidth, upperBandwidth);

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

      // A = 1 2 7x7 [ -0.250345   0.1820724  -0.3249055  0           0           0          0
      //                0.5763928 -0.04433365 -0.2796601  0.2204346   0           0          0 
      //                0          0.07382878  0.5600654 -0.07704575 -0.9237711   0          0  
      //                0          0           0.306021  -0.1677762   0.6205534  -0.8629085  0 
      //                0          0           0         -0.2589606   0.7733703   0.2029475 -0.7045102  
      //                0          0           0          0           0.3579623  -0.5847561  0.1235993 
      //                0          0           0          0           0           0.6222407  0.870272 ]

      // Construct a band factorization class.
      FloatBandFact fact = new FloatBandFact(A);
      // Check to see if A is singular:
      string isSingularString = fact.IsSingular ? "A is singular" : "A is NOT singular";
      Console.WriteLine();
      Console.WriteLine(isSingularString);

      // Retrieve information about the matrix A:
      float det = fact.Determinant();
      float rcond = fact.ConditionNumber();
      FloatMatrix 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.
      FloatVector y0 = new FloatVector(fact.Cols, rng);
      FloatVector y1 = new FloatVector(fact.Cols, rng);
      FloatVector x0 = fact.Solve(y0);
      FloatVector 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.
      FloatMatrix Y = new FloatMatrix(y1.Length, 2);
      Y.Col(0)[Slice.All] = y0;
      Y.Col(1)[Slice.All] = y1;
      FloatMatrix 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.
      FloatBandMatrix B = 1.2F * 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();
    }
  }
}