C# Tridiagonal Matrix Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the tridiagonal matrix classes.
  /// </summary>
  class TridiagonalMatrixExample
  {

    static void Main( string[] args )
    {
      // Set up the parameters that describe the shape of a tridiagonal matrix.
      int rows = 8;
      int cols = 8;

      // Set up a tridiagonal matrix B by setting all the diagonals within the matrix 
      // bandwidth.
      var B = new FloatComplexTriDiagMatrix( rows, cols );
      for ( int i = -1; i <= 1; ++i )
      {
        B.Diagonal( i ).Set( Slice.All, i + 2 );
      }

      Console.WriteLine();

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

      // B =
      // (2,0)   (3,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
      // (1,0)   (2,0)   (3,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
      // (0,0)   (1,0)   (2,0)   (3,0)   (0,0)   (0,0)   (0,0)   (0,0)
      // (0,0)   (0,0)   (1,0)   (2,0)   (3,0)   (0,0)   (0,0)   (0,0)
      // (0,0)   (0,0)   (0,0)   (1,0)   (2,0)   (3,0)   (0,0)   (0,0)
      // (0,0)   (0,0)   (0,0)   (0,0)   (1,0)   (2,0)   (3,0)   (0,0)
      // (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (1,0)   (2,0)   (3,0)
      // (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (1,0)   (2,0)

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

      // You can set the values of elements in the main, super, and sub diagonals
      // of a tridiagonal matrix using the indexer.
      B[2, 1] = new FloatComplex( -100, 99 );
      Console.WriteLine( "B[2,1] = {0}", B[2, 1] );

      // But setting an element that would destroy the tridiagonal structure
      // of the matrix raises a NonModifiableElementException exception.
      try
      {
        B[7, 0] = 21;
      }
      catch ( NonModifiableElementException e )
      {
        Console.WriteLine();
        Console.WriteLine( "NonModifiableElementException: {0}", e.Message );
      }

      // Scalar addition/subtractions/multiplication and matrix 
      // addition/subtraction are supported.
      float s = -.123F;
      FloatComplexTriDiagMatrix C2 = s * B;
      FloatComplexTriDiagMatrix C = 3 - B;
      FloatComplexTriDiagMatrix D = C2 + B;
      Console.WriteLine();
      Console.WriteLine( "D =" );
      Console.WriteLine( D.ToTabDelimited( "F2" ) );

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

      // You can transform the non-zero elements of a banded matrix object by using
      // the Transform() method on its data vector.
      C2.DataVector.Transform( NMathFunctions.FloatComplexPowFunc, 2 ); // Square every element of C2
      Console.WriteLine();
      Console.WriteLine( "C2^2 =" );
      Console.WriteLine( C2.ToTabDelimited( "F2" ));

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

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

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

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