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

    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 an Hermitian banded matrix B by creating a banded complex matrix and setting
      // B equal to the product of the conjugate transpose of that matrix with itself.
      var A = new FloatComplexBandMatrix( order, order, halfBandwidth, halfBandwidth );
      for ( int i = -A.LowerBandwidth; i <= A.UpperBandwidth; ++i )
      {
        A.Diagonal( i ).Set( Slice.All, new FloatComplex( i + 2, i - 1 ) );
      }
      var B = new FloatHermitianBandMatrix( MatrixFunctions.Product( A.ConjTranspose(), A ) );

      Console.WriteLine();

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

      // B =
      // (10,0)  (10,6)  (3,6)   (0,0)   (0,0)   (0,0)   (0,0)
      // (10,-6) (19,0)  (10,6)  (3,6)   (0,0)   (0,0)   (0,0)
      // (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)   (0,0)   (0,0)
      // (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)   (0,0)
      // (0,0)   (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)  (3,6)
      // (0,0)   (0,0)   (0,0)   (3,-6)  (10,-6) (19,0)  (10,6)
      // (0,0)   (0,0)   (0,0)   (0,0)   (3,-6)  (10,-6) (14,0)

      // 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 an Hermitian 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 complex conjugate of that value.
      var z = new FloatComplex( 99, -99 );
      B[2, 1] = z;
      Console.WriteLine( "B[2,1] = {0}", B[2, 1] ); // (99, -99)
      Console.WriteLine( "B[1,2] = {0}", B[1, 2] ); // (99, 99 )

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

      // Scalar multiplication and matrix addition/subtraction are supported.
      z = new FloatComplex( .0001, .0001 );
      FloatHermitianBandMatrix C = z * B;
      FloatHermitianBandMatrix D = C - B;
      Console.WriteLine();
      Console.WriteLine( "D =" );
      Console.WriteLine( D.ToTabDelimited( "F5" ) );

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

      // You can transform the non-zero elements of a banded matrix object by using
      // the Transform() method on its data vector.
      C.DataVector.Transform( NMathFunctions.FloatComplexExpFunc );
      Console.WriteLine();
      Console.WriteLine( "exp(C) = " );
      Console.WriteLine( C.ToTabDelimited( "F5" ) );

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

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

      // You can calculate inverses too.
      FloatComplexMatrix BInv = MatrixFunctions.Inverse( B );
      Console.WriteLine( "BInv = " );
      Console.WriteLine( BInv.ToTabDelimited( "F5" ) );

      // If your Hermitian 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 FloatHermitianBandMatrix( MatrixFunctions.Product( A.ConjTranspose(), A ) );
      y = MatrixFunctions.Product( B, x );
      x2 = MatrixFunctions.Solve( B, y, true ); // 3rd parameter isPositiveDefinite set to true.

      // Lets see how close Bx is to y by computing the l2 norm of their difference.
      residual = x - x2;
      residualL2Norm = Math.Sqrt( NMathFunctions.ConjDot( residual, residual ).Real );
      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();
    }
  }
}