using System; using CenterSpace.NMath.Core; namespace CenterSpace.NMath.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. var 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 ) ); } var B = new FloatHermitianBandMatrix( MatrixFunctions.Product( A.ConjTranspose(), A ) ); Console.WriteLine(); Console.WriteLine( "B = " ); Console.WriteLine( B.ToTabDelimited( "G3" ) ); // B = // (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. var 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 =" ); Console.WriteLine( D.ToTabDelimited( "F5" ) ); // Matrix/vector inner products. var rng = new RandGenUniform( -1, 1 ); rng.Reset( 0x124 ); var 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.FloatComplexExpFunc ); Console.WriteLine(); Console.WriteLine( "exp(C) = " ); Console.WriteLine( C.ToTabDelimited( "F5" ) ); // You can also solve linear systems. FloatComplexVector x2 = MatrixFunctions.Solve( B, y ); // x and x2 should be the same. Lets 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 ); Console.WriteLine(); // You can calculate inverses too. FloatComplexMatrix BInv = MatrixFunctions.Inverse( B ); Console.WriteLine( "BInv = " ); Console.WriteLine( BInv.ToTabDelimited( "F5" ) ); // 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. // Lets 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(); } } }← All NMath Code Examples