C# Sym 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
  /// symmetric matrices.
  /// </summary>
  class SymFactExample
  {

    static void Main(string[] args)
    {
      // Construct a symmetric matrix as the product of the transpose of a 
      // matrix with itself.
      int rows = 5;
      int cols = 5;
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      DoubleMatrix A = new DoubleMatrix(rows, cols, rng);
      DoubleMatrix ATA = NMathFunctions.TransposeProduct(A, A);
      DoubleSymmetricMatrix S = new DoubleSymmetricMatrix(ATA);

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

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

      // Construct a symmetric factorization class.
      DoubleSymFact fact = new DoubleSymFact(S);

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

      // Retrieve information about the matrix S.
      double det = fact.Determinant();
      double rcond = fact.ConditionNumber();
      DoubleSymmetricMatrix SInv = fact.Inverse();

      Console.WriteLine();
      Console.WriteLine("Determinant of S = {0}", det);

      Console.WriteLine();
      Console.WriteLine("Reciprocal condition number = {0}", rcond);

      Console.WriteLine();
      Console.WriteLine("S inverse = {0}", SInv.ToString());

      // Use the factorization to solve some linear systems Ax = y.
      DoubleVector y0 = new DoubleVector(fact.Cols, rng);
      DoubleVector y1 = new DoubleVector(fact.Cols, rng);
      DoubleVector x0 = fact.Solve(y0);
      DoubleVector 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(S, x0)).ToString());

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

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

      // You can also solve for multiple right-hand sides.
      DoubleMatrix Y = new DoubleMatrix(y1.Length, 2);
      Y.Col(0)[Slice.All] = y0;
      Y.Col(1)[Slice.All] = y1;
      DoubleMatrix 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.
      DoubleSymmetricMatrix B = 1.2 * S;
      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();
    }
  }
}