C# 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
  /// tridiagonal matrices.
  /// </summary>
  class TriDiagFactExample
  {

    static void Main(string[] args)
    {
      // Construct a tridiagonal matrix with random entries.
      int rows = 5;
      int cols = 5;
      RandGenUniform rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      FloatComplexVector data1 = new FloatComplexVector( cols, rng );
      FloatComplexVector data2 = new FloatComplexVector( cols - 1, rng );
      FloatComplexVector data3 = new FloatComplexVector( cols - 1, rng );
      FloatComplexTriDiagMatrix A = new FloatComplexTriDiagMatrix( rows, cols );
      A.Diagonal()[Slice.All] = data1;
      A.Diagonal(1)[Slice.All] = data2;
      A.Diagonal(-1)[Slice.All] = data3;

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

      // A = 5x5 [ (-0.4974,0.3315) (0.5601,0.3060)  (0.0000,0.0000)  (0.0000,0.0000)   (0.0000,0.0000) 
      //           (0.7734,0.3580)  (-0.2503,0.5764) (0.2204,-0.0770) (0.0000,0.0000)   (0.0000,0.0000)
      //           (0.0000,0.0000)  (-0.8629,0.2029) (0.1961,0.1821)  (-0.1678,-0.2590) (0.0000,0.0000)  
      //           (0.0000,0.0000)  (0.0000,0.0000)  (-0.5848,0.6222) (-0.0443,0.0738)  (-0.9238,0.6206)
      //           (0.0000,0.0000)  (0.0000,0.0000)  (0.0000,0.0000)  (-0.7045,0.1236)  (-0.3249,-0.2797) ]

      // Construct a tridiagonal factorization class.
      FloatComplexTriDiagFact fact = new FloatComplexTriDiagFact( A );

      // Check to see if A is singular.
      string isSingularString = fact.IsSingular ? "A is singular" : "A is NOT singular";
      Console.WriteLine( isSingularString );

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