C# Hermitian Matrix Example

using System;

using CenterSpace.NMath.Core;


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

    static void Main( string[] args )
    {
      int order = 5;

      // Set up a Hermitian matrix S as the conjugate transpose product of a general 
      // matrix with itself (which is Hermitian).
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      var A = new DoubleComplexMatrix( order, order, rng );

      var S = new DoubleHermitianMatrix( NMathFunctions.ConjTransposeProduct( A, A ) );

      Console.WriteLine();

      Console.WriteLine( "S = " );
      Console.WriteLine( S.ToTabDelimited( "F5" ) );

      // S =
      // (3.12186,0.00000)       (0.12935,0.76321)       (-0.59263,-0.51912)     (1.01693,-0.48541)      (-0.62109,-0.74390)
      // (0.12935,-0.76321)      (1.01859,0.00000)       (-0.61581,0.58225)      (0.34714,-1.17980)      (-0.37649,0.35263)
      // (-0.59263,0.51912)      (-0.61581,-0.58225)     (2.66911,0.00000)       (-0.78612,2.23106)      (0.09417,0.18527)
      // (1.01693,0.48541)       (0.34714,1.17980)       (-0.78612,-2.23106)     (4.10411,0.00000)       (0.59635,0.69605)
      // (-0.62109,0.74390)      (-0.37649,-0.35263)     (0.09417,-0.18527)      (0.59635,-0.69605)      (3.97634,0.00000)

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

      // You can set the values of elements in a Hermitian 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.
      S[2, 1] = new DoubleComplex( 100, -99 );
      Console.WriteLine( "S[2,1] = {0}", S[2, 1] ); // (100, -99)
      Console.WriteLine( "S[1,2] = {0}", S[1, 2] ); // (100, 99)

      // Scalar multiplication and matrix addition/subtraction are supported.
      Double a = -.123;
      DoubleHermitianMatrix C2 = a * S;
      DoubleHermitianMatrix D = C2 + S;
      Console.WriteLine();
      Console.WriteLine( "D = " );
      Console.WriteLine( D.ToTabDelimited( "F5" ) );

      // Matrix/vector products too.
      var x = new DoubleComplexVector( S.Cols, rng ); // vector of random deviates
      DoubleComplexVector y = MatrixFunctions.Product( S, x );
      Console.WriteLine( "Sx = {0}", y.ToString( "G3" ) );

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

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

      // You can transform the elements of a Hermitian matrix object by using
      // the Transform() method.
      C2.DataVector.Transform( NMathFunctions.DoubleComplexCoshFunc );
      Console.WriteLine();
      Console.WriteLine( "cosh(C2) = " );
      Console.WriteLine( C2.ToTabDelimited( "G5" ) );

      // For a matrix to satisfy the strict definition of a Hermitian matrix,
      // its diagonal elements must be real. The Hermitian matrix classes provide
      // a MakeDigaonalReal() method to ensure that your matrix satisfies the
      // the strict definition of Hermitian.
      C2.MakeDiagonalReal();
      Console.WriteLine( "Diagonal element is real: {0}", C2[3, 3].Imag == 0.0 ); // True

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

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