Band Matrices in C# QuickStart Sample

Illustrates how to work with the BandMatrix class in C#.

This sample is also available in: Visual Basic, F#, IronPython.

Overview

This QuickStart sample demonstrates how to work with band matrices in Numerics.NET using the BandMatrix class.

Band matrices are matrices whose elements are nonzero only in a diagonal band around the main diagonal. They appear frequently in scientific computing and numerical analysis, particularly when solving systems of linear equations that arise from finite difference methods.

The sample shows how to:

Create different types of band matrices including general, upper/lower triangular, and symmetric
Extract band matrices from existing matrices
Access and modify band matrix elements efficiently
Work with matrix properties and structural information
Use specialized decompositions for band matrices including LU and Cholesky
Solve linear systems involving band matrices
Handle row and column views of band matrices

The code illustrates the memory-efficient storage scheme used for band matrices and demonstrates how to properly allocate space for decompositions that require extra bandwidth.

The code

using System;

using Numerics.NET;
// The BandMatrix class resides in the Numerics.NET.LinearAlgebra
// namespace.
using Numerics.NET.LinearAlgebra;

// Illustrates the use of the BandMatrix class in the
// Numerics.NET.LinearAlgebra namespace of Numerics.NET.

// The license is verified at runtime. We're using
// a 30 day trial key here. For more information, see
//     https://numerics.net/trial-key
Numerics.NET.License.Verify("your-trial-key-here");

// Band matrices are matrices whose elements
// are nonzero only in a diagonal band around
// the main diagonal.
//
// General band matrices, upper and lower band
// matrices, and symmetric band matrices are all
// represented by a single class: BandMatrix.

//
// Constructing band matrices
//

// Constructing band matrices is similar to
// constructing general matrices. It is done by
// calling a factory method on the Matrix class.
// See theBasicMatrices QuickStart samples
// for a more complete discussion.

// The following creates a 7x5 band matrix with
// upper bandwidth 1 and lower bandwidth 2:
var b1 = Matrix.CreateBanded<double>(7, 5, 2, 1);

// Once the upper and lower bandwidth are set,
// it cannot be changed. Elements that are outside
// the band cannot be set.

// A second factory method lets you create upper
// or lower band matrices. The following constructs
// an 11x11 upper band matrix with unit diagonal
// and three non-zero upper diagonals.
var b2 = Matrix.CreateUpperBanded<double>(11, 11, 3, MatrixDiagonal.UnitDiagonal);

// To create a symmetric band matrix, you only need
// the size and the bandwith. The following creates
// a 6x6 symmetric tri-diagonal matrix:
var b3 = Matrix.CreateSymmetricBanded<double>(7, 1);

// We can assign values to the components by using
// the GetDiagonal method.
b3.GetDiagonal(0).SetValue(2);
b3.GetDiagonal(1).SetValue(-1);

// Extracting band matrices

// Another way to construct a band matrix is by
// extracting them from an existing matrix.
var m = Matrix.CreateFromArray(3, 4, new double[]
    {
        1, 2, 3,
        2, 3, 4,
        3, 4, 5,
        4, 5, 7
    }, MatrixElementOrder.ColumnMajor);
// To get the lower band part of m with bandwidth 2:
var b4 = BandMatrix<double>.Extract(m, 2, 0);

//
// BandMatrix properties
//

// A number of properties are available to determine
// whether a BandMatrix has a special structure:
Console.WriteLine($"b2 is upper? {b2.IsUpperTriangular}");
Console.WriteLine($"b2 is lower? {b2.IsUpperTriangular}");
Console.WriteLine($"b2 is unit diagonal? {b2.IsUnitDiagonal}");
Console.WriteLine($"b2 is symmetrical? {b2.IsSymmetrical}");

//
// BandMatrix methods
//

// You can get and set matrix elements:
b3[2, 3] = 55;
Console.WriteLine("b3[2, 3] = {0:F0}", b3[2, 3]);
// And the change will automatically be reflected
// in the symmetric element:
Console.WriteLine("b3[3, 2] = {0:F0}", b3[3, 2]);

//
// Row and column views
//

// The GetRow and GetColumn methods are
// available.
var row = b2.GetRow(1);
Console.WriteLine($"row 1 of b2 = {row:F0}");
var column = b2.GetColumn(2, 3, 4);
Console.WriteLine("column 3 of b2 from row 4 to ");
Console.WriteLine($"  row 5 = {column:F0}");

//
// Band matrix decompositions
//

// Specialized classes exist to represent the
// LU decomposition of a general band matrix
// and the Cholesky decomposition of a
// symmetric band matrix.

// Because of pivoting, the upper band matrix of
// the LU decomposition has larger bandwidth.
// You need to allocate extra space to be able to
// overwrite a matrix with its LU decomposition.

// The following creates a 7x5 band matrix with
// upper bandwidth 1 and lower bandwidth 2.
var b5 = Matrix.CreateBanded<double>(7, 7, 2, 1, true);
b5.GetDiagonal(0).SetValue(2.0);
b5.GetDiagonal(-2).SetValue(-1.0);
b5.GetDiagonal(1).SetValue(-1.0);

// Other than that, the API is the same as
// other decomposition classes.
var blu = b5.GetLUDecomposition(true);
var solution = blu.Solve(Vector.CreateConstant(b5.ColumnCount, 1.0));
Console.WriteLine($"  solution of b5*x = ones: {solution:F4}");

Console.Write("Press Enter key to exit...");
Console.ReadLine();