C# Herm PD Tri Diag 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
  /// 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;
      var data1 = new DoubleComplexVector( cols, 2 );
      var data2 = new DoubleComplexVector( cols - 1, -1 );
      var 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 = " );
      Console.WriteLine( A.ToTabDelimited( "G3" ) );

      // 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.
      var 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();

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

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

      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 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 = " );
      Console.WriteLine( X.ToTabDelimited( "G3" ) );

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

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

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