C# LU Factorization Example

using System;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.Core.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 LUFactorizationExampple
  {
    static void Main(string[] args)
    {
      Console.WriteLine();

      DoubleMatrix A = new DoubleMatrix("3x3 [2 1 1  4 1 0 -2 2 1]");
      DoubleLUFact lu = new DoubleLUFact(A);

      // 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,
      DoubleVector 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.
      DoubleMatrix 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