C# Band Fact Example

using System;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.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;
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      var data = new FloatVector( ( lowerBandwidth + upperBandwidth + 1 ) * cols, rng );
      var A = new FloatBandMatrix( data, rows, cols, lowerBandwidth, upperBandwidth );

      Console.WriteLine();
      
      Console.WriteLine( "A = " );
      
      // Display only three digits for added readability
      Console.WriteLine( A.ToTabDelimited( "G3" ) );

      // A =
      // -0.25   0.182   -0.325  0       0       0       0
      // 0.576   -0.0443 -0.28   0.22    0       0       0
      // 0       0.0738  0.56    -0.077  -0.924  0       0
      // 0       0       0.306   -0.168  0.621   -0.863  0
      // 0       0       0       -0.259  0.773   0.203   -0.705
      // 0       0       0       0       0.358   -0.585  0.124
      // 0       0       0       0       0       0.622   0.87

      // Construct a band factorization class.
      var 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();

      // In order to get condition number, factor with estimateCondition = true
      fact.Factor( A, true );
      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 =" );
      Console.WriteLine( AInv.ToTabDelimited( "G5" ) );

      // Use the factorization to solve some linear systems Ax = y.
      var y0 = new FloatVector( fact.Cols, rng );
      var 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( "G3" ) );

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

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

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

      // You can also solve for multiple right hand sides.
      var 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 =" );
      Console.WriteLine( X.ToTabDelimited( "G5" ) );

      // 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( "G3" ) );

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