C# Symmetric Band Matrix Example

using System;

using CenterSpace.NMath.Core;


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

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

      // Set up a symmetric banded matrix B by creating a banded matrix, A, 
      // with random entries, and setting B equal to the product of the transpose 
      // of that matrix with itself.
      var A = new DoubleBandMatrix( order, order, halfBandwidth, halfBandwidth );
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      DoubleVector diag;

      // Fill the non-zero entries of A, a diagonal at a time, with random numbers.
      for ( int i = -A.LowerBandwidth; i <= A.UpperBandwidth; ++i )
      {
        diag = A.Diagonal( i );
        for ( int j = 0; j < diag.Length; ++j )
        {
          diag[j] = rng.Next();
        }
      }
      var B = new DoubleSymBandMatrix( MatrixFunctions.TransposeProduct( A, A ) );

      Console.WriteLine();

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

      // B =
      // 0.249   -0.0333 0.0835  0       0       0       0
      // -0.0333 0.121   -0.12   -0.0858 0       0       0
      // 0.0835  -0.12   0.196   0.154   0.231   0       0
      // 0       -0.0858 0.154   0.477   0.581   0.358   0
      // 0       0       0.231   0.581   1.21    0.408   0.152
      // 0       0       0       0.358   0.408   0.512   0.277
      // 0       0       0       0       0.152   0.277   0.647

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

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

      // But setting an element outside the bandwidth of the
      // matrix raises a NonModifiableElementException exception
      try
      {
        B[6, 0] = a;
      }
      catch ( NonModifiableElementException e )
      {
        Console.WriteLine();
        Console.WriteLine( "NonModifiableElementException: {0}", e.Message );
      }

      // Scalar multiplication and matrix addition/subtraction are supported.
      DoubleSymBandMatrix C = a * B;
      DoubleSymBandMatrix D = C + B;
      Console.WriteLine();
      Console.WriteLine( "D =" );
      Console.WriteLine( D.ToTabDelimited( "G3" ) );

      // Matrix/vector inner products too.
      var x = new DoubleVector( B.Cols, rng );
      DoubleVector 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.
      // Change every element of C to its natural log.
      C.DataVector.Transform( NMathFunctions.LogFunc );
      Console.WriteLine();
      Console.WriteLine( "ln(C) =" );
      Console.WriteLine( C.ToTabDelimited( "G3" ) );

      // You can also solve linear systems.
      DoubleVector x2 = MatrixFunctions.Solve( B, 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( B );
      Console.WriteLine();
      Console.WriteLine( "Determinant of B = {0}", det );

      // If your symmetric banded matrix is positive definite, you can invoke
      // the Solve function with a third - the isPositiveDefinite 
      // parameter - set to true. 
      // First make sure B is positive definite.
      B = new DoubleSymBandMatrix( MatrixFunctions.TransposeProduct( A, A ) );
      y = MatrixFunctions.Product( B, x );
      x2 = MatrixFunctions.Solve( B, y, true ); // 3rd parameter isPositiveDefinite set to true.

      // See how close Bx is to y by computing the l2 norm of their difference.
      residual = x - x2;
      residualL2Norm = Math.Sqrt( NMathFunctions.Dot( residual, residual ) );
      Console.WriteLine();
      Console.WriteLine( "PD ||x - x2|| = {0}", residualL2Norm );

      // You can use the Resize() method to change the bandwidths.
      D.Resize( D.Order, 2 );

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