C# Triangular 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 triangular matrix classes.
  /// </summary>
  class TriangularMatrixExample
  {

    static void Main(string[] args)
    {
      // Set up the parameters that describe the shape of a Hermitian banded matrix.
      int order = 5;

      // Create upper and lower triangular matrices.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      DoubleMatrix A = new DoubleMatrix(order, order, rng);
      DoubleLUFact lu = new DoubleLUFact(A);
      DoubleUpperTriMatrix U = new DoubleUpperTriMatrix(lu.U);
      DoubleLowerTriMatrix L = new DoubleLowerTriMatrix(lu.L);

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

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

      // U = 5x5 [ 0.5764 -0.3249 -0.0770  0.7734  0.6222
      //           0.0000 -0.1691 -0.1416  0.0948 -0.9162  
      //           0.0000  0.0000  0.5759  0.3533  0.2067 
      //           0.0000  0.0000  0.0000 -1.4304 -1.0101 
      //           0.0000  0.0000  0.0000  0.0000 -0.5806 ]
      //
      // L = 5x5 [ 1.0000  0.0000 0.0000  0.0000 0.0000
      //           0.3403  1.0000 0.0000  0.0000 0.0000 
      //          -0.8630  0.5814 1.0000  0.0000 0.0000 
      //           0.5751 -0.8428 0.4012  1.0000 0.0000  
      //          -0.4343  0.3979 0.4225 -0.5379 1.0000 ]

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

      // You can set the values of elements in the upper triangular part
      // of a upper triangular matrix and the lower part of a lower
      // triangular matrix.
      double a = 99;
      L[2, 1] = a;
      Console.WriteLine("L[2,1] = {0}", L[2, 1]); // 99
      U[0, 2] = a + 1;
      Console.WriteLine("U[0,2] = {0}", U[0, 2]); // 100

      // But setting the values of elements in the lower triangular part
      // of a upper triangular matrix or the upper part of a lower
      // triangular matrix raises a NonModifiableElementException
      try
      {
        U[3, 0] = a;
      }
      catch (NonModifiableElementException e)
      {
        Console.WriteLine();
        Console.WriteLine("NonModifiableElementException: {0}", e.Message);
      }

      try
      {
        L[0, 3] = a;
      }
      catch (NonModifiableElementException e)
      {
        Console.WriteLine();
        Console.WriteLine("NonModifiableElementException: {0}", e.Message);
      }

      // Scalar multiplication and matrix addition/subtraction are supported.
      DoubleUpperTriMatrix C = a * U;
      DoubleUpperTriMatrix D = C + U;
      Console.WriteLine();
      Console.WriteLine("D = {0}", D.ToString("F4"));

      // Matrix/vector inner products too.
      DoubleVector x = new DoubleVector(L.Cols, rng);
      DoubleVector y = MatrixFunctions.Product(L, x);
      Console.WriteLine();
      Console.WriteLine("Lx = {0}", y.ToString());

      // You can transform the non-zero elements of a triangular matrix object
      // by using the Transform() method on its data vector.
      // Change every element of C to its absolute value.
      C.DataVector.Transform(NMathFunctions.AbsFunction);
      Console.WriteLine();
      Console.WriteLine("abs(C) = {0}", C.ToString("F4"));

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

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

      // You can calculate the determinant too.
      double det = MatrixFunctions.Determinant(U);
      Console.WriteLine();
      Console.WriteLine("Determinant of U = {0}", det);

      // You can use the Resize() method to change the bandwidths.
      D.Resize(6);
      Console.WriteLine();
      Console.WriteLine("D resized = {0}", D.ToString("F4"));

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