C# Symmetric Matrix 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 symmetric matrix classes.
  /// </summary>
  class SymmetricMatrixExample
  {

    static void Main(string[] args)
    {
      int order = 6;

      // Set up a symmetric matrix S as the transpose product of a general 
      // matrix with itself (which is symmetric).
      RandGenUniform rng = new RandGenUniform(-2, 2);
      rng.Reset(0x124);
      FloatMatrix A = new FloatMatrix(order, order, rng);

      FloatSymmetricMatrix S = new FloatSymmetricMatrix(NMathFunctions.TransposeProduct(A, A));
      Console.WriteLine();
      Console.WriteLine("S = {0}", S.ToString());

      // S = 6x6 [  6.511227   0.7702899  2.189394  -2.113362   2.268176  -4.007094  
      //            0.7702899  3.701933  -2.814527  -0.8870333 -3.837288  -2.518311  
      //            2.189394  -2.814527   5.310741  -2.078353   3.270024  -1.440627  
      //           -2.113362  -0.8870333 -2.078353   6.054404   0.336647   5.635363
      //            2.268176  -3.837288   3.270024   0.336647   13.44877   2.772367  
      //           -4.007094  -2.518311  -1.440627   5.635363   2.772367   7.547361 ]

      // Indexer accessor works just like it does for general matrices. 
      Console.WriteLine();
      Console.WriteLine("S[2,2] = {0}", S[2, 2]);
      Console.WriteLine("S[3,0] = {0}", S[3, 0]);

      // You can set the values of elements in a symmetric matrix using the 
      // indexer. Note that setting the element in row i and column j implicitly
      // sets the element in column j and row i to the same value
      S[2, 1] = 100;
      Console.WriteLine("S[2,1] = {0}", S[2, 1]);
      Console.WriteLine("S[1,2] = {0}", S[1, 2]);

      // Scalar addition/subtractions/multiplication and matrix 
      // addition/subtraction are supported.
      float a = -.123F;
      FloatSymmetricMatrix C2 = a * S;
      FloatSymmetricMatrix C = 3 - S;
      FloatSymmetricMatrix D = C2 + S;
      Console.WriteLine();
      Console.WriteLine("D = {0}", D.ToString());

      // Matrix/vector products too.
      rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      FloatVector x = new FloatVector(S.Cols, rng); // vector of random deviates
      FloatVector y = MatrixFunctions.Product(S, x);
      Console.WriteLine();
      Console.WriteLine("Sx = {0}", y.ToString());

      // You can transform the elements of a symmetric matrix object by using
      // the Transform() method.
      C2.Transform(NMathFunctions.FloatPowFunction, 2); // Square every element of C2
      Console.WriteLine();
      Console.WriteLine("C2^2 = {0}", C2.ToString());

      // You can also solve linear systems.
      FloatVector x2 = MatrixFunctions.Solve(S, y);

      // x and x2 should be about the same. Let's look at the l2 norm of 
      // their difference.
      FloatVector residual = x - x2;
      float residualL2Norm = (float)Math.Sqrt(NMathFunctions.Dot(residual, residual));
      Console.WriteLine();
      Console.WriteLine("||x - x2|| = {0}", residualL2Norm);

      // Compute condition number
      float rcond = MatrixFunctions.ConditionNumber(S);
      Console.WriteLine();
      Console.WriteLine("Reciprocal condition number = {0}", rcond);

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