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

    static void Main(string[] args)
    {
      // Set up the parameters that describe the shape of a tridiagonal matrix.
      int rows = 8;
      int cols = 8;

      // Set up a tridiagonal matrix B by setting all the diagonals within the matrix 
      // bandwidth.
      FloatComplexTriDiagMatrix B = new FloatComplexTriDiagMatrix(rows, cols);
      for (int i = -1; i <= 1; ++i)
      {
        B.Diagonal(i).Set(Slice.All, i + 2);
      }
      Console.WriteLine();
      Console.WriteLine("B = {0}", B.ToString());

      // B = 8x8 [ (2,0) (3,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)  
      //           (1,0) (2,0) (3,0) (0,0) (0,0) (0,0) (0,0) (0,0)  
      //           (0,0) (1,0) (2,0) (3,0) (0,0) (0,0) (0,0) (0,0)  
      //           (0,0) (0,0) (1,0) (2,0) (3,0) (0,0) (0,0) (0,0)  
      //           (0,0) (0,0) (0,0) (1,0) (2,0) (3,0) (0,0) (0,0)  
      //           (0,0) (0,0) (0,0) (0,0) (1,0) (2,0) (3,0) (0,0)  
      //           (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (2,0) (3,0)  
      //           (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (1,0) (2,0)]

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

      // You can set the values of elements in the main, super, and sub diagonals
      // of a tridiagonal matrix using the indexer.
      B[2, 1] = new FloatComplex(-100, 99);
      Console.WriteLine("B[2,1] = {0}", B[2, 1]);

      // But setting an element that would destroy the tridiagonal structure
      // of the matrix raises a NonModifiableElementException exception.
      try
      {
        B[7, 0] = 21;
      }
      catch (NonModifiableElementException e)
      {
        Console.WriteLine();
        Console.WriteLine("NonModifiableElementException: {0}", e.Message);
      }

      // Scalar addition/subtractions/multiplication and matrix 
      // addition/subtraction are supported.
      float s = -.123F;
      FloatComplexTriDiagMatrix C2 = s * B;
      FloatComplexTriDiagMatrix C = 3 - B;
      FloatComplexTriDiagMatrix D = C2 + B;
      Console.WriteLine();
      Console.WriteLine("D = {0}", D.ToString());

      // Matrix/vector products too.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      FloatComplexVector x = new FloatComplexVector(B.Cols, rng); // vector of random deviates
      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.
      C2.DataVector.Transform(NMathFunctions.FloatComplexPowFunction, 2); // Square every element of C2
      Console.WriteLine();
      Console.WriteLine("C2^2 = {0}", C2.ToString());

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

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

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

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