C# Triangular Matrix Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.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.
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      var A = new DoubleMatrix( order, order, rng );
      var lu = new DoubleLUFact( A );
      var U = new DoubleUpperTriMatrix( lu.U );
      var L = new DoubleLowerTriMatrix( lu.L );

      Console.WriteLine();

      Console.WriteLine( "U =" );
      Console.WriteLine( U.ToTabDelimited( "G3" ) );

      Console.WriteLine( "L =" );
      Console.WriteLine( L.ToTabDelimited( "G3" ) );

      // U =
      // 0.576   -0.325  -0.077  0.773   0.622
      // 0       -0.169  -0.142  0.0948  -0.916
      // 0       0       0.576   0.353   0.207
      // 0       0       0       -1.43   -1.01
      // 0       0       0       0       -0.581

      // L =
      // 1       0       0       0       0
      // 0.34    1       0       0       0
      // -0.863  0.581   1       0       0
      // 0.575   -0.843  0.401   1       0
      // -0.434  0.398   0.422   -0.538  1

      // 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 =" );
      Console.WriteLine( D.ToTabDelimited( "G3" ) );

      // Matrix/vector inner products too.
      var x = new DoubleVector( L.Cols, rng );
      DoubleVector y = MatrixFunctions.Product( L, x );
      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.AbsFunc );
      Console.WriteLine();
      Console.WriteLine( "abs(C) =" );
      Console.WriteLine( C.ToTabDelimited( "G3" ) );

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

      // x and x2 should be the same. Lets look at the l2 norm of
      // their difference.
      DoubleVector residual = x - x2;
      double residualL2Norm = Math.Sqrt( NMathFunctions.Dot( residual, residual ) );
      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 =" );
      Console.WriteLine( D.ToTabDelimited( "G3" ) );

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