C# 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 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.
      DoubleBandMatrix 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 = {0}", B.ToString());

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

      // D = 1 2 7x7 [ 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 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.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      DoubleVector x = new DoubleVector(B.Cols, rng); // vector of random deviates
      DoubleVector y = MatrixFunctions.Product(B, x);

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

      DoubleBandMatrix BTB = MatrixFunctions.Product(B.Transpose(), B);
      Console.WriteLine();
      Console.WriteLine("BTB = {0}", BTB.ToString());

      // 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.CosFunction); // cosine of Cpi
      Console.WriteLine();
      Console.WriteLine("cos(Cpi) = {0}", Cpi.ToString());

      // 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();
      Console.WriteLine("Banded matrix product = {0}", P.ToString());

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


      // x and x2 should be the same. Let's look at the l2 norm of
      // their difference.
      DoubleVector residual = x - x2;
      double residualL2Norm = Math.Sqrt(NMathFunctions.Dot(residual, residual));
      Console.WriteLine();
      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.
      DoubleTriDiagMatrix T = new DoubleTriDiagMatrix(Cpi);

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