C# LU Factorization Example

using System;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# showing how to use the LU factorization class to solve linear
  /// systems, compute matrix inverses, condition numbers and determinants.
  /// </summary>
  class LUFactorizationExample
  {
    static void Main( string[] args )
    {
      var A = new DoubleMatrix( "3x3 [2 1 1  4 1 0 -2 2 1]" );
      var lu = new DoubleLUFact( A );

      Console.WriteLine();

      // Check to see if the input matrix is singular
      if ( lu.IsSingular )
      {
        Console.WriteLine( "Sorry, A is singular." );
        return;
      }

      // A is not singular, so we can look at the LU components,
      Console.WriteLine( "L..." );
      Console.WriteLine( lu.L.ToTabDelimited() );

      Console.WriteLine( "U..." );
      Console.WriteLine( lu.U.ToTabDelimited() );

      // and the permutation matrix P.
      Console.WriteLine( "P..." );
      Console.WriteLine( lu.P.ToTabDelimited() );

      // We can also compute the determinant and condition number.
      Console.WriteLine( "Determinant of A = {0}", lu.Determinant() );

      // We can choose to estimate the condition number (quick), or
      // compute it directly (accurate). For small matrices these
      // are usually the same.
      double estCond = lu.ConditionNumber( NormType.InfinityNorm, true );
      double computedCond = lu.ConditionNumber( NormType.OneNorm, false );
      Console.WriteLine( "Estimated condition number, infinity-norm = {0}", estCond.ToString( "F3" ) );
      Console.WriteLine( "Computed condition number, 1-norm = {0}", computedCond.ToString( "F3" ) );

      // Finally, we can compute the inverse of A
      Console.WriteLine();
      Console.WriteLine( "A inverse..." );
      Console.WriteLine( lu.Inverse().ToTabDelimited() );
      Console.WriteLine();

      // We can use the LU factorization to solve for one right-hand side,
      var v = new DoubleVector( "[8 11 3]" );
      DoubleVector u = lu.Solve( v );
      Console.WriteLine();
      Console.WriteLine( "The solution, u, Au=v is..." );
      Console.WriteLine( u );
      Console.WriteLine();

      // or we can solve for multiple right hand sides.
      var B = new DoubleMatrix( 3, 3 );
      // Set the all the columns of B to be the vector v.
      for ( int i = 0; i < 3; ++i )
      {
        B.Col( i )[Range.All] = v;
      }

      // Now solve for AX=B. The columns of A should have the same values as
      // the vector u from above.
      DoubleMatrix X = lu.Solve( B );
      Console.WriteLine( "The solution, X, to AX=B is..." );
      Console.WriteLine( X.ToTabDelimited() );

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

    } // Main

  }// class

}// namespace