diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index a810bd2f..64b4eeaf 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex using System.Numerics; using Generic; using Properties; + using Storage; using Threading; /// @@ -44,17 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex [Serializable] public class SparseMatrix : Matrix { - /// - /// 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 Complex[] _nonZeroValues = new Complex[0]; + readonly SparseCompressedRowMatrixStorage _storage; /// /// Gets the number of non zero elements in the matrix. @@ -62,16 +53,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The number of non zero elements. public int NonZerosCount { - get; - private set; + get { return _storage.ValueCount; } } - /// - /// 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. /// @@ -84,7 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public SparseMatrix(int rows, int columns) : base(rows, columns) { - _rowIndex = new int[rows]; + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex.Zero); } /// @@ -95,8 +79,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If is less than one. /// public SparseMatrix(int order) - : this(order, order) + : base(order, order) { + _storage = new SparseCompressedRowMatrixStorage(order, order, Complex.Zero); } /// @@ -110,18 +95,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The value which we assign to each element of the matrix. public SparseMatrix(int rows, int columns, Complex value) - : this(rows, columns) + : base(rows, columns) { - if (value == 0.0) + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex.Zero); + + if (value.IsZero()) { return; } - NonZerosCount = rows * columns; - _nonZeroValues = new Complex[NonZerosCount]; - _columnIndices = new int[NonZerosCount]; + var rowPointers = _storage.RowPointers; + var valueCount = _storage.ValueCount = rows * columns; + var columnIndices = _storage.ColumnIndices = new int[valueCount]; + var values = _storage.Values = new Complex[valueCount]; - for (int i = 0, j = 0; i < _nonZeroValues.Length; i++, j++) + for (int i = 0, j = 0; i < values.Length; i++, j++) { // Reset column position to "0" if (j == columns) @@ -129,14 +117,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex j = 0; } - _nonZeroValues[i] = value; - _columnIndices[i] = j; + values[i] = value; + columnIndices[i] = j; } // Set proper row pointers - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { - _rowIndex[i] = ((i + 1) * columns) - columns; + rowPointers[i] = ((i + 1) * columns) - columns; } } @@ -149,8 +137,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If length is less than * . /// public SparseMatrix(int rows, int columns, Complex[] array) - : this(rows, columns) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex.Zero); + if (rows * columns > array.Length) { throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); @@ -160,7 +150,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, array[i + (j * rows)]); + _storage.SetValueAt(i, j, array[i + (j * rows)]); } } } @@ -170,16 +160,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The 2D array to create this matrix from. public SparseMatrix(Complex[,] array) - : this(array.GetLength(0), array.GetLength(1)) + : base(array.GetLength(0), array.GetLength(1)) { - var rows = array.GetLength(0); - var columns = array.GetLength(1); + _storage = new SparseCompressedRowMatrixStorage(array.GetLength(0), array.GetLength(1), Complex.Zero); - for (var i = 0; i < rows; i++) + for (var i = 0; i < _storage.RowCount; i++) { - for (var j = 0; j < columns; j++) + for (var j = 0; j < _storage.ColumnCount; j++) { - SetValueAt(i, j, array[i, j]); + _storage.SetValueAt(i, j, array[i, j]); } } } @@ -190,8 +179,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// /// The matrix to copy. public SparseMatrix(Matrix matrix) - : this(matrix.RowCount, matrix.ColumnCount) + : base(matrix.RowCount, matrix.ColumnCount) { + _storage = new SparseCompressedRowMatrixStorage(matrix.RowCount, matrix.ColumnCount, Complex.Zero); + var sparseMatrix = matrix as SparseMatrix; var rows = matrix.RowCount; @@ -203,23 +194,28 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, matrix.At(i, j)); + _storage.SetValueAt(i, j, matrix.At(i, j)); } } } else { - NonZerosCount = sparseMatrix.NonZerosCount; - _rowIndex = new int[rows]; - _columnIndices = new int[NonZerosCount]; - _nonZeroValues = new Complex[NonZerosCount]; + var matrixStorage = sparseMatrix.Storage; + var valueCount = _storage.ValueCount = matrixStorage.ValueCount; + _storage.ColumnIndices = new int[valueCount]; + _storage.Values = new Complex[valueCount]; - Array.Copy(sparseMatrix._nonZeroValues, _nonZeroValues, NonZerosCount); - Array.Copy(sparseMatrix._columnIndices, _columnIndices, NonZerosCount); - Array.Copy(sparseMatrix._rowIndex, _rowIndex, rows); + Array.Copy(matrixStorage.Values, _storage.Values, valueCount); + Array.Copy(matrixStorage.ColumnIndices, _storage.ColumnIndices, valueCount); + Array.Copy(matrixStorage.RowPointers, _storage.RowPointers, rows); } } + internal SparseCompressedRowMatrixStorage Storage + { + get { return _storage; } + } + /// /// Creates a SparseMatrix for the given number of rows and columns. /// @@ -297,15 +293,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Where to store the lower triangle. private void LowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row >= _columnIndices[j]) + if (row >= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -359,15 +360,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Where to store the lower triangle. private void UpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row <= _columnIndices[j]) + if (row <= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -426,17 +432,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var result = (SparseMatrix)CreateMatrix(rowCount, columnCount); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (int i = rowIndex, row = 0; i < rowMax; i++, row++) { - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; + for (int j = startIndex; j < endIndex; j++) { // check if the column index is in the range - if ((_columnIndices[j] >= columnIndex) && (_columnIndices[j] < columnIndex + columnCount)) + if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount)) { - var column = _columnIndices[j] - columnIndex; - result.SetValueAt(row, column, _nonZeroValues[j]); + var column = columnIndices[j] - columnIndex; + result._storage.SetValueAt(row, column, values[j]); } } } @@ -493,15 +505,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Where to store the lower triangle. private void StrictlyLowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row > _columnIndices[j]) + if (row > columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -556,15 +573,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Where to store the lower triangle. private void StrictlyUpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row < _columnIndices[j]) + if (row < columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -581,13 +603,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// An array containing the matrix's elements. public override Complex[] ToColumnWiseArray() { + var values = _storage.Values; var ret = new Complex[RowCount * ColumnCount]; for (var j = 0; j < ColumnCount; j++) { for (var i = 0; i < RowCount; i++) { - var index = FindItem(i, j); - ret[(j * RowCount) + i] = index >= 0 ? _nonZeroValues[index] : 0.0; + var index = _storage.FindItem(i, j); + ret[(j * RowCount) + i] = index >= 0 ? values[index] : 0.0; } } @@ -608,8 +631,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// public override Complex At(int row, int column) { - var index = FindItem(row, column); - return index >= 0 ? _nonZeroValues[index] : 0.0; + return _storage.GetValueAt(row, column); } /// @@ -626,166 +648,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// public override void At(int row, int column, Complex value) { - 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, Complex value) - { - var index = FindItem(row, column); - if (index >= 0) - { - // Non-zero item found in matrix - if (value == 0.0) - { - // Delete existing item - DeleteItemByIndex(index, row); - } - else - { - // Update item - _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 ((NonZerosCount == _nonZeroValues.Length) && (NonZerosCount < ((long)RowCount * 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(_nonZeroValues.Length + GrowthSize(), (long)RowCount * ColumnCount); - if (size > int.MaxValue) - { - throw new NotSupportedException(Resources.TooManyElements); - } - - Array.Resize(ref _nonZeroValues, (int)size); - Array.Resize(ref _columnIndices, (int)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 = NonZerosCount - 1; i > index - 1; i--) - { - _nonZeroValues[i + 1] = _nonZeroValues[i]; - _columnIndices[i + 1] = _columnIndices[i]; - } - - // Add the value and the column index - _nonZeroValues[index] = value; - _columnIndices[index] = column; - - // increase the number of non-zero numbers by one - 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 < _rowIndex.Length; i++) - { - _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 < NonZerosCount; i++) - { - _nonZeroValues[i - 1] = _nonZeroValues[i]; - _columnIndices[i - 1] = _columnIndices[i]; - } - - // Decrease value in Row - for (var i = row + 1; i < _rowIndex.Length; i++) - { - _rowIndex[i] -= 1; - } - - 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 ((NonZerosCount > 1024) && (NonZerosCount < _nonZeroValues.Length / 2)) - { - Array.Resize(ref _nonZeroValues, NonZerosCount); - Array.Resize(ref _columnIndices, NonZerosCount); - } + _storage.SetValueAt(row, column, value); } - /// - /// 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 = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; - return Array.BinarySearch(_columnIndices, startIndex, endIndex - startIndex, column); - } - - /// - /// Calculates the amount with which to grow the storage array's if they need to be - /// increased in size. - /// - /// The amount grown. - private int GrowthSize() - { - int delta; - if (_nonZeroValues.Length > 1024) - { - delta = _nonZeroValues.Length / 4; - } - else - { - if (_nonZeroValues.Length > 256) - { - delta = 512; - } - else - { - delta = _nonZeroValues.Length > 64 ? 128 : 32; - } - } - - return delta; - } - #endregion - /// /// Sets all values to zero. /// public override void Clear() { - NonZerosCount = 0; - Array.Clear(_rowIndex, 0, _rowIndex.Length); + _storage.Clear(); } /// @@ -819,16 +690,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex throw DimensionsDontMatch(this, target, "target"); } + var targetStorage = sparseTarget.Storage; + // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value - sparseTarget._nonZeroValues = new Complex[NonZerosCount]; - sparseTarget._columnIndices = new int[NonZerosCount]; - sparseTarget.NonZerosCount = NonZerosCount; + targetStorage.ValueCount = _storage.ValueCount; + targetStorage.Values = new Complex[_storage.ValueCount]; + targetStorage.ColumnIndices = new int[_storage.ValueCount]; - if (NonZerosCount != 0) + if (_storage.ValueCount != 0) { - Array.Copy(_nonZeroValues, sparseTarget._nonZeroValues, NonZerosCount); - Buffer.BlockCopy(_columnIndices, 0, sparseTarget._columnIndices, 0, NonZerosCount * Constants.SizeOfInt); - Buffer.BlockCopy(_rowIndex, 0, sparseTarget._rowIndex, 0, RowCount * Constants.SizeOfInt); + Array.Copy(_storage.Values, targetStorage.Values, _storage.ValueCount); + Buffer.BlockCopy(_storage.ColumnIndices, 0, targetStorage.ColumnIndices, 0, _storage.ValueCount * Constants.SizeOfInt); + Buffer.BlockCopy(_storage.RowPointers, 0, targetStorage.RowPointers, 0, RowCount * Constants.SizeOfInt); } } @@ -840,14 +713,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// public override int GetHashCode() { - var hashNum = Math.Min(NonZerosCount, 25); + var values = _storage.Values; + var hashNum = Math.Min(_storage.ValueCount, 25); long hash = 0; for (var i = 0; i < hashNum; i++) { #if PORTABLE - hash ^= Precision.DoubleToInt64Bits(_nonZeroValues[i].Magnitude); + hash ^= Precision.DoubleToInt64Bits(values[i].Magnitude); #else - hash ^= BitConverter.DoubleToInt64Bits(_nonZeroValues[i].Magnitude); + hash ^= BitConverter.DoubleToInt64Bits(values[i].Magnitude); #endif } @@ -860,18 +734,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The transpose of this matrix. public override Matrix Transpose() { - var ret = new SparseMatrix(ColumnCount, RowCount) - { - _columnIndices = new int[NonZerosCount], - _nonZeroValues = new Complex[NonZerosCount] - }; + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var ret = new SparseMatrix(ColumnCount, RowCount); + var retStorage = ret.Storage; + retStorage.ColumnIndices = new int[valueCount]; + retStorage.Values = new Complex[valueCount]; // Do an 'inverse' CopyTo iterate over the rows - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : NonZerosCount; // Get the values for the current row if (startIndex == endIndex) @@ -882,7 +760,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var j = startIndex; j < endIndex; j++) { - ret.SetValueAt(_columnIndices[j], i, _nonZeroValues[j]); + retStorage.SetValueAt(columnIndices[j], i, values[j]); } } @@ -894,15 +772,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Complex FrobeniusNorm() { var transpose = (SparseMatrix)Transpose(); - var aat = this * transpose; + var aat = (this * transpose).Storage; - var norm = 0.0; + var norm = 0d; - for (var i = 0; i < aat._rowIndex.Length; i++) + for (var i = 0; i < aat.RowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = aat._rowIndex[i]; - var endIndex = i < aat._rowIndex.Length - 1 ? aat._rowIndex[i + 1] : aat.NonZerosCount; + var startIndex = aat.RowPointers[i]; + var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount; // Get the values for the current row if (startIndex == endIndex) @@ -913,9 +791,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var j = startIndex; j < endIndex; j++) { - if (i == aat._columnIndices[j]) + if (i == aat.ColumnIndices[j]) { - norm += aat._nonZeroValues[j].Magnitude; + norm += aat.Values[j].Magnitude; } } } @@ -928,12 +806,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The infinity norm of this matrix. public override Complex InfinityNorm() { - var norm = 0.0; - for (var i = 0; i < _rowIndex.Length; i++) + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var norm = 0d; + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -945,7 +827,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var s = 0.0; for (var j = startIndex; j < endIndex; j++) { - s += _nonZeroValues[j].Magnitude; + s += values[j].Magnitude; } norm = Math.Max(norm, s); @@ -1002,9 +884,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + // Determine bounds in columnIndices array where this item should be searched (using rowIndex) - var startIndex = _rowIndex[rowIndex]; - var endIndex = rowIndex < _rowIndex.Length - 1 ? _rowIndex[rowIndex + 1] : NonZerosCount; + var startIndex = rowPointers[rowIndex]; + var endIndex = rowIndex < rowPointers.Length - 1 ? rowPointers[rowIndex + 1] : valueCount; if (startIndex == endIndex) { @@ -1016,8 +902,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++) { // Copy code from At(row, column) to avoid unnecessary lock - var index = FindItem(rowIndex, i); - result[j] = index >= 0 ? _nonZeroValues[index] : 0.0; + var index = _storage.FindItem(rowIndex, i); + result[j] = index >= 0 ? values[index] : Complex.Zero; } } } @@ -1041,30 +927,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } else { + var resultStorage = resultSparseMatrix.Storage; + var lowerStorage = lowerSparseMatrix.Storage; + if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) { throw DimensionsDontMatch(this, resultSparseMatrix, lowerSparseMatrix); } - resultSparseMatrix.NonZerosCount = NonZerosCount + lowerSparseMatrix.NonZerosCount; - resultSparseMatrix._nonZeroValues = new Complex[resultSparseMatrix.NonZerosCount]; - resultSparseMatrix._columnIndices = new int[resultSparseMatrix.NonZerosCount]; + resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount; + resultStorage.Values = new Complex[resultStorage.ValueCount]; + resultStorage.ColumnIndices = new int[resultStorage.ValueCount]; - Array.Copy(_nonZeroValues, 0, resultSparseMatrix._nonZeroValues, 0, NonZerosCount); - Array.Copy(lowerSparseMatrix._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, NonZerosCount, lowerSparseMatrix.NonZerosCount); + Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount); + Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount); - Array.Copy(_columnIndices, 0, resultSparseMatrix._columnIndices, 0, NonZerosCount); - Array.Copy(_rowIndex, 0, resultSparseMatrix._rowIndex, 0, RowCount); + Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount); + Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount); // Copy and adjust lower column indices and rowIndex - for (int i = NonZerosCount, j = 0; i < resultSparseMatrix.NonZerosCount; i++, j++) + for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++) { - resultSparseMatrix._columnIndices[i] = lowerSparseMatrix._columnIndices[j] + ColumnCount; + resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount; } - for (int i = RowCount, j = 0; i < resultSparseMatrix.RowCount; i++, j++) + for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++) { - resultSparseMatrix._rowIndex[i] = lowerSparseMatrix._rowIndex[j] + NonZerosCount; + resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount; } } } @@ -1080,18 +969,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// public static SparseMatrix Identity(int order) { - var m = new SparseMatrix(order) - { - NonZerosCount = order, - _nonZeroValues = new Complex[order], - _columnIndices = new int[order] - }; + var m = new SparseMatrix(order); + var mStorage = m.Storage; + + mStorage.ValueCount = order; + mStorage.Values = new Complex[order]; + mStorage.ColumnIndices = new int[order]; for (var i = 0; i < order; i++) { - m._nonZeroValues[i] = 1.0; - m._columnIndices[i] = i; - m._rowIndex[i] = i; + mStorage.Values[i] = 1d; + mStorage.ColumnIndices[i] = i; + mStorage.RowPointers[i] = i; } return m; @@ -1126,21 +1015,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } var sparseMatrix = other as SparseMatrix; - if (sparseMatrix == null) { return base.Equals(other); } - if (NonZerosCount != sparseMatrix.NonZerosCount) + var otherStorage = sparseMatrix.Storage; + if (_storage.ValueCount != otherStorage.ValueCount) { return false; } // If all else fails, perform element wise comparison. - for (var index = 0; index < NonZerosCount; index++) + for (var index = 0; index < _storage.ValueCount; index++) { - if (!_nonZeroValues[index].AlmostEqual(sparseMatrix._nonZeroValues[index]) || _columnIndices[index] != sparseMatrix._columnIndices[index]) + if (!_storage.Values[index].AlmostEqual(otherStorage.Values[index]) || _storage.ColumnIndices[index] != otherStorage.ColumnIndices[index]) { return false; } @@ -1173,7 +1062,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex CopyTo(result); } - Control.LinearAlgebraProvider.ScaleArray(2.0, _nonZeroValues, _nonZeroValues); + Control.LinearAlgebraProvider.ScaleArray(2.0, _storage.Values, _storage.Values); return; } @@ -1193,16 +1082,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex left = sparseOther; } - for (var i = 0; i < left.RowCount; i++) + var leftStorage = left.Storage; + for (var i = 0; i < leftStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = left._rowIndex[i]; - var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount; + var startIndex = leftStorage.RowPointers[i]; + var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = left._columnIndices[j]; - var resVal = left._nonZeroValues[j] + result.At(i, columnIndex); + var columnIndex = leftStorage.ColumnIndices[j]; + var resVal = leftStorage.Values[j] + result.At(i, columnIndex); result.At(i, columnIndex, resVal); } } @@ -1219,6 +1109,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { var sparseOther = other as SparseMatrix; var sparseResult = result as SparseMatrix; + if (sparseOther == null || sparseResult == null) { base.DoSubtract(other, result); @@ -1231,18 +1122,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } + var otherStorage = sparseOther.Storage; + if (ReferenceEquals(this, sparseResult)) { - for (var i = 0; i < sparseOther.RowCount; i++) + for (var i = 0; i < otherStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = sparseOther._rowIndex[i]; - var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount; + var startIndex = otherStorage.RowPointers[i]; + var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = sparseOther._columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j]; + var columnIndex = otherStorage.ColumnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j]; result.At(i, columnIndex, resVal); } } @@ -1256,16 +1149,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex sparseResult.Negate(sparseResult); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = _columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j]; + var columnIndex = columnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) + values[j]; result.At(i, columnIndex, resVal); } } @@ -1296,10 +1194,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var row = 0; row < RowCount; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1308,8 +1210,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var index = start; index < end; index++) { - var column = _columnIndices[index]; - result.At(row, column, _nonZeroValues[index] * scalar); + var column = columnIndices[index]; + result.At(row, column, values[index] * scalar); } } } @@ -1320,7 +1222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._nonZeroValues[index] *= scalar); + CommonParallel.For(0, NonZerosCount, index => sparseResult.Storage.Values[index] *= scalar); } } @@ -1331,12 +1233,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The result of the multiplication. protected override void DoMultiply(Matrix other, Matrix result) { + result.Clear(); var columnVector = new DenseVector(other.RowCount); + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1350,7 +1259,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var sum = Complex.Zero; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * columnVector[_columnIndices[index]]; + sum += values[index] * columnVector[columnIndices[index]]; } result.At(row, column, sum); @@ -1365,11 +1274,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1378,7 +1292,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var sum = Complex.Zero; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * rightSide[_columnIndices[index]]; + sum += values[index] * rightSide[columnIndices[index]]; } result[row] = sum; @@ -1402,11 +1316,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } resultSparse.Clear(); + + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var otherStorage = otherSparse.Storage; + for (var j = 0; j < RowCount; j++) { // Get the begin / end index for the row - var startIndexOther = otherSparse._rowIndex[j]; - var endIndexOther = j < otherSparse._rowIndex.Length - 1 ? otherSparse._rowIndex[j + 1] : otherSparse.NonZerosCount; + var startIndexOther = otherStorage.RowPointers[j]; + var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount; if (startIndexOther == endIndexOther) { continue; @@ -1416,8 +1337,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { // Multiply row of matrix A on row of matrix B // Get the begin / end index for the row - var startIndexThis = _rowIndex[i]; - var endIndexThis = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndexThis = rowPointers[i]; + var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; if (startIndexThis == endIndexThis) { continue; @@ -1426,14 +1347,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var sum = Complex.Zero; for (var index = startIndexOther; index < endIndexOther; index++) { - var ind = FindItem(i, otherSparse._columnIndices[index]); + var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]); if (ind >= 0) { - sum += otherSparse._nonZeroValues[index] * _nonZeroValues[ind]; + sum += otherStorage.Values[index] * values[ind]; } } - resultSparse.SetValueAt(i, j, sum + result.At(i, j)); + resultSparse.Storage.SetValueAt(i, j, sum + result.At(i, j)); } } } @@ -1457,18 +1378,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] * other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] * other.At(i, columnIndices[j]); + if (!resVal.IsZero()) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1483,18 +1409,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] / other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] / other.At(i, columnIndices[j]); + if (!resVal.IsZero()) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1506,10 +1437,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The value at the current iteration along with its position (row, column, value). public override IEnumerable> IndexedEnumerator() { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1518,17 +1454,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var index = start; index < end; index++) { - yield return new Tuple(row, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(row, columnIndices[index], values[index]); } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < valueCount) { - for (var index = _rowIndex[lastRow]; index < NonZerosCount; index++) + for (var index = rowPointers[lastRow]; index < valueCount; index++) { - yield return new Tuple(lastRow, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(lastRow, columnIndices[index], values[index]); } } } @@ -1546,10 +1482,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } // todo: we might be able to speed this up by caching one half of the matrix + var rowPointers = _storage.RowPointers; for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1562,11 +1499,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < NonZerosCount) { - if (!CheckIfOppositesAreEqual(_rowIndex[lastRow], NonZerosCount, lastRow)) + if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow)) { return false; } @@ -1585,11 +1522,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If the values are equal or not. private bool CheckIfOppositesAreEqual(int start, int end, int row) { + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var index = start; index < end; index++) { - var column = _columnIndices[index]; + var column = columnIndices[index]; var opposite = At(column, row); - if (!_nonZeroValues[index].Equals(opposite)) + if (!values[index].Equals(opposite)) { return false; } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index a668ed54..76f09b49 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 using Generic; using Numerics; using Properties; + using Storage; using Threading; /// @@ -44,17 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 [Serializable] public class SparseMatrix : Matrix { - /// - /// 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 Complex32[] _nonZeroValues = new Complex32[0]; + readonly SparseCompressedRowMatrixStorage _storage; /// /// Gets the number of non zero elements in the matrix. @@ -62,16 +53,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The number of non zero elements. public int NonZerosCount { - get; - private set; + get { return _storage.ValueCount; } } - /// - /// 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. /// @@ -84,7 +68,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public SparseMatrix(int rows, int columns) : base(rows, columns) { - _rowIndex = new int[rows]; + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex32.Zero); } /// @@ -95,8 +79,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If is less than one. /// public SparseMatrix(int order) - : this(order, order) + : base(order, order) { + _storage = new SparseCompressedRowMatrixStorage(order, order, Complex32.Zero); } /// @@ -110,18 +95,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The value which we assign to each element of the matrix. public SparseMatrix(int rows, int columns, Complex32 value) - : this(rows, columns) + : base(rows, columns) { - if (value == 0.0f) + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex32.Zero); + + if (value.IsZero()) { return; } - NonZerosCount = rows * columns; - _nonZeroValues = new Complex32[NonZerosCount]; - _columnIndices = new int[NonZerosCount]; + var rowPointers = _storage.RowPointers; + var valueCount = _storage.ValueCount = rows * columns; + var columnIndices = _storage.ColumnIndices = new int[valueCount]; + var values = _storage.Values = new Complex32[valueCount]; - for (int i = 0, j = 0; i < _nonZeroValues.Length; i++, j++) + for (int i = 0, j = 0; i < values.Length; i++, j++) { // Reset column position to "0" if (j == columns) @@ -129,14 +117,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 j = 0; } - _nonZeroValues[i] = value; - _columnIndices[i] = j; + values[i] = value; + columnIndices[i] = j; } // Set proper row pointers - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { - _rowIndex[i] = ((i + 1) * columns) - columns; + rowPointers[i] = ((i + 1) * columns) - columns; } } @@ -149,8 +137,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If length is less than * . /// public SparseMatrix(int rows, int columns, Complex32[] array) - : this(rows, columns) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, Complex32.Zero); + if (rows * columns > array.Length) { throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); @@ -160,7 +150,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, array[i + (j * rows)]); + _storage.SetValueAt(i, j, array[i + (j * rows)]); } } } @@ -170,16 +160,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The 2D array to create this matrix from. public SparseMatrix(Complex32[,] array) - : this(array.GetLength(0), array.GetLength(1)) + : base(array.GetLength(0), array.GetLength(1)) { - var rows = array.GetLength(0); - var columns = array.GetLength(1); + _storage = new SparseCompressedRowMatrixStorage(array.GetLength(0), array.GetLength(1), Complex32.Zero); - for (var i = 0; i < rows; i++) + for (var i = 0; i < _storage.RowCount; i++) { - for (var j = 0; j < columns; j++) + for (var j = 0; j < _storage.ColumnCount; j++) { - SetValueAt(i, j, array[i, j]); + _storage.SetValueAt(i, j, array[i, j]); } } } @@ -190,8 +179,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// /// The matrix to copy. public SparseMatrix(Matrix matrix) - : this(matrix.RowCount, matrix.ColumnCount) + : base(matrix.RowCount, matrix.ColumnCount) { + _storage = new SparseCompressedRowMatrixStorage(matrix.RowCount, matrix.ColumnCount, Complex32.Zero); + var sparseMatrix = matrix as SparseMatrix; var rows = matrix.RowCount; @@ -203,23 +194,28 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, matrix.At(i, j)); + _storage.SetValueAt(i, j, matrix.At(i, j)); } } } else { - NonZerosCount = sparseMatrix.NonZerosCount; - _rowIndex = new int[rows]; - _columnIndices = new int[NonZerosCount]; - _nonZeroValues = new Complex32[NonZerosCount]; + var matrixStorage = sparseMatrix.Storage; + var valueCount = _storage.ValueCount = matrixStorage.ValueCount; + _storage.ColumnIndices = new int[valueCount]; + _storage.Values = new Complex32[valueCount]; - Array.Copy(sparseMatrix._nonZeroValues, _nonZeroValues, NonZerosCount); - Array.Copy(sparseMatrix._columnIndices, _columnIndices, NonZerosCount); - Array.Copy(sparseMatrix._rowIndex, _rowIndex, rows); + Array.Copy(matrixStorage.Values, _storage.Values, valueCount); + Array.Copy(matrixStorage.ColumnIndices, _storage.ColumnIndices, valueCount); + Array.Copy(matrixStorage.RowPointers, _storage.RowPointers, rows); } } + internal SparseCompressedRowMatrixStorage Storage + { + get { return _storage; } + } + /// /// Creates a SparseMatrix for the given number of rows and columns. /// @@ -297,15 +293,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Where to store the lower triangle. private void LowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row >= _columnIndices[j]) + if (row >= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -359,15 +360,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Where to store the lower triangle. private void UpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row <= _columnIndices[j]) + if (row <= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -426,18 +432,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var result = (SparseMatrix)CreateMatrix(rowCount, columnCount); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (int i = rowIndex, row = 0; i < rowMax; i++, row++) { - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (int j = startIndex; j < endIndex; j++) { // check if the column index is in the range - if ((_columnIndices[j] >= columnIndex) && (_columnIndices[j] < columnIndex + columnCount)) + if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount)) { - var column = _columnIndices[j] - columnIndex; - result.SetValueAt(row, column, _nonZeroValues[j]); + var column = columnIndices[j] - columnIndex; + result._storage.SetValueAt(row, column, values[j]); } } } @@ -494,15 +505,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Where to store the lower triangle. private void StrictlyLowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row > _columnIndices[j]) + if (row > columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -557,15 +573,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Where to store the lower triangle. private void StrictlyUpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row < _columnIndices[j]) + if (row < columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -582,13 +603,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// An array containing the matrix's elements. public override Complex32[] ToColumnWiseArray() { + var values = _storage.Values; var ret = new Complex32[RowCount * ColumnCount]; for (var j = 0; j < ColumnCount; j++) { for (var i = 0; i < RowCount; i++) { - var index = FindItem(i, j); - ret[(j * RowCount) + i] = index >= 0 ? _nonZeroValues[index] : 0.0f; + var index = _storage.FindItem(i, j); + ret[(j * RowCount) + i] = index >= 0 ? values[index] : 0.0f; } } @@ -609,8 +631,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// public override Complex32 At(int row, int column) { - var index = FindItem(row, column); - return index >= 0 ? _nonZeroValues[index] : 0.0f; + return _storage.GetValueAt(row, column); } /// @@ -627,166 +648,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// public override void At(int row, int column, Complex32 value) { - 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, Complex32 value) - { - var index = FindItem(row, column); - if (index >= 0) - { - // Non-zero item found in matrix - if (value == 0.0f) - { - // Delete existing item - DeleteItemByIndex(index, row); - } - else - { - // Update item - _nonZeroValues[index] = value; - } - } - else - { - // Item not found. Add new value - if (value == 0.0f) - { - return; - } - - index = ~index; - - // Check if the storage needs to be increased - if ((NonZerosCount == _nonZeroValues.Length) && (NonZerosCount < ((long)RowCount * 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(_nonZeroValues.Length + GrowthSize(), (long)RowCount * ColumnCount); - if (size > int.MaxValue) - { - throw new NotSupportedException(Resources.TooManyElements); - } - - Array.Resize(ref _nonZeroValues, (int)size); - Array.Resize(ref _columnIndices, (int)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 = NonZerosCount - 1; i > index - 1; i--) - { - _nonZeroValues[i + 1] = _nonZeroValues[i]; - _columnIndices[i + 1] = _columnIndices[i]; - } - - // Add the value and the column index - _nonZeroValues[index] = value; - _columnIndices[index] = column; - - // increase the number of non-zero numbers by one - 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 < _rowIndex.Length; i++) - { - _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 < NonZerosCount; i++) - { - _nonZeroValues[i - 1] = _nonZeroValues[i]; - _columnIndices[i - 1] = _columnIndices[i]; - } - - // Decrease value in Row - for (var i = row + 1; i < _rowIndex.Length; i++) - { - _rowIndex[i] -= 1; - } - - 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 ((NonZerosCount > 1024) && (NonZerosCount < _nonZeroValues.Length / 2)) - { - Array.Resize(ref _nonZeroValues, NonZerosCount); - Array.Resize(ref _columnIndices, 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 = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; - return Array.BinarySearch(_columnIndices, startIndex, endIndex - startIndex, column); - } - - /// - /// Calculates the amount with which to grow the storage array's if they need to be - /// increased in size. - /// - /// The amount grown. - private int GrowthSize() - { - int delta; - if (_nonZeroValues.Length > 1024) - { - delta = _nonZeroValues.Length / 4; - } - else - { - if (_nonZeroValues.Length > 256) - { - delta = 512; - } - else - { - delta = _nonZeroValues.Length > 64 ? 128 : 32; - } - } - - return delta; + _storage.SetValueAt(row, column, value); } - #endregion /// /// Sets all values to zero. /// public override void Clear() { - NonZerosCount = 0; - Array.Clear(_rowIndex, 0, _rowIndex.Length); + _storage.Clear(); } /// @@ -820,16 +690,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw DimensionsDontMatch(this, target, "target"); } + + var targetStorage = sparseTarget.Storage; + // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value - sparseTarget._nonZeroValues = new Complex32[NonZerosCount]; - sparseTarget._columnIndices = new int[NonZerosCount]; - sparseTarget.NonZerosCount = NonZerosCount; + targetStorage.ValueCount = _storage.ValueCount; + targetStorage.Values = new Complex32[_storage.ValueCount]; + targetStorage.ColumnIndices = new int[_storage.ValueCount]; - if (NonZerosCount != 0) + if (_storage.ValueCount != 0) { - Array.Copy(_nonZeroValues, sparseTarget._nonZeroValues, NonZerosCount); - Buffer.BlockCopy(_columnIndices, 0, sparseTarget._columnIndices, 0, NonZerosCount * Constants.SizeOfInt); - Buffer.BlockCopy(_rowIndex, 0, sparseTarget._rowIndex, 0, RowCount * Constants.SizeOfInt); + Array.Copy(_storage.Values, targetStorage.Values, _storage.ValueCount); + Buffer.BlockCopy(_storage.ColumnIndices, 0, targetStorage.ColumnIndices, 0, _storage.ValueCount * Constants.SizeOfInt); + Buffer.BlockCopy(_storage.RowPointers, 0, targetStorage.RowPointers, 0, RowCount * Constants.SizeOfInt); } } @@ -841,14 +714,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// public override int GetHashCode() { - var hashNum = Math.Min(NonZerosCount, 25); + var values = _storage.Values; + var hashNum = Math.Min(_storage.ValueCount, 25); long hash = 0; for (var i = 0; i < hashNum; i++) { #if PORTABLE - hash ^= Precision.DoubleToInt64Bits(_nonZeroValues[i].Magnitude); + hash ^= Precision.DoubleToInt64Bits(values[i].Magnitude); #else - hash ^= BitConverter.DoubleToInt64Bits(_nonZeroValues[i].Magnitude); + hash ^= BitConverter.DoubleToInt64Bits(values[i].Magnitude); #endif } @@ -861,18 +735,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The transpose of this matrix. public override Matrix Transpose() { - var ret = new SparseMatrix(ColumnCount, RowCount) - { - _columnIndices = new int[NonZerosCount], - _nonZeroValues = new Complex32[NonZerosCount] - }; + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var ret = new SparseMatrix(ColumnCount, RowCount); + var retStorage = ret.Storage; + retStorage.ColumnIndices = new int[valueCount]; + retStorage.Values = new Complex32[valueCount]; // Do an 'inverse' CopyTo iterate over the rows - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : NonZerosCount; // Get the values for the current row if (startIndex == endIndex) @@ -883,7 +761,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var j = startIndex; j < endIndex; j++) { - ret.SetValueAt(_columnIndices[j], i, _nonZeroValues[j]); + retStorage.SetValueAt(columnIndices[j], i, values[j]); } } @@ -895,15 +773,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Complex32 FrobeniusNorm() { var transpose = (SparseMatrix)Transpose(); - var aat = this * transpose; + var aat = (this * transpose).Storage; - var norm = 0.0f; + var norm = 0f; - for (var i = 0; i < aat._rowIndex.Length; i++) + for (var i = 0; i < aat.RowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = aat._rowIndex[i]; - var endIndex = i < aat._rowIndex.Length - 1 ? aat._rowIndex[i + 1] : aat.NonZerosCount; + var startIndex = aat.RowPointers[i]; + var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount; // Get the values for the current row if (startIndex == endIndex) @@ -914,9 +792,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var j = startIndex; j < endIndex; j++) { - if (i == aat._columnIndices[j]) + if (i == aat.ColumnIndices[j]) { - norm += aat._nonZeroValues[j].Magnitude; + norm += aat.Values[j].Magnitude; } } } @@ -929,12 +807,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The infinity norm of this matrix. public override Complex32 InfinityNorm() { - var norm = 0.0f; - for (var i = 0; i < _rowIndex.Length; i++) + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var norm = 0f; + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -946,7 +828,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var s = 0.0f; for (var j = startIndex; j < endIndex; j++) { - s += _nonZeroValues[j].Magnitude; + s += values[j].Magnitude; } norm = Math.Max(norm, s); @@ -1003,9 +885,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + // Determine bounds in columnIndices array where this item should be searched (using rowIndex) - var startIndex = _rowIndex[rowIndex]; - var endIndex = rowIndex < _rowIndex.Length - 1 ? _rowIndex[rowIndex + 1] : NonZerosCount; + var startIndex = rowPointers[rowIndex]; + var endIndex = rowIndex < rowPointers.Length - 1 ? rowPointers[rowIndex + 1] : valueCount; if (startIndex == endIndex) { @@ -1017,8 +903,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++) { // Copy code from At(row, column) to avoid unnecessary lock - var index = FindItem(rowIndex, i); - result[j] = index >= 0 ? _nonZeroValues[index] : 0.0f; + var index = _storage.FindItem(rowIndex, i); + result[j] = index >= 0 ? values[index] : Complex32.Zero; } } } @@ -1042,30 +928,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } else { + var resultStorage = resultSparseMatrix.Storage; + var lowerStorage = lowerSparseMatrix.Storage; + if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) { - throw DimensionsDontMatch(this, lowerSparseMatrix, resultSparseMatrix); + throw DimensionsDontMatch(this, resultSparseMatrix, lowerSparseMatrix); } - resultSparseMatrix.NonZerosCount = NonZerosCount + lowerSparseMatrix.NonZerosCount; - resultSparseMatrix._nonZeroValues = new Complex32[resultSparseMatrix.NonZerosCount]; - resultSparseMatrix._columnIndices = new int[resultSparseMatrix.NonZerosCount]; + resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount; + resultStorage.Values = new Complex32[resultStorage.ValueCount]; + resultStorage.ColumnIndices = new int[resultStorage.ValueCount]; - Array.Copy(_nonZeroValues, 0, resultSparseMatrix._nonZeroValues, 0, NonZerosCount); - Array.Copy(lowerSparseMatrix._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, NonZerosCount, lowerSparseMatrix.NonZerosCount); + Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount); + Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount); - Array.Copy(_columnIndices, 0, resultSparseMatrix._columnIndices, 0, NonZerosCount); - Array.Copy(_rowIndex, 0, resultSparseMatrix._rowIndex, 0, RowCount); + Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount); + Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount); // Copy and adjust lower column indices and rowIndex - for (int i = NonZerosCount, j = 0; i < resultSparseMatrix.NonZerosCount; i++, j++) + for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++) { - resultSparseMatrix._columnIndices[i] = lowerSparseMatrix._columnIndices[j] + ColumnCount; + resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount; } - for (int i = RowCount, j = 0; i < resultSparseMatrix.RowCount; i++, j++) + for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++) { - resultSparseMatrix._rowIndex[i] = lowerSparseMatrix._rowIndex[j] + NonZerosCount; + resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount; } } } @@ -1081,18 +970,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// public static SparseMatrix Identity(int order) { - var m = new SparseMatrix(order) - { - NonZerosCount = order, - _nonZeroValues = new Complex32[order], - _columnIndices = new int[order] - }; + var m = new SparseMatrix(order); + var mStorage = m.Storage; + + mStorage.ValueCount = order; + mStorage.Values = new Complex32[order]; + mStorage.ColumnIndices = new int[order]; for (var i = 0; i < order; i++) { - m._nonZeroValues[i] = 1.0f; - m._columnIndices[i] = i; - m._rowIndex[i] = i; + mStorage.Values[i] = 1f; + mStorage.ColumnIndices[i] = i; + mStorage.RowPointers[i] = i; } return m; @@ -1127,21 +1016,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } var sparseMatrix = other as SparseMatrix; - if (sparseMatrix == null) { return base.Equals(other); } - if (NonZerosCount != sparseMatrix.NonZerosCount) + var otherStorage = sparseMatrix.Storage; + if (_storage.ValueCount != otherStorage.ValueCount) { return false; } // If all else fails, perform element wise comparison. - for (var index = 0; index < NonZerosCount; index++) + for (var index = 0; index < _storage.ValueCount; index++) { - if (!_nonZeroValues[index].AlmostEqual(sparseMatrix._nonZeroValues[index]) || _columnIndices[index] != sparseMatrix._columnIndices[index]) + if (!_storage.Values[index].AlmostEqual(otherStorage.Values[index]) || _storage.ColumnIndices[index] != otherStorage.ColumnIndices[index]) { return false; } @@ -1174,7 +1063,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 CopyTo(result); } - Control.LinearAlgebraProvider.ScaleArray(2.0f, _nonZeroValues, _nonZeroValues); + Control.LinearAlgebraProvider.ScaleArray(2.0f, _storage.Values, _storage.Values); return; } @@ -1194,16 +1083,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 left = sparseOther; } - for (var i = 0; i < left.RowCount; i++) + var leftStorage = left.Storage; + for (var i = 0; i < leftStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = left._rowIndex[i]; - var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount; + var startIndex = leftStorage.RowPointers[i]; + var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = left._columnIndices[j]; - var resVal = left._nonZeroValues[j] + result.At(i, columnIndex); + var columnIndex = leftStorage.ColumnIndices[j]; + var resVal = leftStorage.Values[j] + result.At(i, columnIndex); result.At(i, columnIndex, resVal); } } @@ -1232,18 +1122,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } + var otherStorage = sparseOther.Storage; + if (ReferenceEquals(this, sparseResult)) { - for (var i = 0; i < sparseOther.RowCount; i++) + for (var i = 0; i < otherStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = sparseOther._rowIndex[i]; - var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount; + var startIndex = otherStorage.RowPointers[i]; + var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = sparseOther._columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j]; + var columnIndex = otherStorage.ColumnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j]; result.At(i, columnIndex, resVal); } } @@ -1257,16 +1149,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 sparseResult.Negate(sparseResult); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = _columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j]; + var columnIndex = columnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) + values[j]; result.At(i, columnIndex, resVal); } } @@ -1297,10 +1194,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var row = 0; row < RowCount; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1309,8 +1210,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var index = start; index < end; index++) { - var column = _columnIndices[index]; - result.At(row, column, _nonZeroValues[index] * scalar); + var column = columnIndices[index]; + result.At(row, column, values[index] * scalar); } } } @@ -1321,7 +1222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._nonZeroValues[index] *= scalar); + CommonParallel.For(0, NonZerosCount, index => sparseResult.Storage.Values[index] *= scalar); } } @@ -1334,11 +1235,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { result.Clear(); var columnVector = new DenseVector(other.RowCount); + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1352,7 +1259,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var sum = Complex32.Zero; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * columnVector[_columnIndices[index]]; + sum += values[index] * columnVector[columnIndices[index]]; } result.At(row, column, sum); @@ -1367,11 +1274,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1380,7 +1292,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var sum = Complex32.Zero; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * rightSide[_columnIndices[index]]; + sum += values[index] * rightSide[columnIndices[index]]; } result[row] = sum; @@ -1404,11 +1316,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } resultSparse.Clear(); + + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var otherStorage = otherSparse.Storage; + for (var j = 0; j < RowCount; j++) { // Get the begin / end index for the row - var startIndexOther = otherSparse._rowIndex[j]; - var endIndexOther = j < otherSparse._rowIndex.Length - 1 ? otherSparse._rowIndex[j + 1] : otherSparse.NonZerosCount; + var startIndexOther = otherStorage.RowPointers[j]; + var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount; if (startIndexOther == endIndexOther) { continue; @@ -1418,8 +1337,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { // Multiply row of matrix A on row of matrix B // Get the begin / end index for the row - var startIndexThis = _rowIndex[i]; - var endIndexThis = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndexThis = rowPointers[i]; + var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; if (startIndexThis == endIndexThis) { continue; @@ -1428,14 +1347,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var sum = Complex32.Zero; for (var index = startIndexOther; index < endIndexOther; index++) { - var ind = FindItem(i, otherSparse._columnIndices[index]); + var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]); if (ind >= 0) { - sum += otherSparse._nonZeroValues[index] * _nonZeroValues[ind]; + sum += otherStorage.Values[index] * values[ind]; } } - resultSparse.SetValueAt(i, j, sum + result.At(i, j)); + resultSparse.Storage.SetValueAt(i, j, sum + result.At(i, j)); } } } @@ -1459,18 +1378,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] * other.At(i, _columnIndices[j]); - if (resVal != 0.0f) + var resVal = values[j] * other.At(i, columnIndices[j]); + if (!resVal.IsZero()) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1485,18 +1409,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] / other.At(i, _columnIndices[j]); - if (resVal != 0.0f) + var resVal = values[j] / other.At(i, columnIndices[j]); + if (!resVal.IsZero()) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1508,10 +1437,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The value at the current iteration along with its position (row, column, value). public override IEnumerable> IndexedEnumerator() { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1520,17 +1454,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var index = start; index < end; index++) { - yield return new Tuple(row, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(row, columnIndices[index], values[index]); } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < valueCount) { - for (var index = _rowIndex[lastRow]; index < NonZerosCount; index++) + for (var index = rowPointers[lastRow]; index < valueCount; index++) { - yield return new Tuple(lastRow, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(lastRow, columnIndices[index], values[index]); } } } @@ -1548,10 +1482,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } // todo: we might be able to speed this up by caching one half of the matrix + var rowPointers = _storage.RowPointers; for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1564,11 +1499,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < NonZerosCount) { - if (!CheckIfOppositesAreEqual(_rowIndex[lastRow], NonZerosCount, lastRow)) + if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow)) { return false; } @@ -1587,11 +1522,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If the values are equal or not. private bool CheckIfOppositesAreEqual(int start, int end, int row) { + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var index = start; index < end; index++) { - var column = _columnIndices[index]; + var column = columnIndices[index]; var opposite = At(column, row); - if (!_nonZeroValues[index].Equals(opposite)) + if (!values[index].Equals(opposite)) { return false; } diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 89d26f2a..6644e608 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System.Collections.Generic; using Generic; using Properties; + using Storage; using Threading; /// @@ -43,17 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double [Serializable] public class SparseMatrix : Matrix { - /// - /// 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]; + readonly SparseCompressedRowMatrixStorage _storage; /// /// Gets the number of non zero elements in the matrix. @@ -61,15 +52,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The number of non zero elements. public int NonZerosCount { - get; - private set; + get { return _storage.ValueCount; } } - - /// - /// 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. @@ -80,9 +64,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// The number of columns. /// - public SparseMatrix(int rows, int columns) : base(rows, columns) + public SparseMatrix(int rows, int columns) + : base(rows, columns) { - _rowIndex = new int[rows]; + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0d); } /// @@ -92,8 +77,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// /// If is less than one. /// - public SparseMatrix(int order) : this(order, order) + public SparseMatrix(int order) + : base(order, order) { + _storage = new SparseCompressedRowMatrixStorage(order, order, 0d); } /// @@ -106,18 +93,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// 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) + public SparseMatrix(int rows, int columns, double value) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0d); + if (value == 0.0) { return; } - NonZerosCount = rows * columns; - _nonZeroValues = new double[NonZerosCount]; - _columnIndices = new int[NonZerosCount]; + var rowPointers = _storage.RowPointers; + var valueCount = _storage.ValueCount = rows * columns; + var columnIndices = _storage.ColumnIndices = new int[valueCount]; + var values = _storage.Values = new double[valueCount]; - for (int i = 0, j = 0; i < _nonZeroValues.Length; i++, j++) + for (int i = 0, j = 0; i < values.Length; i++, j++) { // Reset column position to "0" if (j == columns) @@ -125,14 +116,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double j = 0; } - _nonZeroValues[i] = value; - _columnIndices[i] = j; + values[i] = value; + columnIndices[i] = j; } // Set proper row pointers - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { - _rowIndex[i] = ((i + 1) * columns) - columns; + rowPointers[i] = ((i + 1) * columns) - columns; } } @@ -144,8 +135,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The one dimensional array to create this matrix from. This array should store the matrix in column-major order. see: http://en.wikipedia.org/wiki/Column-major_order /// If length is less than * . /// - public SparseMatrix(int rows, int columns, double[] array) : this(rows, columns) + public SparseMatrix(int rows, int columns, double[] array) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0d); + if (rows * columns > array.Length) { throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); @@ -155,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, array[i + (j * rows)]); + _storage.SetValueAt(i, j, array[i + (j * rows)]); } } } @@ -164,16 +158,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// 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)) + public SparseMatrix(double[,] array) + : base(array.GetLength(0), array.GetLength(1)) { - var rows = array.GetLength(0); - var columns = array.GetLength(1); + _storage = new SparseCompressedRowMatrixStorage(array.GetLength(0), array.GetLength(1), 0d); - for (var i = 0; i < rows; i++) + for (var i = 0; i < _storage.RowCount; i++) { - for (var j = 0; j < columns; j++) + for (var j = 0; j < _storage.ColumnCount; j++) { - SetValueAt(i, j, array[i, j]); + _storage.SetValueAt(i, j, array[i, j]); } } } @@ -183,8 +177,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// the values from the given matrix. /// /// The matrix to copy. - public SparseMatrix(Matrix matrix) : this(matrix.RowCount, matrix.ColumnCount) + public SparseMatrix(Matrix matrix) + : base(matrix.RowCount, matrix.ColumnCount) { + _storage = new SparseCompressedRowMatrixStorage(matrix.RowCount, matrix.ColumnCount, 0d); + var sparseMatrix = matrix as SparseMatrix; var rows = matrix.RowCount; @@ -196,23 +193,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, matrix.At(i, j)); + _storage.SetValueAt(i, j, matrix.At(i, j)); } } } else { - NonZerosCount = sparseMatrix.NonZerosCount; - _rowIndex = new int[rows]; - _columnIndices = new int[NonZerosCount]; - _nonZeroValues = new double[NonZerosCount]; + var matrixStorage = sparseMatrix.Storage; + var valueCount = _storage.ValueCount = matrixStorage.ValueCount; + _storage.ColumnIndices = new int[valueCount]; + _storage.Values = new double[valueCount]; - Buffer.BlockCopy(sparseMatrix._nonZeroValues, 0, _nonZeroValues, 0, NonZerosCount * Constants.SizeOfDouble); - Buffer.BlockCopy(sparseMatrix._columnIndices, 0, _columnIndices, 0, NonZerosCount * Constants.SizeOfInt); - Buffer.BlockCopy(sparseMatrix._rowIndex, 0, _rowIndex, 0, rows * Constants.SizeOfInt); + Buffer.BlockCopy(matrixStorage.Values, 0, _storage.Values, 0, valueCount * Constants.SizeOfDouble); + Buffer.BlockCopy(matrixStorage.ColumnIndices, 0, _storage.ColumnIndices, 0, valueCount * Constants.SizeOfInt); + Buffer.BlockCopy(matrixStorage.RowPointers, 0, _storage.RowPointers, 0, rows * Constants.SizeOfInt); } } + internal SparseCompressedRowMatrixStorage Storage + { + get { return _storage; } + } + /// /// Creates a SparseMatrix for the given number of rows and columns. /// @@ -290,15 +292,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Where to store the lower triangle. private void LowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row >= _columnIndices[j]) + if (row >= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -352,15 +359,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Where to store the lower triangle. private void UpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row <= _columnIndices[j]) + if (row <= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -419,18 +431,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double var result = (SparseMatrix)CreateMatrix(rowCount, columnCount); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (int i = rowIndex, row = 0; i < rowMax; i++, row++) { - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (int j = startIndex; j < endIndex; j++) { // check if the column index is in the range - if ((_columnIndices[j] >= columnIndex) && (_columnIndices[j] < columnIndex + columnCount)) + if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount)) { - var column = _columnIndices[j] - columnIndex; - result.SetValueAt(row, column, _nonZeroValues[j]); + var column = columnIndices[j] - columnIndex; + result._storage.SetValueAt(row, column, values[j]); } } } @@ -487,15 +504,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Where to store the lower triangle. private void StrictlyLowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row > _columnIndices[j]) + if (row > columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -550,15 +572,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Where to store the lower triangle. private void StrictlyUpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row < _columnIndices[j]) + if (row < columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -575,13 +602,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// An array containing the matrix's elements. public override double[] ToColumnWiseArray() { + var values = _storage.Values; var ret = new double[RowCount * ColumnCount]; for (var j = 0; j < ColumnCount; j++) { for (var i = 0; i < RowCount; i++) { - var index = FindItem(i, j); - ret[(j * RowCount) + i] = index >= 0 ? _nonZeroValues[index] : 0.0; + var index = _storage.FindItem(i, j); + ret[(j * RowCount) + i] = index >= 0 ? values[index] : 0.0; } } @@ -602,7 +630,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override double At(int row, int column) { - return GetValueAt(row, column); + return _storage.GetValueAt(row, column); } /// @@ -619,185 +647,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override void At(int row, int column, double value) { - SetValueAt(row, column, value); - } - - #region Internal methods - CRS storage implementation - - /// - /// Retrieves the requested element without range checking. - /// - /// - /// The row of the element. - /// - /// - /// The column of the element. - /// - /// - /// The requested element. - /// - private double GetValueAt(int row, int column) - { - var index = FindItem(row, column); - return index >= 0 ? _nonZeroValues[index] : 0.0; - } - - /// - /// 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 = FindItem(row, column); - if (index >= 0) - { - // Non-zero item found in matrix - if (value == 0.0) - { - // Delete existing item - DeleteItemByIndex(index, row); - } - else - { - // Update item - _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 ((NonZerosCount == _nonZeroValues.Length) && (NonZerosCount < ((long)RowCount * 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(_nonZeroValues.Length + GrowthSize(), (long)RowCount * ColumnCount); - if (size > int.MaxValue) - { - throw new NotSupportedException(Resources.TooManyElements); - } - - Array.Resize(ref _nonZeroValues, (int)size); - Array.Resize(ref _columnIndices, (int)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 = NonZerosCount - 1; i > index - 1; i--) - { - _nonZeroValues[i + 1] = _nonZeroValues[i]; - _columnIndices[i + 1] = _columnIndices[i]; - } - - // Add the value and the column index - _nonZeroValues[index] = value; - _columnIndices[index] = column; - - // increase the number of non-zero numbers by one - 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 < _rowIndex.Length; i++) - { - _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 < NonZerosCount; i++) - { - _nonZeroValues[i - 1] = _nonZeroValues[i]; - _columnIndices[i - 1] = _columnIndices[i]; - } - - // Decrease value in Row - for (var i = row + 1; i < _rowIndex.Length; i++) - { - _rowIndex[i] -= 1; - } - - 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 ((NonZerosCount > 1024) && (NonZerosCount < _nonZeroValues.Length / 2)) - { - Array.Resize(ref _nonZeroValues, NonZerosCount); - Array.Resize(ref _columnIndices, NonZerosCount); - } + _storage.SetValueAt(row, column, value); } - - /// - /// 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 = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; - return Array.BinarySearch(_columnIndices, startIndex, endIndex - startIndex, column); - } - - /// - /// Calculates the amount with which to grow the storage array's if they need to be - /// increased in size. - /// - /// The amount grown. - private int GrowthSize() - { - int delta; - if (_nonZeroValues.Length > 1024) - { - delta = _nonZeroValues.Length / 4; - } - else - { - if (_nonZeroValues.Length > 256) - { - delta = 512; - } - else - { - delta = _nonZeroValues.Length > 64 ? 128 : 32; - } - } - - return delta; - } - #endregion /// /// Sets all values to zero. /// public override void Clear() { - NonZerosCount = 0; - Array.Clear(_rowIndex, 0, _rowIndex.Length); + _storage.Clear(); } /// @@ -831,16 +689,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw DimensionsDontMatch(this, target, "target"); } + var targetStorage = sparseTarget.Storage; + // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value - sparseTarget._nonZeroValues = new double[NonZerosCount]; - sparseTarget._columnIndices = new int[NonZerosCount]; - sparseTarget.NonZerosCount = NonZerosCount; + targetStorage.ValueCount = _storage.ValueCount; + targetStorage.Values = new double[_storage.ValueCount]; + targetStorage.ColumnIndices = new int[_storage.ValueCount]; - if (NonZerosCount != 0) + if (_storage.ValueCount != 0) { - Buffer.BlockCopy(_nonZeroValues, 0, sparseTarget._nonZeroValues, 0, NonZerosCount*Constants.SizeOfDouble); - Buffer.BlockCopy(_columnIndices, 0, sparseTarget._columnIndices, 0, NonZerosCount*Constants.SizeOfInt); - Buffer.BlockCopy(_rowIndex, 0, sparseTarget._rowIndex, 0, RowCount*Constants.SizeOfInt); + Buffer.BlockCopy(_storage.Values, 0, targetStorage.Values, 0, _storage.ValueCount * Constants.SizeOfDouble); + Buffer.BlockCopy(_storage.ColumnIndices, 0, targetStorage.ColumnIndices, 0, _storage.ValueCount * Constants.SizeOfInt); + Buffer.BlockCopy(_storage.RowPointers, 0, targetStorage.RowPointers, 0, RowCount * Constants.SizeOfInt); } } @@ -852,14 +712,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override int GetHashCode() { - var hashNum = Math.Min(NonZerosCount, 25); + var values = _storage.Values; + var hashNum = Math.Min(_storage.ValueCount, 25); long hash = 0; for (var i = 0; i < hashNum; i++) { #if PORTABLE - hash ^= Precision.DoubleToInt64Bits(_nonZeroValues[i]); + hash ^= Precision.DoubleToInt64Bits(values[i]); #else - hash ^= BitConverter.DoubleToInt64Bits(_nonZeroValues[i]); + hash ^= BitConverter.DoubleToInt64Bits(values[i]); #endif } @@ -872,18 +733,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The transpose of this matrix. public override Matrix Transpose() { - var ret = new SparseMatrix(ColumnCount, RowCount) - { - _columnIndices = new int[NonZerosCount], - _nonZeroValues = new double[NonZerosCount] - }; + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var ret = new SparseMatrix(ColumnCount, RowCount); + var retStorage = ret.Storage; + retStorage.ColumnIndices = new int[valueCount]; + retStorage.Values = new double[valueCount]; // Do an 'inverse' CopyTo iterate over the rows - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -894,7 +759,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var j = startIndex; j < endIndex; j++) { - ret.SetValueAt(_columnIndices[j], i, _nonZeroValues[j]); + retStorage.SetValueAt(columnIndices[j], i, values[j]); } } @@ -906,15 +771,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override double FrobeniusNorm() { var transpose = (SparseMatrix)Transpose(); - var aat = this * transpose; + var aat = (this * transpose).Storage; - var norm = 0.0; + var norm = 0d; - for (var i = 0; i < aat._rowIndex.Length; i++) + for (var i = 0; i < aat.RowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = aat._rowIndex[i]; - var endIndex = i < aat._rowIndex.Length - 1 ? aat._rowIndex[i + 1] : aat.NonZerosCount; + var startIndex = aat.RowPointers[i]; + var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount; // Get the values for the current row if (startIndex == endIndex) @@ -925,27 +790,30 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var j = startIndex; j < endIndex; j++) { - if (i == aat._columnIndices[j]) + if (i == aat.ColumnIndices[j]) { - norm += Math.Abs(aat._nonZeroValues[j]); + norm += Math.Abs(aat.Values[j]); } } } - norm = Math.Sqrt(norm); - return norm; + return Math.Sqrt(norm); } /// Calculates the infinity norm of this matrix. /// The infinity norm of this matrix. public override double InfinityNorm() { - var norm = 0.0; - for (var i = 0; i < _rowIndex.Length; i++) + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var norm = 0d; + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -954,10 +822,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double continue; } - var s = 0.0; + var s = 0d; for (var j = startIndex; j < endIndex; j++) { - s += Math.Abs(_nonZeroValues[j]); + s += Math.Abs(values[j]); } norm = Math.Max(norm, s); @@ -1014,9 +882,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + // Determine bounds in columnIndices array where this item should be searched (using rowIndex) - var startIndex = _rowIndex[rowIndex]; - var endIndex = rowIndex < _rowIndex.Length - 1 ? _rowIndex[rowIndex + 1] : NonZerosCount; + var startIndex = rowPointers[rowIndex]; + var endIndex = rowIndex < rowPointers.Length - 1 ? rowPointers[rowIndex + 1] : valueCount; if (startIndex == endIndex) { @@ -1028,8 +900,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++) { // Copy code from At(row, column) to avoid unnecessary lock - var index = FindItem(rowIndex, i); - result[j] = index >= 0 ? _nonZeroValues[index] : 0.0; + var index = _storage.FindItem(rowIndex, i); + result[j] = index >= 0 ? values[index] : 0d; } } } @@ -1053,30 +925,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else { - if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) + var resultStorage = resultSparseMatrix.Storage; + var lowerStorage = lowerSparseMatrix.Storage; + + if (resultStorage.RowCount != RowCount + lowerStorage.RowCount || resultStorage.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) { throw DimensionsDontMatch(this, lowerSparseMatrix, resultSparseMatrix); } - resultSparseMatrix.NonZerosCount = NonZerosCount + lowerSparseMatrix.NonZerosCount; - resultSparseMatrix._nonZeroValues = new double[resultSparseMatrix.NonZerosCount]; - resultSparseMatrix._columnIndices = new int[resultSparseMatrix.NonZerosCount]; - - Array.Copy(_nonZeroValues, 0, resultSparseMatrix._nonZeroValues, 0, NonZerosCount); - Array.Copy(lowerSparseMatrix._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, NonZerosCount, lowerSparseMatrix.NonZerosCount); + resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount; + resultStorage.Values = new double[resultStorage.ValueCount]; + resultStorage.ColumnIndices = new int[resultStorage.ValueCount]; + + Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount); + Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount); - Array.Copy(_columnIndices, 0, resultSparseMatrix._columnIndices, 0, NonZerosCount); - Array.Copy(_rowIndex, 0, resultSparseMatrix._rowIndex, 0, RowCount); + Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount); + Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount); // Copy and adjust lower column indices and rowIndex - for (int i = NonZerosCount, j = 0; i < resultSparseMatrix.NonZerosCount; i++, j++) + for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++) { - resultSparseMatrix._columnIndices[i] = lowerSparseMatrix._columnIndices[j] + ColumnCount; + resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount; } - for (int i = RowCount, j = 0; i < resultSparseMatrix.RowCount; i++, j++) + for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++) { - resultSparseMatrix._rowIndex[i] = lowerSparseMatrix._rowIndex[j] + NonZerosCount; + resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount; } } } @@ -1092,18 +967,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public static SparseMatrix Identity(int order) { - var m = new SparseMatrix(order) - { - NonZerosCount = order, - _nonZeroValues = new double[order], - _columnIndices = new int[order] - }; + var m = new SparseMatrix(order); + var mStorage = m.Storage; + + mStorage.ValueCount = order; + mStorage.Values = new double[order]; + mStorage.ColumnIndices = new int[order]; for (var i = 0; i < order; i++) { - m._nonZeroValues[i] = 1.0; - m._columnIndices[i] = i; - m._rowIndex[i] = i; + mStorage.Values[i] = 1d; + mStorage.ColumnIndices[i] = i; + mStorage.RowPointers[i] = i; } return m; @@ -1138,21 +1013,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var sparseMatrix = other as SparseMatrix; - if (sparseMatrix == null) { return base.Equals(other); } - if (NonZerosCount != sparseMatrix.NonZerosCount) + var otherStorage = sparseMatrix.Storage; + if (_storage.ValueCount != otherStorage.ValueCount) { return false; } // If all else fails, perform element wise comparison. - for (var index = 0; index < NonZerosCount; index++) + for (var index = 0; index < _storage.ValueCount; index++) { - if (!_nonZeroValues[index].AlmostEqual(sparseMatrix._nonZeroValues[index]) || _columnIndices[index] != sparseMatrix._columnIndices[index]) + if (!_storage.Values[index].AlmostEqual(otherStorage.Values[index]) || _storage.ColumnIndices[index] != otherStorage.ColumnIndices[index]) { return false; } @@ -1185,7 +1060,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CopyTo(result); } - Control.LinearAlgebraProvider.ScaleArray(2.0, _nonZeroValues, _nonZeroValues); + Control.LinearAlgebraProvider.ScaleArray(2.0, _storage.Values, _storage.Values); return; } @@ -1205,16 +1080,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double left = sparseOther; } - for (var i = 0; i < left.RowCount; i++) + var leftStorage = left.Storage; + for (var i = 0; i < leftStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = left._rowIndex[i]; - var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount; + var startIndex = leftStorage.RowPointers[i]; + var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = left._columnIndices[j]; - var resVal = left._nonZeroValues[j] + result.At(i, columnIndex); + var columnIndex = leftStorage.ColumnIndices[j]; + var resVal = leftStorage.Values[j] + result.At(i, columnIndex); result.At(i, columnIndex, resVal); } } @@ -1231,6 +1107,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { var sparseOther = other as SparseMatrix; var sparseResult = result as SparseMatrix; + if (sparseOther == null || sparseResult == null) { base.DoSubtract(other, result); @@ -1243,18 +1120,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } + var otherStorage = sparseOther.Storage; + if (ReferenceEquals(this, sparseResult)) { - for (var i = 0; i < sparseOther.RowCount; i++) + for (var i = 0; i < otherStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = sparseOther._rowIndex[i]; - var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount; + var startIndex = otherStorage.RowPointers[i]; + var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = sparseOther._columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j]; + var columnIndex = otherStorage.ColumnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j]; result.At(i, columnIndex, resVal); } } @@ -1268,16 +1147,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double sparseResult.Negate(sparseResult); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = _columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j]; + var columnIndex = columnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) + values[j]; result.At(i, columnIndex, resVal); } } @@ -1297,7 +1181,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - if (scalar == 0.0 || NonZerosCount == 0) + if (scalar == 0.0 || _storage.ValueCount == 0) { result.Clear(); return; @@ -1308,10 +1192,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var row = 0; row < RowCount; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1320,8 +1208,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var index = start; index < end; index++) { - var column = _columnIndices[index]; - result.At(row, column, _nonZeroValues[index] * scalar); + var column = columnIndices[index]; + result.At(row, column, values[index] * scalar); } } } @@ -1332,7 +1220,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._nonZeroValues[index] *= scalar); + CommonParallel.For(0, _storage.ValueCount, index => sparseResult.Storage.Values[index] *= scalar); } } @@ -1345,11 +1233,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double { result.Clear(); var columnVector = new DenseVector(other.RowCount); + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1360,10 +1254,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double // Multiply row of matrix A on column of matrix B other.Column(column, columnVector); - var sum = 0.0; + var sum = 0d; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * columnVector[_columnIndices[index]]; + sum += values[index] * columnVector[columnIndices[index]]; } result.At(row, column, sum); @@ -1378,20 +1272,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; } - var sum = 0.0; + var sum = 0d; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * rightSide[_columnIndices[index]]; + sum += values[index] * rightSide[columnIndices[index]]; } result[row] = sum; @@ -1415,11 +1314,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double } resultSparse.Clear(); + + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var otherStorage = otherSparse.Storage; + for (var j = 0; j < RowCount; j++) { // Get the begin / end index for the row - var startIndexOther = otherSparse._rowIndex[j]; - var endIndexOther = j < otherSparse._rowIndex.Length - 1 ? otherSparse._rowIndex[j + 1] : otherSparse.NonZerosCount; + var startIndexOther = otherStorage.RowPointers[j]; + var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount; if (startIndexOther == endIndexOther) { continue; @@ -1429,28 +1335,28 @@ namespace MathNet.Numerics.LinearAlgebra.Double { // Multiply row of matrix A on row of matrix B // Get the begin / end index for the row - var startIndexThis = _rowIndex[i]; - var endIndexThis = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndexThis = rowPointers[i]; + var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; if (startIndexThis == endIndexThis) { continue; } - var sum = 0.0; + var sum = 0d; for (var index = startIndexOther; index < endIndexOther; index++) { - var ind = FindItem(i, otherSparse._columnIndices[index]); + var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]); if (ind >= 0) { - sum += otherSparse._nonZeroValues[index] * _nonZeroValues[ind]; + sum += otherStorage.Values[index]*values[ind]; } } - resultSparse.SetValueAt(i, j, sum + result.At(i, j)); + resultSparse.Storage.SetValueAt(i, j, sum + result.At(i, j)); } } } - + /// /// Negate each element of this matrix and place the results into the result matrix. /// @@ -1470,18 +1376,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] * other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] * other.At(i, columnIndices[j]); + if (resVal != 0d) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1496,18 +1407,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] / other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] / other.At(i, columnIndices[j]); + if (resVal != 0d) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1520,23 +1436,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoModulus(double divisor, Matrix result) { - var denseResult = result as SparseMatrix; - - if (denseResult == null) + var sparseResult = result as SparseMatrix; + if (sparseResult == null) { base.DoModulus(divisor, result); + return; } - else + + if (!ReferenceEquals(this, result)) { - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } + CopyTo(result); + } - for (var index = 0; index < denseResult._nonZeroValues.Length; index++) - { - denseResult._nonZeroValues[index] %= divisor; - } + var resultStorage = sparseResult.Storage; + for (var index = 0; index < resultStorage.Values.Length; index++) + { + resultStorage.Values[index] %= divisor; } } @@ -1546,10 +1461,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The value at the current iteration along with its position (row, column, value). public override IEnumerable> IndexedEnumerator() { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1558,17 +1478,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var index = start; index < end; index++) { - yield return new Tuple(row, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(row, columnIndices[index], values[index]); } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < valueCount) { - for (var index = _rowIndex[lastRow]; index < NonZerosCount; index++) + for (var index = rowPointers[lastRow]; index < valueCount; index++) { - yield return new Tuple(lastRow, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(lastRow, columnIndices[index], values[index]); } } } @@ -1586,10 +1506,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double } // todo: we might be able to speed this up by caching one half of the matrix + var rowPointers = _storage.RowPointers; for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1602,11 +1523,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < _storage.ValueCount) { - if (!CheckIfOppositesAreEqual(_rowIndex[lastRow], NonZerosCount, lastRow)) + if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow)) { return false; } @@ -1625,11 +1546,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the values are equal or not. private bool CheckIfOppositesAreEqual(int start, int end, int row) { + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var index = start; index < end; index++) { - var column = _columnIndices[index]; + var column = columnIndices[index]; var opposite = At(column, row); - if (!_nonZeroValues[index].Equals(opposite)) + if (!values[index].Equals(opposite)) { return false; } diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index c88589a2..9bda0ac1 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single using System.Collections.Generic; using Generic; using Properties; + using Storage; using Threading; /// @@ -41,35 +42,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Wikipedia - CSR. /// [Serializable] - public class SparseMatrix : Matrix + public class SparseMatrix : Matrix { - /// - /// 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]; + readonly SparseCompressedRowMatrixStorage _storage; - /// - /// 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 float[] _nonZeroValues = new float[0]; - /// /// Gets the number of non zero elements in the matrix. /// /// The number of non zero elements. public int NonZerosCount { - get; - private set; + get { return _storage.ValueCount; } } - - /// - /// 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. @@ -80,9 +64,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// The number of columns. /// - public SparseMatrix(int rows, int columns) : base(rows, columns) + public SparseMatrix(int rows, int columns) + : base(rows, columns) { - _rowIndex = new int[rows]; + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0f); } /// @@ -92,8 +77,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// /// If is less than one. /// - public SparseMatrix(int order) : this(order, order) + public SparseMatrix(int order) + : base(order, order) { + _storage = new SparseCompressedRowMatrixStorage(order, order, 0f); } /// @@ -106,18 +93,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The number of columns. /// /// The value which we assign to each element of the matrix. - public SparseMatrix(int rows, int columns, float value) : this(rows, columns) + public SparseMatrix(int rows, int columns, float value) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0f); + if (value == 0.0) { return; } - NonZerosCount = rows * columns; - _nonZeroValues = new float[NonZerosCount]; - _columnIndices = new int[NonZerosCount]; + var rowPointers = _storage.RowPointers; + var valueCount = _storage.ValueCount = rows * columns; + var columnIndices = _storage.ColumnIndices = new int[valueCount]; + var values = _storage.Values = new float[valueCount]; - for (int i = 0, j = 0; i < _nonZeroValues.Length; i++, j++) + for (int i = 0, j = 0; i < values.Length; i++, j++) { // Reset column position to "0" if (j == columns) @@ -125,14 +116,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single j = 0; } - _nonZeroValues[i] = value; - _columnIndices[i] = j; + values[i] = value; + columnIndices[i] = j; } - + // Set proper row pointers - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { - _rowIndex[i] = ((i + 1) * columns) - columns; + rowPointers[i] = ((i + 1) * columns) - columns; } } @@ -144,8 +135,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The one dimensional array to create this matrix from. This array should store the matrix in column-major order. see: http://en.wikipedia.org/wiki/Column-major_order /// If length is less than * . /// - public SparseMatrix(int rows, int columns, float[] array) : this(rows, columns) + public SparseMatrix(int rows, int columns, float[] array) + : base(rows, columns) { + _storage = new SparseCompressedRowMatrixStorage(rows, columns, 0f); + if (rows * columns > array.Length) { throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions); @@ -155,7 +149,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, array[i + (j * rows)]); + _storage.SetValueAt(i, j, array[i + (j * rows)]); } } } @@ -164,16 +158,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Initializes a new instance of the class from a 2D array. /// /// The 2D array to create this matrix from. - public SparseMatrix(float[,] array) : this(array.GetLength(0), array.GetLength(1)) + public SparseMatrix(float[,] array) + : base(array.GetLength(0), array.GetLength(1)) { - var rows = array.GetLength(0); - var columns = array.GetLength(1); + _storage = new SparseCompressedRowMatrixStorage(array.GetLength(0), array.GetLength(1), 0f); - for (var i = 0; i < rows; i++) + for (var i = 0; i < _storage.RowCount; i++) { - for (var j = 0; j < columns; j++) + for (var j = 0; j < _storage.ColumnCount; j++) { - SetValueAt(i, j, array[i, j]); + _storage.SetValueAt(i, j, array[i, j]); } } } @@ -186,6 +180,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single public SparseMatrix(Matrix matrix) : this(matrix.RowCount, matrix.ColumnCount) { + _storage = new SparseCompressedRowMatrixStorage(matrix.RowCount, matrix.ColumnCount, 0f); + var sparseMatrix = matrix as SparseMatrix; var rows = matrix.RowCount; @@ -197,23 +193,28 @@ namespace MathNet.Numerics.LinearAlgebra.Single { for (var j = 0; j < columns; j++) { - SetValueAt(i, j, matrix.At(i, j)); + _storage.SetValueAt(i, j, matrix.At(i, j)); } } } else { - NonZerosCount = sparseMatrix.NonZerosCount; - _rowIndex = new int[rows]; - _columnIndices = new int[NonZerosCount]; - _nonZeroValues = new float[NonZerosCount]; + var matrixStorage = sparseMatrix.Storage; + var valueCount = _storage.ValueCount = matrixStorage.ValueCount; + _storage.ColumnIndices = new int[valueCount]; + _storage.Values = new float[valueCount]; - Buffer.BlockCopy(sparseMatrix._nonZeroValues, 0, _nonZeroValues, 0, NonZerosCount * Constants.SizeOfFloat); - Buffer.BlockCopy(sparseMatrix._columnIndices, 0, _columnIndices, 0, NonZerosCount * Constants.SizeOfInt); - Buffer.BlockCopy(sparseMatrix._rowIndex, 0, _rowIndex, 0, rows * Constants.SizeOfInt); + Buffer.BlockCopy(matrixStorage.Values, 0, _storage.Values, 0, valueCount * Constants.SizeOfFloat); + Buffer.BlockCopy(matrixStorage.ColumnIndices, 0, _storage.ColumnIndices, 0, valueCount * Constants.SizeOfInt); + Buffer.BlockCopy(matrixStorage.RowPointers, 0, _storage.RowPointers, 0, rows * Constants.SizeOfInt); } } + internal SparseCompressedRowMatrixStorage Storage + { + get { return _storage; } + } + /// /// Creates a SparseMatrix for the given number of rows and columns. /// @@ -291,15 +292,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Where to store the lower triangle. private void LowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row >= _columnIndices[j]) + if (row >= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -353,15 +359,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Where to store the lower triangle. private void UpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row <= _columnIndices[j]) + if (row <= columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -420,18 +431,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single var result = (SparseMatrix)CreateMatrix(rowCount, columnCount); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (int i = rowIndex, row = 0; i < rowMax; i++, row++) { - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (int j = startIndex; j < endIndex; j++) { // check if the column index is in the range - if ((_columnIndices[j] >= columnIndex) && (_columnIndices[j] < columnIndex + columnCount)) + if ((columnIndices[j] >= columnIndex) && (columnIndices[j] < columnIndex + columnCount)) { - var column = _columnIndices[j] - columnIndex; - result.SetValueAt(row, column, _nonZeroValues[j]); + var column = columnIndices[j] - columnIndex; + result._storage.SetValueAt(row, column, values[j]); } } } @@ -488,15 +504,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Where to store the lower triangle. private void StrictlyLowerTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row > _columnIndices[j]) + if (row > columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -551,15 +572,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Where to store the lower triangle. private void StrictlyUpperTriangleImpl(Matrix result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < result.RowCount; row++) { - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - if (row < _columnIndices[j]) + if (row < columnIndices[j]) { - result.At(row, _columnIndices[j], _nonZeroValues[j]); + result.At(row, columnIndices[j], values[j]); } } } @@ -576,13 +602,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// An array containing the matrix's elements. public override float[] ToColumnWiseArray() { + var values = _storage.Values; var ret = new float[RowCount * ColumnCount]; for (var j = 0; j < ColumnCount; j++) { for (var i = 0; i < RowCount; i++) { - var index = FindItem(i, j); - ret[(j * RowCount) + i] = index >= 0 ? _nonZeroValues[index] : 0.0f; + var index = _storage.FindItem(i, j); + ret[(j * RowCount) + i] = index >= 0 ? values[index] : 0.0f; } } @@ -603,8 +630,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// public override float At(int row, int column) { - var index = FindItem(row, column); - return index >= 0 ? _nonZeroValues[index] : 0.0f; + return _storage.GetValueAt(row, column); } /// @@ -621,166 +647,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// public override void At(int row, int column, float value) { - 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, float value) - { - var index = FindItem(row, column); - if (index >= 0) - { - // Non-zero item found in matrix - if (value == 0.0) - { - // Delete existing item - DeleteItemByIndex(index, row); - } - else - { - // Update item - _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 ((NonZerosCount == _nonZeroValues.Length) && (NonZerosCount < ((long)RowCount * 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(_nonZeroValues.Length + GrowthSize(), (long)RowCount * ColumnCount); - if (size > int.MaxValue) - { - throw new NotSupportedException(Resources.TooManyElements); - } - - Array.Resize(ref _nonZeroValues, (int)size); - Array.Resize(ref _columnIndices, (int)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 = NonZerosCount - 1; i > index - 1; i--) - { - _nonZeroValues[i + 1] = _nonZeroValues[i]; - _columnIndices[i + 1] = _columnIndices[i]; - } - - // Add the value and the column index - _nonZeroValues[index] = value; - _columnIndices[index] = column; - - // increase the number of non-zero numbers by one - 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 < _rowIndex.Length; i++) - { - _rowIndex[i] += 1; - } - } + _storage.SetValueAt(row, column, value); } - /// - /// 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 < NonZerosCount; i++) - { - _nonZeroValues[i - 1] = _nonZeroValues[i]; - _columnIndices[i - 1] = _columnIndices[i]; - } - - // Decrease value in Row - for (var i = row + 1; i < _rowIndex.Length; i++) - { - _rowIndex[i] -= 1; - } - - 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 ((NonZerosCount > 1024) && (NonZerosCount < _nonZeroValues.Length / 2)) - { - Array.Resize(ref _nonZeroValues, NonZerosCount); - Array.Resize(ref _columnIndices, 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 = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; - return Array.BinarySearch(_columnIndices, startIndex, endIndex - startIndex, column); - } - - /// - /// Calculates the amount with which to grow the storage array's if they need to be - /// increased in size. - /// - /// The amount grown. - private int GrowthSize() - { - int delta; - if (_nonZeroValues.Length > 1024) - { - delta = _nonZeroValues.Length / 4; - } - else - { - if (_nonZeroValues.Length > 256) - { - delta = 512; - } - else - { - delta = _nonZeroValues.Length > 64 ? 128 : 32; - } - } - - return delta; - } - #endregion - /// /// Sets all values to zero. /// public override void Clear() { - NonZerosCount = 0; - Array.Clear(_rowIndex, 0, _rowIndex.Length); + _storage.Clear(); } /// @@ -814,16 +689,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw DimensionsDontMatch(this, target, "target"); } + var targetStorage = sparseTarget.Storage; + // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value - sparseTarget._nonZeroValues = new float[NonZerosCount]; - sparseTarget._columnIndices = new int[NonZerosCount]; - sparseTarget.NonZerosCount = NonZerosCount; + targetStorage.ValueCount = _storage.ValueCount; + targetStorage.Values = new float[_storage.ValueCount]; + targetStorage.ColumnIndices = new int[_storage.ValueCount]; - if (NonZerosCount != 0) + if (_storage.ValueCount != 0) { - Buffer.BlockCopy(_nonZeroValues, 0, sparseTarget._nonZeroValues, 0, NonZerosCount*Constants.SizeOfFloat); - Buffer.BlockCopy(_columnIndices, 0, sparseTarget._columnIndices, 0, NonZerosCount*Constants.SizeOfInt); - Buffer.BlockCopy(_rowIndex, 0, sparseTarget._rowIndex, 0, RowCount*Constants.SizeOfInt); + Buffer.BlockCopy(_storage.Values, 0, targetStorage.Values, 0, _storage.ValueCount * Constants.SizeOfFloat); + Buffer.BlockCopy(_storage.ColumnIndices, 0, targetStorage.ColumnIndices, 0, _storage.ValueCount * Constants.SizeOfInt); + Buffer.BlockCopy(_storage.RowPointers, 0, targetStorage.RowPointers, 0, RowCount * Constants.SizeOfInt); } } @@ -835,14 +712,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// public override int GetHashCode() { - var hashNum = Math.Min(NonZerosCount, 25); + var values = _storage.Values; + var hashNum = Math.Min(_storage.ValueCount, 25); long hash = 0; for (var i = 0; i < hashNum; i++) { #if PORTABLE - hash ^= Precision.DoubleToInt64Bits(_nonZeroValues[i]); + hash ^= Precision.DoubleToInt64Bits(values[i]); #else - hash ^= BitConverter.DoubleToInt64Bits(_nonZeroValues[i]); + hash ^= BitConverter.DoubleToInt64Bits(values[i]); #endif } @@ -855,18 +733,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The transpose of this matrix. public override Matrix Transpose() { - var ret = new SparseMatrix(ColumnCount, RowCount) - { - _columnIndices = new int[NonZerosCount], - _nonZeroValues = new float[NonZerosCount] - }; + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var ret = new SparseMatrix(ColumnCount, RowCount); + var retStorage = ret.Storage; + retStorage.ColumnIndices = new int[valueCount]; + retStorage.Values = new float[valueCount]; // Do an 'inverse' CopyTo iterate over the rows - for (var i = 0; i < _rowIndex.Length; i++) + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -877,7 +759,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var j = startIndex; j < endIndex; j++) { - ret.SetValueAt(_columnIndices[j], i, _nonZeroValues[j]); + retStorage.SetValueAt(columnIndices[j], i, values[j]); } } @@ -889,15 +771,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override float FrobeniusNorm() { var transpose = (SparseMatrix)Transpose(); - var aat = this * transpose; + var aat = (this * transpose).Storage; - var norm = 0.0f; + var norm = 0f; - for (var i = 0; i < aat._rowIndex.Length; i++) + for (var i = 0; i < aat.RowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = aat._rowIndex[i]; - var endIndex = i < aat._rowIndex.Length - 1 ? aat._rowIndex[i + 1] : aat.NonZerosCount; + var startIndex = aat.RowPointers[i]; + var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount; // Get the values for the current row if (startIndex == endIndex) @@ -908,27 +790,30 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var j = startIndex; j < endIndex; j++) { - if (i == aat._columnIndices[j]) + if (i == aat.ColumnIndices[j]) { - norm += Math.Abs(aat._nonZeroValues[j]); + norm += Math.Abs(aat.Values[j]); } } } - norm = Convert.ToSingle(Math.Sqrt(norm)); - return norm; + return Convert.ToSingle(Math.Sqrt(norm)); } /// Calculates the infinity norm of this matrix. /// The infinity norm of this matrix. public override float InfinityNorm() { - var norm = 0.0f; - for (var i = 0; i < _rowIndex.Length; i++) + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var norm = 0f; + for (var i = 0; i < rowPointers.Length; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; // Get the values for the current row if (startIndex == endIndex) @@ -937,10 +822,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single continue; } - var s = 0.0f; + var s = 0f; for (var j = startIndex; j < endIndex; j++) { - s += Math.Abs(_nonZeroValues[j]); + s += Math.Abs(values[j]); } norm = Math.Max(norm, s); @@ -997,9 +882,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result"); } + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + // Determine bounds in columnIndices array where this item should be searched (using rowIndex) - var startIndex = _rowIndex[rowIndex]; - var endIndex = rowIndex < _rowIndex.Length - 1 ? _rowIndex[rowIndex + 1] : NonZerosCount; + var startIndex = rowPointers[rowIndex]; + var endIndex = rowIndex < rowPointers.Length - 1 ? rowPointers[rowIndex + 1] : valueCount; if (startIndex == endIndex) { @@ -1011,8 +900,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (int i = columnIndex, j = 0; i < columnIndex + length; i++, j++) { // Copy code from At(row, column) to avoid unnecessary lock - var index = FindItem(rowIndex, i); - result[j] = index >= 0 ? _nonZeroValues[index] : 0.0f; + var index = _storage.FindItem(rowIndex, i); + result[j] = index >= 0 ? values[index] : 0f; } } } @@ -1036,30 +925,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single } else { - if (resultSparseMatrix.RowCount != RowCount + lowerSparseMatrix.RowCount || resultSparseMatrix.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) + var resultStorage = resultSparseMatrix.Storage; + var lowerStorage = lowerSparseMatrix.Storage; + + if (resultStorage.RowCount != RowCount + lowerStorage.RowCount || resultStorage.ColumnCount != ColumnCount + lowerSparseMatrix.ColumnCount) { - throw DimensionsDontMatch(this, lowerSparseMatrix, resultSparseMatrix, "result"); + throw DimensionsDontMatch(this, lowerSparseMatrix, resultSparseMatrix); } - resultSparseMatrix.NonZerosCount = NonZerosCount + lowerSparseMatrix.NonZerosCount; - resultSparseMatrix._nonZeroValues = new float[resultSparseMatrix.NonZerosCount]; - resultSparseMatrix._columnIndices = new int[resultSparseMatrix.NonZerosCount]; - - Array.Copy(_nonZeroValues, 0, resultSparseMatrix._nonZeroValues, 0, NonZerosCount); - Array.Copy(lowerSparseMatrix._nonZeroValues, 0, resultSparseMatrix._nonZeroValues, NonZerosCount, lowerSparseMatrix.NonZerosCount); + resultStorage.ValueCount = _storage.ValueCount + lowerStorage.ValueCount; + resultStorage.Values = new float[resultStorage.ValueCount]; + resultStorage.ColumnIndices = new int[resultStorage.ValueCount]; - Array.Copy(_columnIndices, 0, resultSparseMatrix._columnIndices, 0, NonZerosCount); - Array.Copy(_rowIndex, 0, resultSparseMatrix._rowIndex, 0, RowCount); + Array.Copy(_storage.Values, 0, resultStorage.Values, 0, _storage.ValueCount); + Array.Copy(lowerStorage.Values, 0, resultStorage.Values, _storage.ValueCount, lowerStorage.ValueCount); + + Array.Copy(_storage.ColumnIndices, 0, resultStorage.ColumnIndices, 0, _storage.ValueCount); + Array.Copy(_storage.RowPointers, 0, resultStorage.RowPointers, 0, RowCount); // Copy and adjust lower column indices and rowIndex - for (int i = NonZerosCount, j = 0; i < resultSparseMatrix.NonZerosCount; i++, j++) + for (int i = _storage.ValueCount, j = 0; i < resultStorage.ValueCount; i++, j++) { - resultSparseMatrix._columnIndices[i] = lowerSparseMatrix._columnIndices[j] + ColumnCount; + resultStorage.ColumnIndices[i] = lowerStorage.ColumnIndices[j] + ColumnCount; } - for (int i = RowCount, j = 0; i < resultSparseMatrix.RowCount; i++, j++) + for (int i = RowCount, j = 0; i < resultStorage.RowCount; i++, j++) { - resultSparseMatrix._rowIndex[i] = lowerSparseMatrix._rowIndex[j] + NonZerosCount; + resultStorage.RowPointers[i] = lowerStorage.RowPointers[j] + _storage.ValueCount; } } } @@ -1075,18 +967,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// public static SparseMatrix Identity(int order) { - var m = new SparseMatrix(order) - { - NonZerosCount = order, - _nonZeroValues = new float[order], - _columnIndices = new int[order] - }; + var m = new SparseMatrix(order); + var mStorage = m.Storage; + + mStorage.ValueCount = order; + mStorage.Values = new float[order]; + mStorage.ColumnIndices = new int[order]; for (var i = 0; i < order; i++) { - m._nonZeroValues[i] = 1.0f; - m._columnIndices[i] = i; - m._rowIndex[i] = i; + mStorage.Values[i] = 1f; + mStorage.ColumnIndices[i] = i; + mStorage.RowPointers[i] = i; } return m; @@ -1121,21 +1013,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single } var sparseMatrix = other as SparseMatrix; - if (sparseMatrix == null) { return base.Equals(other); } - if (NonZerosCount != sparseMatrix.NonZerosCount) + var otherStorage = sparseMatrix.Storage; + if (_storage.ValueCount != otherStorage.ValueCount) { return false; } // If all else fails, perform element wise comparison. - for (var index = 0; index < NonZerosCount; index++) + for (var index = 0; index < _storage.ValueCount; index++) { - if (!_nonZeroValues[index].AlmostEqual(sparseMatrix._nonZeroValues[index]) || _columnIndices[index] != sparseMatrix._columnIndices[index]) + if (!_storage.Values[index].AlmostEqual(otherStorage.Values[index]) || _storage.ColumnIndices[index] != otherStorage.ColumnIndices[index]) { return false; } @@ -1168,7 +1060,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single CopyTo(result); } - Control.LinearAlgebraProvider.ScaleArray(2.0f, _nonZeroValues, _nonZeroValues); + Control.LinearAlgebraProvider.ScaleArray(2.0f, _storage.Values, _storage.Values); return; } @@ -1188,16 +1080,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single left = sparseOther; } - for (var i = 0; i < left.RowCount; i++) + var leftStorage = left.Storage; + for (var i = 0; i < leftStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = left._rowIndex[i]; - var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount; + var startIndex = leftStorage.RowPointers[i]; + var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = left._columnIndices[j]; - var resVal = left._nonZeroValues[j] + result.At(i, columnIndex); + var columnIndex = leftStorage.ColumnIndices[j]; + var resVal = leftStorage.Values[j] + result.At(i, columnIndex); result.At(i, columnIndex, resVal); } } @@ -1226,18 +1119,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single return; } + var otherStorage = sparseOther.Storage; + if (ReferenceEquals(this, sparseResult)) { - for (var i = 0; i < sparseOther.RowCount; i++) + for (var i = 0; i < otherStorage.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = sparseOther._rowIndex[i]; - var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount; + var startIndex = otherStorage.RowPointers[i]; + var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = sparseOther._columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j]; + var columnIndex = otherStorage.ColumnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j]; result.At(i, columnIndex, resVal); } } @@ -1251,16 +1146,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single sparseResult.Negate(sparseResult); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var columnIndex = _columnIndices[j]; - var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j]; + var columnIndex = columnIndices[j]; + var resVal = sparseResult.At(i, columnIndex) + values[j]; result.At(i, columnIndex, resVal); } } @@ -1291,10 +1191,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var row = 0; row < RowCount; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1303,8 +1207,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var index = start; index < end; index++) { - var column = _columnIndices[index]; - result.At(row, column, _nonZeroValues[index] * scalar); + var column = columnIndices[index]; + result.At(row, column, values[index] * scalar); } } } @@ -1315,7 +1219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._nonZeroValues[index] *= scalar); + CommonParallel.For(0, NonZerosCount, index => sparseResult.Storage.Values[index] *= scalar); } } @@ -1328,11 +1232,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single { result.Clear(); var columnVector = new DenseVector(other.RowCount); + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; @@ -1342,11 +1252,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single { // Multiply row of matrix A on column of matrix B other.Column(column, columnVector); - - var sum = 0.0f; + + var sum = 0f; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * columnVector[_columnIndices[index]]; + sum += values[index] * columnVector[columnIndices[index]]; } result.At(row, column, sum); @@ -1361,20 +1271,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector result) { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount; row++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[row]; - var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + var startIndex = rowPointers[row]; + var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount; if (startIndex == endIndex) { continue; } - var sum = 0.0f; + var sum = 0f; for (var index = startIndex; index < endIndex; index++) { - sum += _nonZeroValues[index] * rightSide[_columnIndices[index]]; + sum += values[index] * rightSide[columnIndices[index]]; } result[row] = sum; @@ -1398,11 +1313,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single } resultSparse.Clear(); + + var rowPointers = _storage.RowPointers; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + var otherStorage = otherSparse.Storage; + for (var j = 0; j < RowCount; j++) { // Get the begin / end index for the row - var startIndexOther = otherSparse._rowIndex[j]; - var endIndexOther = j < otherSparse._rowIndex.Length - 1 ? otherSparse._rowIndex[j + 1] : otherSparse.NonZerosCount; + var startIndexOther = otherStorage.RowPointers[j]; + var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount; if (startIndexOther == endIndexOther) { continue; @@ -1412,24 +1334,24 @@ namespace MathNet.Numerics.LinearAlgebra.Single { // Multiply row of matrix A on row of matrix B // Get the begin / end index for the row - var startIndexThis = _rowIndex[i]; - var endIndexThis = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndexThis = rowPointers[i]; + var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; if (startIndexThis == endIndexThis) { continue; } - var sum = 0.0f; + var sum = 0f; for (var index = startIndexOther; index < endIndexOther; index++) { - var ind = FindItem(i, otherSparse._columnIndices[index]); + var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]); if (ind >= 0) { - sum += otherSparse._nonZeroValues[index] * _nonZeroValues[ind]; + sum += otherStorage.Values[index]*values[ind]; } } - resultSparse.SetValueAt(i, j, sum + result.At(i, j)); + resultSparse.Storage.SetValueAt(i, j, sum + result.At(i, j)); } } } @@ -1453,18 +1375,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] * other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] * other.At(i, columnIndices[j]); + if (resVal != 0f) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1479,18 +1406,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single { result.Clear(); + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var i = 0; i < other.RowCount; i++) { // Get the begin / end index for the current row - var startIndex = _rowIndex[i]; - var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount; + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; for (var j = startIndex; j < endIndex; j++) { - var resVal = _nonZeroValues[j] / other.At(i, _columnIndices[j]); - if (resVal != 0.0) + var resVal = values[j] / other.At(i, columnIndices[j]); + if (resVal != 0f) { - result.At(i, _columnIndices[j], resVal); + result.At(i, columnIndices[j], resVal); } } } @@ -1503,23 +1435,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoModulus(float divisor, Matrix result) { - var denseResult = result as SparseMatrix; - - if (denseResult == null) + var sparseResult = result as SparseMatrix; + if (sparseResult == null) { base.DoModulus(divisor, result); + return; } - else + + if (!ReferenceEquals(this, result)) { - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } + CopyTo(result); + } - for (var index = 0; index < denseResult._nonZeroValues.Length; index++) - { - denseResult._nonZeroValues[index] %= divisor; - } + var resultStorage = sparseResult.Storage; + for (var index = 0; index < resultStorage.Values.Length; index++) + { + resultStorage.Values[index] %= divisor; } } @@ -1529,10 +1460,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The value at the current iteration along with its position (row, column, value). public override IEnumerable> IndexedEnumerator() { + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1541,17 +1477,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var index = start; index < end; index++) { - yield return new Tuple(row, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(row, columnIndices[index], values[index]); } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < NonZerosCount) { - for (var index = _rowIndex[lastRow]; index < NonZerosCount; index++) + for (var index = rowPointers[lastRow]; index < NonZerosCount; index++) { - yield return new Tuple(lastRow, _columnIndices[index], _nonZeroValues[index]); + yield return new Tuple(lastRow, columnIndices[index], values[index]); } } } @@ -1569,10 +1505,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single } // todo: we might be able to speed this up by caching one half of the matrix + var rowPointers = _storage.RowPointers; for (var row = 0; row < RowCount - 1; row++) { - var start = _rowIndex[row]; - var end = _rowIndex[row + 1]; + var start = rowPointers[row]; + var end = rowPointers[row + 1]; if (start == end) { @@ -1585,11 +1522,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single } } - var lastRow = _rowIndex.Length - 1; + var lastRow = rowPointers.Length - 1; - if (_rowIndex[lastRow] < NonZerosCount) + if (rowPointers[lastRow] < NonZerosCount) { - if (!CheckIfOppositesAreEqual(_rowIndex[lastRow], NonZerosCount, lastRow)) + if (!CheckIfOppositesAreEqual(rowPointers[lastRow], NonZerosCount, lastRow)) { return false; } @@ -1608,11 +1545,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// If the values are equal or not. private bool CheckIfOppositesAreEqual(int start, int end, int row) { + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + for (var index = start; index < end; index++) { - var column = _columnIndices[index]; + var column = columnIndices[index]; var opposite = At(column, row); - if (!_nonZeroValues[index].Equals(opposite)) + if (!values[index].Equals(opposite)) { return false; } diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs index 70be6cf5..9141a0fd 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs @@ -6,14 +6,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage internal class DenseColumnMajorMatrixStorage where T : struct, IEquatable, IFormattable { - public int RowCount { get; private set; } - public int ColumnCount { get; private set; } - public T[] Data { get; private set; } + // [ruegg] public fields are OK here + + public readonly int RowCount; + public readonly int ColumnCount; + + public readonly T[] Data; internal DenseColumnMajorMatrixStorage(int rows, int columns) { RowCount = rows; ColumnCount = columns; + Data = new T[rows * columns]; } @@ -21,6 +25,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { RowCount = rows; ColumnCount = columns; + Data = data; } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs new file mode 100644 index 00000000..8d983193 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -0,0 +1,229 @@ +using System; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.LinearAlgebra.Storage +{ + internal class SparseCompressedRowMatrixStorage + where T : struct, IEquatable, IFormattable + { + // [ruegg] public fields are OK here + + public readonly int RowCount; + public readonly int ColumnCount; + readonly T _zero; + + /// + /// 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" + /// + public readonly int[] RowPointers; + + /// + /// 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. + /// + public int[] ColumnIndices; + + /// + /// 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. + /// + public T[] Values; + + /// + /// Gets the number of non zero elements in the matrix. + /// + /// The number of non zero elements. + public int ValueCount; + + internal SparseCompressedRowMatrixStorage(int rows, int columns, T zero) + { + RowCount = rows; + ColumnCount = columns; + _zero = zero; + + RowPointers = new int[rows]; + ColumnIndices = new int[0]; + Values = new T[0]; + ValueCount = 0; + } + + public T this[int row, int column] + { + get { return GetValueAt(row, column); } + set { SetValueAt(row, column, value); } + } + + public void Clear() + { + ValueCount = 0; + Array.Clear(RowPointers, 0, RowPointers.Length); + } + + /// + /// Retrieves the requested element without range checking. + /// + /// + /// The row of the element. + /// + /// + /// The column of the element. + /// + /// + /// The requested element. + /// + public T GetValueAt(int row, int column) + { + var index = FindItem(row, column); + return index >= 0 ? Values[index] : _zero; + } + + /// + /// 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 + public void SetValueAt(int row, int column, T value) + { + var index = FindItem(row, column); + if (index >= 0) + { + // Non-zero item found in matrix + if (_zero.Equals(value)) + { + // Delete existing item + DeleteItemByIndex(index, row); + } + else + { + // Update item + Values[index] = value; + } + } + else + { + // Item not found. Add new value + if (_zero.Equals(value)) + { + return; + } + + index = ~index; + + // Check if the storage needs to be increased + if ((ValueCount == Values.Length) && (ValueCount < ((long)RowCount * 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(Values.Length + GrowthSize(), (long)RowCount * ColumnCount); + if (size > int.MaxValue) + { + throw new NotSupportedException(Resources.TooManyElements); + } + + Array.Resize(ref Values, (int)size); + Array.Resize(ref ColumnIndices, (int)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 = ValueCount - 1; i > index - 1; i--) + { + Values[i + 1] = Values[i]; + ColumnIndices[i + 1] = ColumnIndices[i]; + } + + // Add the value and the column index + Values[index] = value; + ColumnIndices[index] = column; + + // increase the number of non-zero numbers by one + ValueCount += 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 < RowPointers.Length; i++) + { + RowPointers[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 < ValueCount; i++) + { + Values[i - 1] = Values[i]; + ColumnIndices[i - 1] = ColumnIndices[i]; + } + + // Decrease value in Row + for (var i = row + 1; i < RowPointers.Length; i++) + { + RowPointers[i] -= 1; + } + + ValueCount -= 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 ((ValueCount > 1024) && (ValueCount < Values.Length / 2)) + { + Array.Resize(ref Values, ValueCount); + Array.Resize(ref ColumnIndices, ValueCount); + } + } + + /// + /// 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 + public int FindItem(int row, int column) + { + // Determin bounds in columnIndices array where this item should be searched (using rowIndex) + var startIndex = RowPointers[row]; + var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount; + return Array.BinarySearch(ColumnIndices, startIndex, endIndex - startIndex, column); + } + + /// + /// Calculates the amount with which to grow the storage array's if they need to be + /// increased in size. + /// + /// The amount grown. + private int GrowthSize() + { + int delta; + if (Values.Length > 1024) + { + delta = Values.Length / 4; + } + else + { + if (Values.Length > 256) + { + delta = 512; + } + else + { + delta = Values.Length > 64 ? 128 : 32; + } + } + + return delta; + } + } +} diff --git a/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs index 29a4f727..5f6779a8 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs @@ -6,25 +6,30 @@ namespace MathNet.Numerics.LinearAlgebra.Storage internal class SparseDiagonalMatrixStorage where T : struct, IEquatable, IFormattable { - public int RowCount { get; private set; } - public int ColumnCount { get; private set; } - public T[] Data { get; private set; } + // [ruegg] public fields are OK here + + public readonly int RowCount; + public readonly int ColumnCount; readonly T _zero; + public readonly T[] Data; + internal SparseDiagonalMatrixStorage(int rows, int columns, T zero) { RowCount = rows; ColumnCount = columns; - Data = new T[Math.Min(rows, columns)]; _zero = zero; + + Data = new T[Math.Min(rows, columns)]; } internal SparseDiagonalMatrixStorage(int rows, int columns, T zero, T[] data) { RowCount = rows; ColumnCount = columns; - Data = data; _zero = zero; + + Data = data; } public T this[int row, int column] diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 6fdfb964..24e2c91b 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -331,6 +331,7 @@ + diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index 3c81beb6..da3a7848 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -891,6 +891,9 @@ LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs + + LinearAlgebra\Storage\SparseCompressedRowMatrixStorage.cs + LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs