diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs new file mode 100644 index 00000000..3ef77e28 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -0,0 +1,845 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Double +{ + using System; + using Properties; + using Threading; + + /// + /// Sparse Matrix implementation + /// + public class SparseMatrix : Matrix + { + /// + /// Object for use in "lock" + /// + private readonly object lockObject = new object(); + + /// + /// The array containing the row indices of the existing rows. Element "j" of the array gives the index of the + /// element in the array that is first non-zero element in a row "j" + /// + private readonly int[] rowIndex = new int[0]; + + /// + /// Array that contains the non-zero elements of matrix. Values of the non-zero elements of matrix are mapped into the values + /// array using the row-major storage mapping described in a compressed sparse row (CSR) format. + /// + private double[] nonZeroValues = new double[0]; + + /// + /// Gets the number of non zero elements in the matrix. + /// + /// The number of non zero elements. + public int NonZerosCount + { + get; + private set; + } + + /// + /// An array containing the column indices of the non-zero values. Element "I" of the array + /// is the number of the column in matrix that contains the I-th value in the array. + /// + private int[] columnIndices = new int[0]; + + /// + /// Initializes a new instance of the class. + /// + /// + /// The number of rows. + /// + /// + /// The number of columns. + /// + public SparseMatrix(int rows, int columns) : base(rows, columns) + { + this.rowIndex = new int[rows]; + } + + /// + /// Initializes a new instance of the class. This matrix is square with a given size. + /// + /// the size of the square matrix. + /// + /// If is less than one. + /// + public SparseMatrix(int order) : this(order, order) + { + } + + /// + /// Initializes a new instance of the class with all entries set to a particular value. + /// + /// + /// The number of rows. + /// + /// + /// The number of columns. + /// + /// The value which we assign to each element of the matrix. + public SparseMatrix(int rows, int columns, double value) : this(rows, columns) + { + if (value == 0.0) + { + return; + } + + this.NonZerosCount = rows * columns; + this.nonZeroValues = new double[this.NonZerosCount]; + this.columnIndices = new int[this.NonZerosCount]; + + for (int i = 0, j = 0; i < this.nonZeroValues.Length; i++, j++) + { + // Reset column position to "0" + if (j == columns) + { + j = 0; + } + + this.nonZeroValues[i] = value; + this.columnIndices[i] = j; + } + + // Set proper row pointers + for (var i = 0; i < this.rowIndex.Length; i++) + { + this.rowIndex[i] = ((i + 1) * columns) - columns; + } + } + + /// + /// Initializes a new instance of the class from a one dimensional array. + /// + /// The number of rows. + /// The number of columns. + /// The one dimensional array to create this matrix from. This array should store the matrix in column-major order. + /// If lenght is less than * . + /// + public SparseMatrix(int rows, int columns, double[] array) : this(rows, columns) + { + if (rows * columns > array.Length) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + for (var i = 0; i < rows; i++) + { + for (var j = 0; j < columns; j++) + { + this.SetValueAt(i, j, array[i + (j * rows)]); + } + } + } + + /// + /// Initializes a new instance of the class from a 2D array. + /// + /// The 2D array to create this matrix from. + public SparseMatrix(double[,] array) : this(array.GetLength(0), array.GetLength(1)) + { + var rows = array.GetLength(0); + var columns = array.GetLength(1); + + for (var i = 0; i < rows; i++) + { + for (var j = 0; j < columns; j++) + { + this.SetValueAt(i, j, array[i, j]); + } + } + } + + /// + /// Creates a SparseMatrix for the given number of rows and columns. + /// + /// + /// The number of rows. + /// + /// + /// The number of columns. + /// + /// + /// A SparseMatrix with the given dimensions. + /// + public override Matrix CreateMatrix(int numberOfRows, int numberOfColumns) + { + return new SparseMatrix(numberOfRows, numberOfColumns); + } + + /// + /// Creates a with a the given dimension. + /// + /// The size of the vector. + /// + /// A with the given dimension. + /// + public override Vector CreateVector(int size) + { + return new SparseVector(size); + } + + /// + /// Retrieves the requested element without range checking. + /// + /// + /// The row of the element. + /// + /// + /// The column of the element. + /// + /// + /// The requested element. + /// + public override double At(int row, int column) + { + lock (this.lockObject) + { + var index = this.FindItem(row, column); + return index >= 0 ? this.nonZeroValues[index] : 0.0; + } + } + + /// + /// Sets the value of the given element. + /// + /// + /// The row of the element. + /// + /// + /// The column of the element. + /// + /// + /// The value to set the element to. + /// + public override void At(int row, int column, double value) + { + lock (this.lockObject) + { + this.SetValueAt(row, column, value); + } + } + + #region Internal methods - CRS storage implementation + + /// + /// Created this method because we cannot call "virtual At" in constructor of the class, but we need to do it + /// + /// The row of the element. + /// The column of the element. + /// The value to set the element to. + /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks + private void SetValueAt(int row, int column, double value) + { + var index = this.FindItem(row, column); + if (index >= 0) + { + // Non-zero item found in matrix + if (value == 0.0) + { + // Delete existing item + this.DeleteItemByIndex(index, row); + } + else + { + // Update item + this.nonZeroValues[index] = value; + } + } + else + { + // Item not found. Add new value + if (value == 0.0) + { + return; + } + + index = ~index; + + // Check if the storage needs to be increased + if ((this.NonZerosCount == this.nonZeroValues.Length) && (this.NonZerosCount < (this.RowCount * this.ColumnCount))) + { + // Value array is completely full so we increase the size + // Determine the increase in size. We will not grow beyond the size of the matrix + var size = Math.Min(this.nonZeroValues.Length + this.GrowthSize(), this.RowCount * this.ColumnCount); + Array.Resize(ref this.nonZeroValues, size); + Array.Resize(ref this.columnIndices, size); + } + + // Move all values (with an position larger than index) in the value array to the next position + // move all values (with an position larger than index) in the columIndices array to the next position + for (var i = this.NonZerosCount - 1; i > index - 1; i--) + { + this.nonZeroValues[i + 1] = this.nonZeroValues[i]; + this.columnIndices[i + 1] = this.columnIndices[i]; + } + + // Add the value and the column index + this.nonZeroValues[index] = value; + this.columnIndices[index] = column; + + // increase the number of non-zero numbers by one + this.NonZerosCount += 1; + + // add 1 to all the row indices for rows bigger than rowIndex + // so that they point to the correct part of the value array again. + for (var i = row + 1; i < this.rowIndex.Length; i++) + { + this.rowIndex[i] += 1; + } + } + } + + /// + /// Delete value from internal storage + /// + /// Index of value in nonZeroValues array + /// Row number of matrix + /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks + private void DeleteItemByIndex(int itemIndex, int row) + { + // Move all values (with an position larger than index) in the value array to the previous position + // move all values (with an position larger than index) in the columIndices array to the previous position + for (var i = itemIndex + 1; i < this.NonZerosCount; i++) + { + this.nonZeroValues[i - 1] = this.nonZeroValues[i]; + this.columnIndices[i - 1] = this.columnIndices[i]; + } + + // Decrease value in Row + for (var i = row + 1; i < this.rowIndex.Length; i++) + { + this.rowIndex[i] -= 1; + } + + this.NonZerosCount -= 1; + + // Check if the storage needs to be shrink. This is reasonable to do if + // there are a lot of non-zero elements and storage is two times bigger + if ((this.NonZerosCount > 1024) && (this.NonZerosCount < this.nonZeroValues.Length / 2)) + { + Array.Resize(ref this.nonZeroValues, this.NonZerosCount); + Array.Resize(ref this.columnIndices, this.NonZerosCount); + } + } + + /// + /// Find item Index in nonZeroValues array + /// + /// Matrix row index + /// Matrix column index + /// Item index + /// WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks + private int FindItem(int row, int column) + { + // Determin bounds in columnIndices array where this item should be searched (using rowIndex) + var startIndex = this.rowIndex[row]; + var endIndex = row < this.rowIndex.Length - 1 ? this.rowIndex[row + 1] : this.NonZerosCount; + return Array.BinarySearch(this.columnIndices, startIndex, endIndex - startIndex, column); + } + + /// + /// Calculate grows size + /// + /// Proposed new size + private int GrowthSize() + { + int delta; + if (this.nonZeroValues.Length > 1024) + { + delta = this.nonZeroValues.Length / 4; + } + else + { + if (this.nonZeroValues.Length > 256) + { + delta = 512; + } + else + { + delta = this.nonZeroValues.Length > 64 ? 128 : 32; + } + } + + return delta; + } + + #endregion + + /// + /// Sets all values to zero. + /// + public override void Clear() + { + this.NonZerosCount = 0; + Array.Clear(this.rowIndex, 0, this.rowIndex.Length); + } + + /// + /// Copies the elements of this matrix to the given matrix. + /// + /// + /// The matrix to copy values into. + /// + /// + /// If target is . + /// + /// + /// If this and the target matrix do not have the same dimensions.. + /// + public override void CopyTo(Matrix target) + { + var sparseTarget = target as SparseMatrix; + + if (sparseTarget == null) + { + base.CopyTo(target); + } + else + { + if (ReferenceEquals(this, target)) + { + return; + } + + if (this.RowCount != target.RowCount || this.ColumnCount != target.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "target"); + } + + // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value + sparseTarget.nonZeroValues = new double[this.NonZerosCount]; + sparseTarget.columnIndices = new int[this.NonZerosCount]; + sparseTarget.NonZerosCount = this.NonZerosCount; + + Buffer.BlockCopy(this.nonZeroValues, 0, sparseTarget.nonZeroValues, 0, this.NonZerosCount * Constants.SizeOfDouble); + Buffer.BlockCopy(this.columnIndices, 0, sparseTarget.columnIndices, 0, this.NonZerosCount * Constants.SizeOfInt); + Buffer.BlockCopy(this.rowIndex, 0, sparseTarget.rowIndex, 0, this.RowCount * Constants.SizeOfInt); + } + } + + /// + /// Indicates whether the current object is equal to another object of the same type. + /// + /// + /// An object to compare with this object. + /// + /// + /// true if the current object is equal to the parameter; otherwise, false. + /// + public override bool Equals(object obj) + { + var sparseMatrix = obj as SparseMatrix; + + if (sparseMatrix == null) + { + return base.Equals(obj); + } + + // Accept if the argument is the same object as this. + if (ReferenceEquals(this, sparseMatrix)) + { + return true; + } + + if (this.ColumnCount != sparseMatrix.ColumnCount || this.RowCount != sparseMatrix.RowCount || this.NonZerosCount != sparseMatrix.NonZerosCount) + { + return false; + } + + // If all else fails, perform element wise comparison. + for (var index = 0; index < this.NonZerosCount; index++) + { + if (!this.nonZeroValues[index].AlmostEqual(sparseMatrix.nonZeroValues[index]) || this.columnIndices[index] != sparseMatrix.columnIndices[index]) + { + return false; + } + } + + return true; + } + + /// + /// Returns a hash code for this instance. + /// + /// + /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table. + /// + public override int GetHashCode() + { + var hashNum = Math.Min(this.NonZerosCount, 25); + long hash = 0; + for (var i = 0; i < hashNum; i++) + { +#if SILVERLIGHT + hash ^= Precision.DoubleToInt64Bits(this.nonZeroValues[i]); +#else + hash ^= BitConverter.DoubleToInt64Bits(this.nonZeroValues[i]); +#endif + } + + return BitConverter.ToInt32(BitConverter.GetBytes(hash), 4); + } + + /// + /// Returns the transpose of this matrix. + /// + /// The transpose of this matrix. + public override Matrix Transpose() + { + var ret = new SparseMatrix(this.ColumnCount, this.RowCount); + + // Do an 'inverse' CopyTo iterate over the rows + for (var i = 0; i < this.rowIndex.Length; i++) + { + // Get the begin / end index for the current row + var startIndex = this.rowIndex[i]; + var endIndex = i < this.rowIndex.Length - 1 ? this.rowIndex[i + 1] : this.NonZerosCount; + + // Get the values for the current row + if (startIndex == endIndex) + { + // Begin and end are equal. There are no values in the row, Move to the next row + continue; + } + + for (var j = startIndex; j < endIndex; j++) + { + ret[this.columnIndices[j], i] = this.nonZeroValues[j]; + } + } + + return ret; + } + + #region Elementary operations + + /// + /// Adds another matrix to this matrix. The result will be written into this matrix. + /// + /// The matrix to add to this matrix. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public override void Add(Matrix other) + { + if (ReferenceEquals(this, other)) + { + this.Multiply(2); + return; + } + + var m = other as SparseMatrix; + if (m == null) + { + base.Add(other); + } + else + { + this.Add(m); + } + } + + /// + /// Adds another to this matrix. The result will be written into this matrix. + /// + /// The to add to this matrix. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public void Add(SparseMatrix other) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (other.RowCount != this.RowCount || other.ColumnCount != this.ColumnCount) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + for (var i = 0; i < other.RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = other.rowIndex[i]; + var endIndex = i < other.rowIndex.Length - 1 ? other.rowIndex[i + 1] : other.NonZerosCount; + + for (var j = startIndex; j < endIndex; j++) + { + var index = this.FindItem(i, other.columnIndices[j]); + if (index >= 0) + { + if (this.nonZeroValues[index] + other.nonZeroValues[j] == 0.0) + { + this.DeleteItemByIndex(index, i); + } + else + { + this.nonZeroValues[index] += other.nonZeroValues[j]; + } + } + else + { + this.SetValueAt(i, other.columnIndices[j], other.nonZeroValues[j]); + } + } + } + } + + /// + /// Subtracts another matrix from this matrix. The result will be written into this matrix. + /// + /// The matrix to subtract. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public override void Subtract(Matrix other) + { + // We are substracting Matrix form itself + if (ReferenceEquals(this, other)) + { + this.Clear(); + return; + } + + var m = other as SparseMatrix; + if (m == null) + { + base.Subtract(other); + } + else + { + this.Subtract(m); + } + } + + /// + /// Subtracts another from this matrix. The result will be written into this matrix. + /// + /// The to subtract. + /// If the other matrix is . + /// If the two matrices don't have the same dimensions. + public void Subtract(SparseMatrix other) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (other.RowCount != this.RowCount || other.ColumnCount != this.ColumnCount) + { + throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); + } + + for (var i = 0; i < other.RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = other.rowIndex[i]; + var endIndex = i < other.rowIndex.Length - 1 ? other.rowIndex[i + 1] : other.NonZerosCount; + + for (var j = startIndex; j < endIndex; j++) + { + var index = this.FindItem(i, other.columnIndices[j]); + if (index >= 0) + { + if (this.nonZeroValues[index] - other.nonZeroValues[j] == 0.0) + { + this.DeleteItemByIndex(index, i); + } + else + { + this.nonZeroValues[index] -= other.nonZeroValues[j]; + } + } + else + { + this.SetValueAt(i, other.columnIndices[j], -other.nonZeroValues[j]); + } + } + } + } + + /// + /// Multiplies each element of this matrix with a scalar. + /// + /// The scalar to multiply with. + public override void Multiply(double scalar) + { + if (1.0.AlmostEqualInDecimalPlaces(scalar, 15)) + { + return; + } + + if (0.0.AlmostEqualInDecimalPlaces(scalar, 15)) + { + this.Clear(); + return; + } + + Control.LinearAlgebraProvider.ScaleArray(scalar, this.nonZeroValues); + } + + /// + /// Multiplies this sparse matrix with another sparse matrix and places the results into the result sparse matrix. + /// + /// The matrix to multiply with. + /// The result of the multiplication. + /// If the other matrix is . + /// If the result matrix is . + /// If this.Columns != other.Rows. + /// If the result matrix's dimensions are not the this.Rows x other.Columns. + public override void Multiply(Matrix other, Matrix result) + { + var otherSparseMatrix = other as SparseMatrix; + var resultSparseMatrix = result as SparseMatrix; + + if (otherSparseMatrix == null || resultSparseMatrix == null) + { + base.Multiply(other, result); + } + else + { + if (this.ColumnCount != otherSparseMatrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (resultSparseMatrix.RowCount != this.RowCount || resultSparseMatrix.ColumnCount != otherSparseMatrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + resultSparseMatrix.Clear(); + + var columnVector = new SparseVector(otherSparseMatrix.RowCount); + for (var row = 0; row < this.RowCount; row++) + { + // Get the begin / end index for the current row + var startIndex = this.rowIndex[row]; + var endIndex = row < this.rowIndex.Length - 1 ? this.rowIndex[row + 1] : this.NonZerosCount; + for (var column = 0; column < otherSparseMatrix.ColumnCount; column++) + { + columnVector.Clear(); + otherSparseMatrix.Column(column, columnVector); + + // Multiply row of matrix A on column of matrix B + var sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => this.nonZeroValues[index] * columnVector[this.columnIndices[index]]); + resultSparseMatrix.SetValueAt(row, column, sum); + } + } + } + } + + /// + /// Multiplies this matrix with another matrix and returns the result. + /// + /// The matrix to multiply with. + /// If this.Columns != other.Rows. + /// If the other matrix is . + /// The result of multiplication. + public override Matrix Multiply(Matrix other) + { + var matrix = other as SparseMatrix; + if (matrix == null) + { + return base.Multiply(other); + } + + if (this.ColumnCount != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + var result = (SparseMatrix)this.CreateMatrix(this.RowCount, matrix.ColumnCount); + Multiply(matrix, result); + return result; + } + + /// + /// Multiplies two sparse matrices. + /// + /// The left matrix to multiply. + /// The right matrix to multiply. + /// The result of multiplication. + /// If or is . + /// If the dimensions of or don't conform. + public static SparseMatrix operator *(SparseMatrix leftSide, SparseMatrix rightSide) + { + if (leftSide == null) + { + throw new ArgumentNullException("leftSide"); + } + + if (rightSide == null) + { + throw new ArgumentNullException("rightSide"); + } + + if (leftSide.ColumnCount != rightSide.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + return (SparseMatrix)leftSide.Multiply(rightSide); + } + + #endregion + + #region Static constructors for special matrices. + + /// + /// Initializes a square with all zero's except for ones on the diagonal. + /// + /// the size of the square matrix. + /// Identity SparseMatrix + /// + /// If is less than one. + /// + public static SparseMatrix Identity(int order) + { + var m = new SparseMatrix(order) + { + NonZerosCount = order, + nonZeroValues = new double[order], + columnIndices = new int[order] + }; + + for (var i = 0; i < order; i++) + { + m.nonZeroValues[i] = 1.0; + m.columnIndices[i] = i; + m.rowIndex[i] = i; + } + + return m; + } + + #endregion + } +} \ No newline at end of file diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs new file mode 100644 index 00000000..38d5d9e5 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs @@ -0,0 +1,170 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double +{ + using System; + using System.Collections.Generic; + using LinearAlgebra.Double; + using MbUnit.Framework; + + public class SparseMatrixTests : MatrixTests + { + protected override Matrix CreateMatrix(int rows, int columns) + { + return new SparseMatrix(rows, columns); + } + + protected override Matrix CreateMatrix(double[,] data) + { + return new SparseMatrix(data); + } + + protected override Vector CreateVector(int size) + { + return new SparseVector(size); + } + + protected override Vector CreateVector(double[] data) + { + return new SparseVector(data); + } + + [Test] + public void CanCreateMatrixFrom1DArray() + { + var testData = new Dictionary + { + { "Singular3x3", new SparseMatrix(3, 3, new double[] { 1, 1, 1, 1, 1, 1, 2, 2, 2 }) }, + { "Square3x3", new SparseMatrix(3, 3, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5, -3.3, 2.2, 6.6 }) }, + { "Square4x4", new SparseMatrix(4, 4, new[] { -1.1, 0.0, 1.0, -4.4, -2.2, 1.1, 2.1, 5.5, -3.3, 2.2, 6.2, 6.6, -4.4, 3.3, 4.3, -7.7 }) }, + { "Tall3x2", new SparseMatrix(3, 2, new[] { -1.1, 0.0, -4.4, -2.2, 1.1, 5.5 }) }, + { "Wide2x3", new SparseMatrix(2, 3, new[] { -1.1, 0.0, -2.2, 1.1, -3.3, 2.2 }) } + }; + + foreach (var name in testData.Keys) + { + Assert.AreEqual(this.testMatrices[name], testData[name]); + } + } + + [Test, Ignore] + public void MatrixFrom1DArrayIsReference() + { + // Sparse Matrix copies values from double[], but no remember reference. + var data = new double[] { 1, 1, 1, 1, 1, 1, 2, 2, 2 }; + var matrix = new SparseMatrix(3, 3, data); + matrix[0, 0] = 10.0; + Assert.AreEqual(10.0, data[0]); + } + + [Test] + public void MatrixFrom2DArrayIsCopy() + { + var matrix = new SparseMatrix(this.testData2D["Singular3x3"]); + matrix[0, 0] = 10.0; + Assert.AreEqual(1.0, this.testData2D["Singular3x3"][0, 0]); + } + + [Test] + [Row("Singular3x3")] + [Row("Singular3x3")] + [Row("Square3x3")] + [Row("Square4x4")] + [Row("Tall3x2")] + [Row("Wide2x3")] + public void CanCreateMatrixFrom2DArray(string name) + { + var matrix = new SparseMatrix(this.testData2D[name]); + for (var i = 0; i < this.testData2D[name].GetLength(0); i++) + { + for (var j = 0; j < this.testData2D[name].GetLength(1); j++) + { + Assert.AreEqual(this.testData2D[name][i, j], matrix[i, j]); + } + } + } + + [Test] + public void CanCreateMatrixWithUniformValues() + { + var matrix = new SparseMatrix(10, 10, 10.0); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreEqual(matrix[i, j], 10.0); + } + } + } + + [Test] + public void CanCreateIdentity() + { + var matrix = SparseMatrix.Identity(5); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreEqual(i == j ? 1.0 : 0.0, matrix[i, j]); + } + } + } + + [Test] + [Row(0)] + [Row(-1)] + [ExpectedArgumentException] + public void IdentityFailsWithZeroOrNegativeOrder(int order) + { + var matrix = SparseMatrix.Identity(order); + } + + [Test] + public void CanCreateHeavySparseMatrix() + { + var matrix = new SparseMatrix(500, 1000); + var nonzero = 0; + var rnd = new Random(0); + + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + var value = rnd.NextDouble(); + if (value != 0) + { + nonzero++; + } + + matrix[i, j] = value; + } + } + + Assert.AreEqual(matrix.NonZerosCount, nonzero); + } + } +} \ No newline at end of file