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

    static void Main( string[] args )
    {
      // Set up the parameters that describe the shape of a banded matrix.
      int upperBandwidth = 2;
      int lowerBandwidth = 1;
      int rows = 7;
      int cols = 7;

      // Set up a banded matrix B by setting all the diagonals within the matrix 
      // bandwidth.
      var B = new DoubleBandMatrix( rows, cols, lowerBandwidth, upperBandwidth );
      for ( int i = -B.LowerBandwidth; i <= B.UpperBandwidth; ++i )
      {
        B.Diagonal( i ).Set( Slice.All, i + 2 );
      }

      Console.WriteLine();

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

      // B =
      // 2       3       4       0       0       0       0
      // 1       2       3       4       0       0       0
      // 0       1       2       3       4       0       0
      // 0       0       1       2       3       4       0
      // 0       0       0       1       2       3       4
      // 0       0       0       0       1       2       3
      // 0       0       0       0       0       1       2

      // 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 banded matrix using the indexer.
      B[2, 1] = 99;
      Console.WriteLine( "B[2,1] = {0}", B[2, 1] );

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

      // Scalar multiplication and addition/subtraction are supported.
      DoubleBandMatrix C = 3 * B;
      DoubleBandMatrix D = C - B;

      Console.WriteLine();
      Console.WriteLine( "D =" );
      Console.WriteLine( D.ToTabDelimited() );

      // D =
      // 4       6       8       0       0       0       0
      // 2       4       6       8       0       0       0
      // 0       198     4       6       8       0       0
      // 0       0       2       4       6       8       0
      // 0       0       0       2       4       6       8
      // 0       0       0       0       2       4       6
      // 0       0       0       0       0       2       4

      // Matrix/vector, Matrix/Matrix products too.
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      var x = new DoubleVector( B.Cols, rng ); // vector of random deviates
      DoubleVector y = MatrixFunctions.Product( B, x );

      Console.WriteLine( "Bx = " );
      Console.WriteLine( y.ToString( "G5" ) );

      DoubleBandMatrix BTB = MatrixFunctions.Product( B.Transpose(), B );
      Console.WriteLine();
      Console.WriteLine( "BTB = " );
      Console.WriteLine( BTB.ToTabDelimited( "G5" ) );

      // You can transform the non-zero elements of a banded matrix object by using
      // the Transform() method on its data vector.
      DoubleBandMatrix Cpi = Math.PI * C;
      Cpi.DataVector.Transform( NMathFunctions.CosFunc ); // cosine of Cpi
      Console.WriteLine( "cos(Cpi) =" );
      Console.WriteLine( Cpi.ToTabDelimited( "G5" ) );

      // Since the inner product of two banded matrices is a banded matrix,
      // matrix inner products are supported.
      DoubleBandMatrix P = MatrixFunctions.Product( Cpi, C );
      Console.WriteLine( "Banded matrix product =" );
      Console.WriteLine( P.ToTabDelimited( "G5" ) );

      // 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 use the Resize() method to change the bandwidths.
      Cpi.Resize( Cpi.Rows, Cpi.Cols, 1, 1 ); // matrix is now tridiagonal.

      // And construct a tridiagonal matrix from it.
      var T = new DoubleTriDiagMatrix( Cpi );

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