My C# Implementation Of Basic Linear Algebra Concepts

Overview

 
Today, I will be sharing with you my C# implementation of basic linear algebra concepts. This code has been posted to GitHub under a MIT license, so feel free to modify and deal with code without any restrictions or limitations (no guarantees of any kind.) And please let me know your feedback, comments, suggestions, and corrections.
 

Introduction

 
The library has been built using C# and has been posted to GitHub ang NuGet.
 
It covers the following,
 
Matrices
  • Matrix addition / subtraction.
  • Matrix concatenation / shrinking / extraction.
  • Matrix reduction and elimination (Gauss and Gauss-Jordan).
  • Matrix rank, solution, and solution state (unique, infinite, and no solution.)
  • Determinant calculation using Laprace Expansion method.
  • Invertibility and inverses.
  • Matrix mirrors.
  • Matrix multiplication and powers.
  • Transposes.
  • Identity matrices.
  • Transformation matrices.
  • Projection matrices.
  • Reflection matrices.
  • 2D and 3D rotation matrices.
  • 2D and 3D shearing matrices.
  • Cloning and equality comparers.
Vectors
  • Vector addition / subtraction.
  • Angel between vectors.
  • Vector normalization and magnitude calculator.
  • Dot product and cross product.
  • Vector scaling.
  • Projection onto other vectors and subspaces.
  • Cloning and equality comparers.
The design is based on raw functions that operate on Array objects. Wrapper classes for matrices and vectors (Matrix and Vector, respectively) have been created to hide those keen functions. Static initializers and operators have been added to simplify the use of the wrapper classes.
 
I will not go into a discussion of how a function mathematically works. Despite this, feel free to correct me whenever you want.
 
Let’s have a brief look at some of the highlights of the code.
 

Matrices

 
Matrix Elimination
 
This function demonstrates Gaussian and Gauss-Jordan elimination. It can be used to reduce a matrix to the row-echelon form (REF / Gaussian) or reduced row-echelon form (RREF / Gauss-Jordan.) It can be used to solve a matrix for a number of augmented columns and to calculate matrix rank.
  1. /// <summary>  
  2. /// Reduces matrix to row-echelon (REF/Gauss) or reduced row-echelon (RREF/Gauss-Jordan) form and solves for augmented columns (if any.)  
  3. /// </summary>  
  4. public static MatrixEliminationResult Eliminate(double[,] input, MatrixReductionForm form, int augmentedCols = 0) {  
  5.   int totalRowCount = input.GetLength(0);  
  6.   int totalColCount = input.GetLength(1);  
  7.   
  8.   if (augmentedCols >= totalColCount)  
  9.     throw new ArgumentException(nameof(augmentedCols),  
  10.       Properties.Resources.Exception_TooMuchAugmentedColumns);  
  11.   
  12.   MatrixEliminationResult result = new MatrixEliminationResult();  
  13.   
  14.   double[,] output = CreateCopy(input);  
  15.   
  16.   // number of pivots found  
  17.   int numPivots = 0;  
  18.   
  19.   // loop through columns, exclude augmented columns  
  20.   for (int col = 0; col < totalColCount - augmentedCols; col++) {  
  21.     int? pivotRow = FindPivot(output, numPivots, col, totalRowCount);  
  22.   
  23.     if (pivotRow == null)  
  24.       continue// no pivots, go to another column  
  25.   
  26.     ReduceRow(output, pivotRow.Value, col, totalColCount);  
  27.   
  28.     SwitchRows(output, pivotRow.Value, numPivots, totalColCount);  
  29.   
  30.     pivotRow = numPivots;  
  31.     numPivots++;  
  32.   
  33.     // Eliminate Previous Rows  
  34.     if (form == MatrixReductionForm.ReducedRowEchelonForm) {  
  35.       for (int tmpRow = 0; tmpRow < pivotRow; tmpRow++)  
  36.         EliminateRow(output, tmpRow, pivotRow.Value, col, totalColCount);  
  37.     }  
  38.   
  39.     // Eliminate Next Rows  
  40.     for (int tmpRow = pivotRow.Value; tmpRow < totalRowCount; tmpRow++)  
  41.       EliminateRow(output, tmpRow, pivotRow.Value, col, totalColCount);  
  42.   
  43.   }  
  44.   
  45.   result.Rank = numPivots;  
  46.   result.FullMatrix = output;  
  47.   result.AugmentedColumnCount = augmentedCols;  
  48.   if (augmentedCols > 0 && form == MatrixReductionForm.ReducedRowEchelonForm) { // matrix has solution   
  49.     result.Solution = ExtractColumns(output, totalColCount - augmentedCols, totalColCount - 1);  
  50.   }  
  51.   
  52.   return result;  
  53. }  
  54.   
  55. private static int? FindPivot(double[,] input, int startRow, int col, int rowCount) {  
  56.   for (int i = startRow; i < rowCount; i++) {  
  57.     if (input[i, col] != 0)  
  58.       return i;  
  59.   }  
  60.   
  61.   return null;  
  62. }  
  63.   
  64. private static void SwitchRows(double[,] input, int row1, int row2, int colCount) {  
  65.   if (row1 == row2)  
  66.     return;  
  67.   
  68.   for (int col = 0; col < colCount; col++) {  
  69.     var tmpVal1 = input[row1, col];  
  70.     input[row1, col] = input[row2, col];  
  71.     input[row2, col] = tmpVal1;  
  72.   }  
  73. }  
  74.   
  75. private static void ReduceRow(double[,] input, int row, int col, int colCount) {  
  76.   var coeffecient = 1.0 / input[row, col];  
  77.   
  78.   if (coeffecient == 1)  
  79.     return;  
  80.   
  81.   for (; col < colCount; col++)  
  82.     input[row, col] *= coeffecient;  
  83. }  
  84.   
  85. /// <summary>  
  86. /// Eliminates row using another pivot row.  
  87. /// </summary>  
  88. private static void EliminateRow(double[,] input, int row, int pivotRow, int pivotCol, int colCount) {  
  89.   if (pivotRow == row)  
  90.     return;  
  91.   
  92.   if (input[row, pivotCol] == 0)  
  93.     return;  
  94.   
  95.   double coeffecient = input[row, pivotCol];  
  96.   
  97.   for (int col = pivotCol; col < colCount; col++) {  
  98.     input[row, col] -= input[pivotRow, col] * coeffecient;  
  99.   }  
  100. }   
This function has been wrapped into our Matrix wrapper with those functions:
  1. public virtual Matrix Reduce(MatrixReductionForm form, int? augmentedColCount = null);  
  2. public virtual MatrixSolutionState GetSolutionState(int? augmentedColCount = null);  
  3. public virtual Matrix Solve(int? augmentedColCount = null);  
  4. public virtual int GetRank(int? augmentedColCount = null);   
Matrix Determinant
 
This function determinant uses the Laprace Expansion method. It works by calculating the cross-product of the matrix. Here’s the implementation,
  1. /// <summary>  
  2. /// Calculates determinant. Internally uses Laprace Expansion method.  
  3. /// </summary>  
  4. /// <remarks>  
  5. /// Returns 1 for an empty matrix. See https://math.stackexchange.com/questions/1762537/what-is-the-determinant-of/1762542  
  6. /// </remarks>  
  7. public static double Determinant(double[,] input) {  
  8.   var results = CrossProduct(input);  
  9.   
  10.   return results.Sum();  
  11. }  
  12.   
  13. public static double[] CrossProduct(double[,] input) {  
  14.   int rowCount = input.GetLength(0);  
  15.   int colCount = input.GetLength(1);  
  16.   
  17.   if (rowCount != colCount)  
  18.     throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);  
  19.   
  20.   if (rowCount == 0)  
  21.     return new double[] { 1 };  
  22.   
  23.   if (rowCount == 1)  
  24.     return new double[] { input[0, 0] };  
  25.   
  26.   if (rowCount == 2)  
  27.     return new double[] { (input[0, 0] * input[1, 1]) - (input[0, 1] * input[1, 0]) };  
  28.   
  29.   double[] results = new double[colCount];  
  30.   
  31.   for (int col = 0; col < colCount; col++) {  
  32.     // checkerboard pattern, even col  = 1, odd col = -1  
  33.     var checkerboardFactor = col % 2.0 == 0 ? 1 : -1;  
  34.     var coeffecient = input[0, col];  
  35.   
  36.     var crossMatrix = GetCrossMatrix(input, 0, col);  
  37.     results[col] = checkerboardFactor * (coeffecient * Determinant(crossMatrix));  
  38.   }  
  39.   
  40.   return results;  
  41. }  
  42.   
  43. /// <summary>  
  44. /// Retrieves all matrix entries except the specified row and col. Used in cross product and determinant functions.  
  45. /// </summary>  
  46. public static double[,] GetCrossMatrix(double[,] input, int skipRow, int skipCol) {  
  47.   int rowCount = input.GetLength(0);  
  48.   int colCount = input.GetLength(1);  
  49.   
  50.   var output = new double[rowCount - 1, colCount - 1];  
  51.   int outputRow = 0;  
  52.   
  53.   for (int row = 0; row < rowCount; row++) {  
  54.     if (row == skipRow)  
  55.       continue;  
  56.   
  57.     int outputCol = 0;  
  58.   
  59.     for (int col = 0; col < colCount; col++) {  
  60.       if (col == skipCol)  
  61.         continue;  
  62.   
  63.       output[outputRow, outputCol] = input[row, col];  
  64.   
  65.       outputCol++;  
  66.     }  
  67.     outputRow++;  
  68.   }  
  69.   
  70.   return output;  
  71. }   
Matrix Inverses
 
The matrix inverse function eliminates a function and determines if all columns are unique by checking matrix rank. Another function checks if the matrix is invertible by checking determinant.
  1. /// <summary>  
  2. /// Returns a value indicates whether the specified matrix is invertible. Internally uses array determinant.  
  3. /// </summary>  
  4. public static bool IsInvertible(double[,] input) {  
  5.   var rowCount = input.GetLength(0);  
  6.   var colCount = input.GetLength(1);  
  7.   
  8.   if (rowCount != colCount)  
  9.     return false;  
  10.   
  11.   return Determinant(input) != 0;  
  12. }  
  13.   
  14. /// <summary>  
  15. /// Calculates the inverse of a matrix. Returns null if non-invertible.  
  16. /// </summary>  
  17. public static double[,] Invert(double[,] matrix) {  
  18.   var rowCount = matrix.GetLength(0);  
  19.   var colCount = matrix.GetLength(1);  
  20.   
  21.   if (rowCount != colCount)  
  22.     throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);  
  23.   
  24.   var newMatrix = ConcatHorizontally(matrix, CreateIdentityMatrix(rowCount));  
  25.   
  26.   var result = Eliminate(newMatrix, MatrixReductionForm.ReducedRowEchelonForm, rowCount);  
  27.   if (result.Rank != colCount)  
  28.     return null;  
  29.   
  30.   return result.Solution;  
  31. }   
Matrix Multiplication
 
Four variations of the Multiply function are available: scalar multiplication, matrix multiplications, unsafe multiplications, and powers.
  1. /// <summary>  
  2. /// Multiplies/Scales a matrix by a scalar input.  
  3. /// </summary>  
  4. public static double[,] Multiply(double scalar, double[,] input) {  
  5.   var rowCount = input.GetLength(0);  
  6.   var colCount = input.GetLength(1);  
  7.   
  8.   // creating the final product matrix  
  9.   double[,] output = new double[rowCount, colCount];  
  10.   
  11.   for (int row = 0; row < rowCount; row++) {  
  12.     for (int col = 0; col < colCount; col++) {  
  13.       output[row, col] = scalar * input[row, col];  
  14.     }  
  15.   }  
  16.   
  17.   return output;  
  18. }  
  19.   
  20. /// <summary>  
  21. /// Multiplies two matrices together.  
  22. /// </summary>  
  23. public static double[,] Multiply(double[,] matrix1, double[,] matrix2) {  
  24.   // cahing matrix lengths for better performance  
  25.   var matrix1Rows = matrix1.GetLength(0);  
  26.   var matrix1Cols = matrix1.GetLength(1);  
  27.   var matrix2Rows = matrix2.GetLength(0);  
  28.   var matrix2Cols = matrix2.GetLength(1);  
  29.   
  30.   // checking if product is defined  
  31.   if (matrix1Cols != matrix2Rows)  
  32.     throw new InvalidOperationException(Properties.Resources.Exception_UndefinedProduct);  
  33.   
  34.   // creating the final product matrix  
  35.   double[,] output = new double[matrix1Rows, matrix2Cols];  
  36.   
  37.   // looping through matrix 1 rows  
  38.   for (int matrix1_row = 0; matrix1_row < matrix1Rows; matrix1_row++) {  
  39.     // for each matrix 1 row, loop through matrix 2 columns  
  40.     for (int matrix2_col = 0; matrix2_col < matrix2Cols; matrix2_col++) {  
  41.       // loop through matrix 1 columns to calculate the dot product  
  42.       for (int matrix1_col = 0; matrix1_col < matrix1Cols; matrix1_col++) {  
  43.         output[matrix1_row, matrix2_col] += matrix1[matrix1_row, matrix1_col] * matrix2[matrix1_col, matrix2_col];  
  44.       }  
  45.     }  
  46.   }  
  47.   
  48.   return output;  
  49. }  
  50.   
  51. /// <summary>  
  52. /// Multiplies two matrices together. Uses unsafe code. Better with extremely large matrices.  
  53. /// </summary>  
  54. public static double[,] MultiplyUnsafe(double[,] matrix1, double[,] matrix2) {  
  55.   // cahing matrix lengths for better performance  
  56.   var matrix1Rows = matrix1.GetLength(0);  
  57.   var matrix1Cols = matrix1.GetLength(1);  
  58.   var matrix2Rows = matrix2.GetLength(0);  
  59.   var matrix2Cols = matrix2.GetLength(1);  
  60.   
  61.   // checking if product is defined  
  62.   if (matrix1Cols != matrix2Rows)  
  63.     throw new InvalidOperationException(Properties.Resources.Exception_UndefinedProduct);  
  64.   
  65.   // creating the final product matrix  
  66.   double[,] output = new double[matrix1Rows, matrix2Cols];  
  67.   
  68.   unsafe  
  69.   {  
  70.     // fixing pointers to matrices  
  71.     fixed (  
  72.       double* pProduct = output,  
  73.       pMatrix1 = matrix1,  
  74.       pMatrix2 = matrix2) {  
  75.   
  76.       int i = 0;  
  77.       // looping through matrix 1 rows  
  78.       for (int matrix1_row = 0; matrix1_row < matrix1Rows; matrix1_row++) {  
  79.         // for each matrix 1 row, loop through matrix 2 columns  
  80.         for (int matrix2_col = 0; matrix2_col < matrix2Cols; matrix2_col++) {  
  81.           // loop through matrix 1 columns to calculate the dot product  
  82.           for (int matrix1_col = 0; matrix1_col < matrix1Cols; matrix1_col++) {  
  83.   
  84.             var val1 = *(pMatrix1 + (matrix1Rows * matrix1_row) + matrix1_col);  
  85.             var val2 = *(pMatrix2 + (matrix2Cols * matrix1_col) + matrix2_col);  
  86.   
  87.             *(pProduct + i) += val1 * val2;  
  88.   
  89.           }  
  90.   
  91.           i++;  
  92.   
  93.         }  
  94.       }  
  95.   
  96.     }  
  97.   }  
  98.   
  99.   return output;  
  100. }  
  101.   
  102. /// <summary>  
  103. /// Raises an input matrix to the nth power.  
  104. /// </summary>  
  105. public static double[,] Power(double[,] input, int power) {  
  106.   if (input.GetLength(0) != input.GetLength(1))  
  107.     throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);  
  108.   if (power < 0)  
  109.     throw new ArgumentException(nameof(power));  
  110.   if (power == 0)  
  111.     return CreateIdentityMatrix(input.GetLength(0));  
  112.   
  113.   double[,] output = CreateCopy(input);  
  114.   for (int i = 0; i < power; i++) {  
  115.     output = Multiply(output, input);  
  116.   }  
  117.   
  118.   return output;  
  119. }  
Subspace Projection Matrices
 
The following function creates projection matrix from a subspace.
  1. /// <summary>  
  2. /// Creates projection matrix for the specified subspace.  
  3. /// </summary>  
  4. public static double[,] CreateProjectionMatrix(double[,] subspace) {  
  5.   var subspaceTranspose = MatrixFunctions.Transpose(subspace);  
  6.   
  7.   double[,] value = MatrixFunctions.Multiply(subspaceTranspose, subspace);  
  8.   
  9.   value = MatrixFunctions.Invert(value);  
  10.   
  11.   value = MatrixFunctions.Multiply(value, subspaceTranspose);  
  12.   
  13.   value = MatrixFunctions.Multiply(subspace, value);  
  14.   
  15.   return value;  
  16. }   
Reflection Matrices
 
The following functions create 2D and 3D reflection matrices,
  1. public static double[,] Create2DReflectionMatrix(MatrixAxis axis) {  
  2.   if (axis == MatrixAxis.Z)  
  3.     throw new InvalidOperationException(Properties.Resources.Exception_InvalidAxis);  
  4.   
  5.   var output = CreateIdentityMatrix(2);  
  6.   
  7.   if (axis == MatrixAxis.X)  
  8.     output[1, 1] *= -1;  
  9.   else  if (axis == MatrixAxis.Y )  
  10.     output[0, 0] *= -1;  
  11.   
  12.   return output;  
  13. }  
  14.   
  15. public static double[,] Create3DReflectionMatrix(Matrix3DReflectionPlane axis) {  
  16.   var output = CreateIdentityMatrix(4);  
  17.   
  18.   if (axis == Matrix3DReflectionPlane.XY)  
  19.     output[2, 2] *= -1;  
  20.   
  21.   else if (axis == Matrix3DReflectionPlane.YZ)  
  22.     output[0,0] *= -1;  
  23.   
  24.   else if (axis == Matrix3DReflectionPlane.ZX)  
  25.     output[1,1] *= -1;  
  26.   
  27.   return output;  
  28. }  
Rotation Matrices
 
The following functions create 2D and 3D rotation matrices. They accept angle, unit (radians / degrees) and direction (clockwise / counter-clockwise),
  1. /// <summary>  
  2. /// Creates 2-dimensional rotation matrix using the specified angle.  
  3. /// </summary>  
  4. public static double[,] Create2DRotationMatrix(double angle, AngleUnit unit, MatrixRotationDirection direction) {  
  5.   // sin and cos accept only radians  
  6.   
  7.   double angleRadians = angle;  
  8.   if (unit == AngleUnit.Degrees)  
  9.     angleRadians = Converters.DegreesToRadians(angleRadians);  
  10.   
  11.   
  12.   double[,] output = new double[2, 2];  
  13.   
  14.   output[0, 0] = Math.Cos(angleRadians);  
  15.   output[1, 0] = Math.Sin(angleRadians);  
  16.   output[0, 1] = Math.Sin(angleRadians);  
  17.   output[1, 1] = Math.Cos(angleRadians);  
  18.   
  19.   if (direction == MatrixRotationDirection.Clockwise)  
  20.     output[1, 0] *= -1;  
  21.   else  
  22.     output[0, 1] *= -1;  
  23.   
  24.   
  25.   return output;  
  26. }  
  27.   
  28. /// <summary>  
  29. /// Creates 2-dimensional rotation matrix using the specified angle and axis.  
  30. /// </summary>  
  31. public static double[,] Create3DRotationMatrix(double angle, AngleUnit unit, MatrixAxis axis) {  
  32.   // sin and cos accept only radians  
  33.   
  34.   double angleRadians = angle;  
  35.   if (unit == AngleUnit.Degrees)  
  36.     angleRadians = Converters.DegreesToRadians(angleRadians);  
  37.   
  38.   
  39.   double[,] output = new double[3, 3];  
  40.   
  41.   if (axis == MatrixAxis.X) {  
  42.     output[0, 0] = 1;  
  43.     output[1, 1] = Math.Cos(angleRadians);  
  44.     output[2, 1] = Math.Sin(angleRadians);  
  45.     output[1, 2] = -1 * Math.Sin(angleRadians);  
  46.     output[2, 2] = Math.Cos(angleRadians);  
  47.   } else if (axis == MatrixAxis.Y) {  
  48.     output[1, 1] = 1;  
  49.     output[0, 0] = Math.Cos(angleRadians);  
  50.     output[2, 0] = -1 * Math.Sin(angleRadians);  
  51.     output[0, 2] = Math.Sin(angleRadians);  
  52.     output[2, 2] = Math.Cos(angleRadians);  
  53.   } else if (axis == MatrixAxis.Z) {  
  54.     output[2, 2] = 1;  
  55.     output[0, 0] = Math.Cos(angleRadians);  
  56.     output[1, 0] = Math.Sin(angleRadians);  
  57.     output[0, 1] = -1 * Math.Sin(angleRadians);  
  58.     output[1, 1] = Math.Cos(angleRadians);  
  59.   }  
  60.   
  61.   
  62.   return output;  
  63. }   
Shearing Matrices
 
The following functions create shearing matrices using a specific factor over the given axis,
  1. public static double[,] Create2DShearingMatrix(double factor, MatrixAxis axis) {  
  2.   if (axis == MatrixAxis.Z)  
  3.     throw new InvalidOperationException(Properties.Resources.Exception_InvalidAxis);  
  4.   
  5.   var output = CreateIdentityMatrix(2);  
  6.   
  7.   if (axis ==  MatrixAxis.X)  
  8.     output[0, 1] = factor;  
  9.   
  10.   else if (axis == MatrixAxis.Y )  
  11.     output[1, 0] = factor;  
  12.   
  13.   return output;  
  14. }  
  15.   
  16. public static double[,] Create3DShearingMatrix(double factor, MatrixAxis axis) {  
  17.   var output = CreateIdentityMatrix(4);  
  18.   
  19.   if (axis == MatrixAxis.X) {   
  20.     output[1, 0] = factor;  
  21.     output[2, 0] = factor;  
  22.   } else if (axis == MatrixAxis.Y) {  
  23.     output[0, 1] = factor;  
  24.     output[2, 1] = factor;  
  25.   } else if (axis == MatrixAxis.Z) {  
  26.     output[0, 2] = factor;  
  27.     output[1, 2] = factor;  
  28.   }  
  29.   
  30.   return output;  
  31. }   

Vectors

 
The following lists some of the major vector functions covered in this library. We are having a look at the raw functions. Wrappers will be covered later.
 
Angels
 
This function calculates the angle between two vectors. It relies on two other functions, DotProduct and GetMagnitude.
  1. public static double AngleBetween(double[] vector1, double[] vector2, AngleUnit unit) {  
  2.   if (vector1.Length != vector2.Length)  
  3.     throw new InvalidOperationException(Properties.Resources.Exception_DimensionsMismatch);  
  4.   
  5.   var dotProduct = DotProduct(vector1, vector2);  
  6.   var lengthProduct = GetMagnitude(vector1) * GetMagnitude(vector2);  
  7.   
  8.   var result = Math.Acos(dotProduct / lengthProduct);  
  9.   if (unit == AngleUnit.Degrees)  
  10.     result = Converters.RadiansToDegrees(result);  
  11.   
  12.   return result;  
  13. }   
Normalization
 
The following code demonstrates two operations: normalization (i.e. converting to unit vector) and magnitude (length) calculator.
  1. /// <summary>  
  2. /// Normalizes a vector.  
  3. /// </summary>  
  4. public static double[] ToUnitVector(double[] input) {  
  5.   var length = input.Length;  
  6.   
  7.   double[] output = new double[length];  
  8.   var coeffecient = 1.0 / GetMagnitude(input);  
  9.   
  10.   for (int i = 0; i < length; i++) {  
  11.     output[i] = coeffecient * input[i];  
  12.   }  
  13.   
  14.   return output;  
  15. }  
  16.   
  17.   
  18. public static double GetMagnitude(double[] input) {  
  19.   double val = 0;  
  20.   
  21.   for (int i = 0; i < input.Length; i++)  
  22.     val += input[i] * input[i];  
  23.   
  24.   val = Math.Sqrt(val);  
  25.   
  26.   return val;  
  27. }   
Dot-Product and Cross-Product
 
The following code demonstrates scaling, dot-product, and cross product.
  1. public static double DotProduct(double[] vector1, double[] vector2) {  
  2.   if (vector1.Length != vector2.Length)  
  3.     throw new InvalidOperationException(Properties.Resources.Exception_DimensionsMismatch);  
  4.   
  5.   double product = 0;  
  6.   for (int i = 0; i < vector1.Length; i++)  
  7.     product += vector1[i] * vector2[i];  
  8.   
  9.   return product;  
  10. }  
  11.   
  12. public static double[] Scale(double scalar, double[] vector) {  
  13.   double[] product = new double[vector.Length];  
  14.   
  15.   for (int i = 0; i < vector.Length; i++) {  
  16.     product[i] = scalar * vector[i];  
  17.   }  
  18.   
  19.   return product;  
  20. }  
  21.   
  22.   
  23. public static double[] CrossProduct(double[] vector1, double[] vector2) {  
  24.   int length = 3;  
  25.   
  26.   if (length != vector1.Length || length != vector2.Length)  
  27.     throw new InvalidOperationException(Properties.Resources.Exception_3DRequired);  
  28.   
  29.   double[,] matrix = new double[length, length];  
  30.   for (int row = 0; row < length; row++) {  
  31.     for (int col = 0; col < length; col++) {  
  32.       if (row == 0)  
  33.         matrix[row, col] = 1;  
  34.       else if (row == 1)  
  35.         matrix[row, col] = vector1[col];  
  36.       else if (row == 2)  
  37.         matrix[row, col] = vector2[col];  
  38.     }  
  39.   }  
  40.   
  41.   
  42.   return MatrixFunctions.CrossProduct(matrix);  
  43. }  
Projection
 
Projection of a vector onto another vector is handled using the following functions:
  1. public static double ProjectionFactor(double[] vector1, double[] vector2) {  
  2.   return DotProduct(vector1, vector2) / DotProduct(vector2, vector2);  
  3. }  
  4. public static double[] Projection(double[] vector1, double[] vector2) {  
  5.   var factor = ProjectionFactor(vector1, vector2);  
  6.   return Scale(factor, vector2);  
  7. }   

Wrappers

 
Working with raw Array objects especially multi-dimensional arrays is fairly cumbersome. To make it easier, we wrapped those raw functions into two wrapper classes, Matrix and Vector. Both classes implement ICloneable and IEquatable. Both have indexers and a direct way to access their inner representation (i.e. array).
 
Properties and Methods
 
The following figures show the exported members of both classes.
 
 
 
Operators
 
We overrided some operators for both classes to allow clean and clear access to their operations. Here’s the overriden operator list for Matrix:
  1. public static Matrix operator +(Matrix m) { return m; }  
  2. public static Matrix operator -(Matrix m) { return m.Scale(-1); }  
  3. public static Matrix operator +(Matrix a, Matrix b) { return a.Add(b); }  
  4. public static Matrix operator -(Matrix a, Matrix b) { return a.Subtract(b); }  
  5. public static Matrix operator +(Matrix a, double[,] b) { return a.Add(b); }  
  6. public static Matrix operator -(Matrix a, double[,] b) { return a.Subtract(b); }  
  7.   
  8. public static Matrix operator *(double a, Matrix m) { return m.Scale(a); }  
  9. public static Matrix operator *(Matrix a, Matrix b) { return a.Multiply(b); }  
  10. public static Matrix operator *(Matrix a, double[,] b) { return a.Multiply(b); }  
  11.   
  12. public static bool operator ==(Matrix a, Matrix b) {  
  13.   if ((a as object) == null || (b as object) == null)  
  14.     return false;  
  15.   return a.Equals(b);  
  16. }  
  17. public static bool operator !=(Matrix a, Matrix b) {  
  18.   if ((a as object) == null || (b as object) == null)  
  19.     return false;  
  20.   return a.Equals(b) == false;  
  21. }  
  22.   
  23. public static Matrix operator ^(Matrix a, int power) { return a.Power(power); }  
  24.   
  25.   
  26. public static implicit operator double[,] (Matrix m) { return m.InnerMatrix; }  
  27. public static explicit operator Matrix(double[,] m) { return new Linears.Matrix(m); }   

Last Word

 
As stated previously, the code is using a MIT license. So please use the code without any restrictions, limitations, or warranties, and feel free to relay your comments, suggestions, and corrections.


Similar Articles