C# Hermitian Band Matrix Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.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.
      FloatComplexBandMatrix 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));
      }
      FloatHermitianBandMatrix B =
        new FloatHermitianBandMatrix(MatrixFunctions.Product(A.ConjTranspose(), A));

      Console.WriteLine();
      Console.WriteLine("B = {0}", B.ToString());
      Console.WriteLine();

      // B = 2 7x7 [ (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.
      FloatComplex 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 = {0}", D.ToString());

      // Matrix/vector inner products.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      FloatComplexVector 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.FloatComplexExpFunction);
      Console.WriteLine();
      Console.WriteLine("exp(C) = {0}", C.ToString());

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

      // x and x2 should be the same. Let's 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);

      // You can calculate inverses too.
      FloatComplexMatrix BInv = MatrixFunctions.Inverse(B);
      Console.WriteLine();
      Console.WriteLine("BInv = {0}", BInv.ToString());

      // 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.

      // Let's 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();
    }
  }
}