diff --git a/src/Numerics.Tests/LinearAlgebraTests/MatrixStructureTheory.cs b/src/Numerics.Tests/LinearAlgebraTests/MatrixStructureTheory.cs index 6a32c910..58ce9e38 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/MatrixStructureTheory.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/MatrixStructureTheory.cs @@ -27,11 +27,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using System.Linq; using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.LinearAlgebra.Storage; using NUnit.Framework; +using System; +using System.Linq; namespace MathNet.Numerics.UnitTests.LinearAlgebraTests { @@ -611,13 +610,49 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests } } - var matrix = Matrix.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues); - Assert.That(matrix.GetType().Name, Is.EqualTo("SparseMatrix")); - Assert.That(matrix.RowCount, Is.EqualTo(3)); - Assert.That(matrix.ColumnCount, Is.EqualTo(4)); + var A = Matrix.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues); + Assert.That(A.GetType().Name, Is.EqualTo("SparseMatrix")); + Assert.That(A.RowCount, Is.EqualTo(3)); + Assert.That(A.ColumnCount, Is.EqualTo(4)); + + cooRowIndices.Reverse(); + cooColumnIndices.Reverse(); + cooValues.Reverse(); + var B = Matrix.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues); + for (int j = 0; j < 4; j++) + { for (int i = 0; i < 3; i++) - Assert.That(matrix[i, j], Is.EqualTo(rows[i][j])); + { + Assert.That(A[i, j], Is.EqualTo(rows[i][j])); + Assert.That(B[i, j], Is.EqualTo(rows[i][j])); + } + } + } + + [Test] + public void CanCreateSparseFromNonOrderedDuplicatedCoordinateFormat() + { + int rowCount = 2, columnCount = 2, valueCount = 5; + var cooRowIndices = new int[5] { 1, 0, 1, 0, 1 }; + var cooColumnIndices = new int[5] { 1, 0, 0, 1, 1 }; + var cooValues = Vector.Build.Random(5, 0).ToArray(); + + var A = Matrix.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues); + + cooRowIndices.Reverse(); + cooColumnIndices.Reverse(); + cooValues.Reverse(); + + var B = Matrix.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues); + + for (int j = 0; j < columnCount; j++) + { + for (int i = 0; i < rowCount; i++) + { + Assert.That(A[i, j], Is.EqualTo(B[i, j])); + } + } } [Test] @@ -654,13 +689,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests csrRowPointers[i] += csrRowPointers[i - 1]; } - var matrix = Matrix.Build.SparseFromCompressedSparseRowFormat(rowCount, columnCount, valueCount, csrRowPointers, csrColumnIndices, csrValues); - Assert.That(matrix.GetType().Name, Is.EqualTo("SparseMatrix")); - Assert.That(matrix.RowCount, Is.EqualTo(3)); - Assert.That(matrix.ColumnCount, Is.EqualTo(4)); + var A = Matrix.Build.SparseFromCompressedSparseRowFormat(rowCount, columnCount, valueCount, csrRowPointers, csrColumnIndices, csrValues); + Assert.That(A.GetType().Name, Is.EqualTo("SparseMatrix")); + Assert.That(A.RowCount, Is.EqualTo(3)); + Assert.That(A.ColumnCount, Is.EqualTo(4)); + + csrRowPointers.Reverse(); + csrColumnIndices.Reverse(); + csrValues.Reverse(); + var B = Matrix.Build.SparseFromCompressedSparseRowFormat(rowCount, columnCount, valueCount, csrRowPointers, csrColumnIndices, csrValues); + for (int j = 0; j < 4; j++) + { for (int i = 0; i < 3; i++) - Assert.That(matrix[i, j], Is.EqualTo(rows[i][j])); + { + Assert.That(A[i, j], Is.EqualTo(rows[i][j])); + Assert.That(B[i, j], Is.EqualTo(rows[i][j])); + } + } } [Test] @@ -697,13 +743,24 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests cscColumnPointers[i] += cscColumnPointers[i - 1]; } - var matrix = Matrix.Build.SparseFromCompressedSparseColumnFormat(rowCount, columnCount, valueCount, cscRowIndices, cscColumnPointers, cscValues); - Assert.That(matrix.GetType().Name, Is.EqualTo("SparseMatrix")); - Assert.That(matrix.RowCount, Is.EqualTo(3)); - Assert.That(matrix.ColumnCount, Is.EqualTo(4)); + var A = Matrix.Build.SparseFromCompressedSparseColumnFormat(rowCount, columnCount, valueCount, cscRowIndices, cscColumnPointers, cscValues); + Assert.That(A.GetType().Name, Is.EqualTo("SparseMatrix")); + Assert.That(A.RowCount, Is.EqualTo(3)); + Assert.That(A.ColumnCount, Is.EqualTo(4)); + + cscRowIndices.Reverse(); + cscColumnPointers.Reverse(); + cscValues.Reverse(); + var B = Matrix.Build.SparseFromCompressedSparseColumnFormat(rowCount, columnCount, valueCount, cscRowIndices, cscColumnPointers, cscValues); + for (int j = 0; j < 4; j++) + { for (int i = 0; i < 3; i++) - Assert.That(matrix[i, j], Is.EqualTo(rows[i][j])); + { + Assert.That(A[i, j], Is.EqualTo(rows[i][j])); + Assert.That(B[i, j], Is.EqualTo(rows[i][j])); + } + } } [Test] diff --git a/src/Numerics/LinearAlgebra/Builder.cs b/src/Numerics/LinearAlgebra/Builder.cs index 12527d8f..0a43c00b 100644 --- a/src/Numerics/LinearAlgebra/Builder.cs +++ b/src/Numerics/LinearAlgebra/Builder.cs @@ -48,8 +48,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double return new DenseMatrix(storage); } - public override Matrix Sparse(SparseCompressedRowMatrixStorage storage) + public override Matrix Sparse(SparseCompressedRowMatrixStorage storage, bool cleanup = false) { + if (cleanup) + { + SumDuplicates(storage); + } return new SparseMatrix(storage); } @@ -73,6 +77,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double new ResidualStopCriterion(1e-12) }; } + + internal override double AddEntries(double x, double y) + { + return x + y; + } } internal class VectorBuilder : VectorBuilder @@ -111,8 +120,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single return new DenseMatrix(storage); } - public override Matrix Sparse(SparseCompressedRowMatrixStorage storage) + public override Matrix Sparse(SparseCompressedRowMatrixStorage storage, bool cleanup = false) { + if (cleanup) + { + SumDuplicates(storage); + } return new SparseMatrix(storage); } @@ -136,6 +149,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single new ResidualStopCriterion(1e-6) }; } + + internal override float AddEntries(float x, float y) + { + return x + y; + } } internal class VectorBuilder : VectorBuilder @@ -176,8 +194,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return new DenseMatrix(storage); } - public override Matrix Sparse(SparseCompressedRowMatrixStorage storage) + public override Matrix Sparse(SparseCompressedRowMatrixStorage storage, bool cleanup = false) { + if (cleanup) + { + SumDuplicates(storage); + } return new SparseMatrix(storage); } @@ -201,6 +223,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex new ResidualStopCriterion(1e-12) }; } + + internal override Complex AddEntries(Complex x, Complex y) + { + return x + y; + } } internal class VectorBuilder : VectorBuilder @@ -239,8 +266,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return new DenseMatrix(storage); } - public override Matrix Sparse(SparseCompressedRowMatrixStorage storage) + public override Matrix Sparse(SparseCompressedRowMatrixStorage storage, bool cleanup = false) { + if (cleanup) + { + SumDuplicates(storage); + } return new SparseMatrix(storage); } @@ -264,6 +295,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 new ResidualStopCriterion(1e-6) }; } + + internal override Numerics.Complex32 AddEntries(Numerics.Complex32 x, Numerics.Complex32 y) + { + return x + y; + } } internal class VectorBuilder : VectorBuilder @@ -787,7 +823,9 @@ namespace MathNet.Numerics.LinearAlgebra /// Intended for advanced scenarios where you're working directly with /// storage for performance or interop reasons. /// - public abstract Matrix Sparse(SparseCompressedRowMatrixStorage storage); + /// The SparseCompressedRowMatrixStorage + /// Remove and sum duplicate entries. + public abstract Matrix Sparse(SparseCompressedRowMatrixStorage storage, bool cleanup = false); /// /// Create a sparse matrix of T with the given number of rows and columns. @@ -1136,22 +1174,22 @@ namespace MathNet.Numerics.LinearAlgebra // [ 0 0 f j l n ] // [ 0 d g k m 0 ] // - // Rows = 4, Columns = 6, NonZeroCount = 14 + // rows = 4, columns = 6, valueCount = 14 // - // (1) COO, Coordinate Format: + // (1) COO, Coordinate, ijv, or triplet format: // cooRowIndices = { 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 } // cooColumnIndices = { 1, 3, 0, 1, 2, 3, 2, 3, 4, 5, 1, 2, 3, 4 } // cooValues = { b, h, a, c, e, i, f, j, l, n, d, g, k, m } // - // (2) CSR, Compressed Sparse Row representation: + // (2) CSR, Compressed Sparse Row or Compressed Row Storage(CRS) or Yale format: // csrRowPointers = { 0, 2, 6, 10, 14 } // csrColumnIndices = { 1, 3, 0, 1, 2, 3, 2, 3, 4, 5, 1, 2, 3, 4 } // csrValues = { b, h, a, c, e, i, f, j, l, n, d, g, k, m } // - // (3) CSC, Compressed Sparse Column representation: - // csrColumnPointers = { 0, 1, 4, 7, 11, 13, 14 } - // csrRowIndices = { 1, 0, 1, 3, 1, 2, 3, 0, 1, 2, 3, 2, 3, 2 } - // csrValues = { a, b, c, d, e, f, g, h, i, j, k, l, m, n } + // (3) CSC, Compressed Sparse Column or Compressed Column Storage(CCS) format: + // cscColumnPointers = { 0, 1, 4, 7, 11, 13, 14 } + // cscRowIndices = { 1, 0, 1, 3, 1, 2, 3, 0, 1, 2, 3, 2, 3, 2 } + // cscValues = { a, b, c, d, e, f, g, h, i, j, k, l, m, n } /// @@ -1159,81 +1197,110 @@ namespace MathNet.Numerics.LinearAlgebra /// This new matrix will be independent from the given arrays. /// A new memory block will be allocated for storing the matrix. /// - public Matrix SparseFromCoordinateFormat(int rows, int columns, int nonZeroCount, int[] cooRowIndices, int[] cooColumnIndices, T[] cooValues) - { - if (cooValues == null) - throw new NullReferenceException(nameof(cooValues)); - if (cooRowIndices == null) - throw new NullReferenceException(nameof(cooRowIndices)); - if (cooColumnIndices == null) - throw new NullReferenceException(nameof(cooColumnIndices)); - - if (cooRowIndices.Length < nonZeroCount || cooColumnIndices.Length < nonZeroCount || cooValues.Length < nonZeroCount) + /// The number of rows. + /// The number of columns. + /// The number of stored values including explicit zeros. + /// The row index array of the coordinate format. + /// The column index array of the coordinate format. + /// The data array of the coordinate format. + /// The sparse matrix from the coordinate format. + /// Duplicate entries will be summed together and + /// explicit zeros will be not intentionally removed. + public Matrix SparseFromCoordinateFormat(int rows, int columns, int valueCount, int[] rowIndices, int[] columnIndices, T[] values) + { + if (values == null) + throw new NullReferenceException(nameof(values)); + if (rowIndices == null) + throw new NullReferenceException(nameof(rowIndices)); + if (columnIndices == null) + throw new NullReferenceException(nameof(columnIndices)); + + if (rowIndices.Length < valueCount || columnIndices.Length < valueCount || values.Length < valueCount) { - throw new Exception($"The given array has the wrong length. Should be {nonZeroCount}."); + throw new Exception($"The given array has the wrong length. Should be {valueCount}."); } // convert from COO to CSR - var csrValues = new T[nonZeroCount]; - var csrColumnIndices = new int[nonZeroCount]; + var csrValues = new T[valueCount]; + var csrColumnIndices = new int[valueCount]; var csrRowPointers = new int[rows + 1]; - for (int i = 0; i < nonZeroCount; i++) + for (int i = 0; i < valueCount; i++) { - csrRowPointers[cooRowIndices[i] + 1]++; + csrRowPointers[rowIndices[i]]++; } - for (int i = 1; i < rows + 1; i++) + for (int i = 0, cumsum = 0; i < rows; i++) { - csrRowPointers[i] += csrRowPointers[i - 1]; + var temp = csrRowPointers[i]; + csrRowPointers[i] = cumsum; + cumsum += temp; } - var curr = new int[rows]; - for (int i = 0; i < nonZeroCount; i++) + csrRowPointers[rows] = valueCount; + + for (int i = 0; i < valueCount; i++) { - int row = cooRowIndices[i]; - var loc = csrRowPointers[row] + curr[row]; - curr[row]++; + var row = rowIndices[i]; + var loc = csrRowPointers[row]; - csrColumnIndices[loc] = cooColumnIndices[i]; - csrValues[loc] = cooValues[i]; + csrColumnIndices[loc] = columnIndices[i]; + csrValues[loc] = values[i]; + + csrRowPointers[row]++; + } + for (int i = 0, last = 0; i <= rows; i++) + { + var temp = csrRowPointers[i]; + csrRowPointers[i] = last; + last = temp; } var storage = new SparseCompressedRowMatrixStorage(rows, columns, csrRowPointers, csrColumnIndices, csrValues); - return Sparse(storage); + return Sparse(storage, true); } - /// /// Create a new sparse matrix from a compressed sparse row format. /// This new matrix will be independent from the given arrays. /// A new memory block will be allocated for storing the matrix. /// - public Matrix SparseFromCompressedSparseRowFormat(int rows, int columns, int nonZeroCount, int[] csrRowPointers, int[] csrColumnIndices, T[] csrValues) - { - if (csrValues == null) - throw new NullReferenceException(nameof(csrValues)); - if (csrColumnIndices == null) - throw new NullReferenceException(nameof(csrColumnIndices)); - if (csrRowPointers == null) - throw new NullReferenceException(nameof(csrRowPointers)); - if (csrRowPointers.Length < rows) + /// The number of rows. + /// The number of columns. + /// The number of stored values including explicit zeros. + /// The row pointer array of the compressed sparse row format. + /// The column index array of the compressed sparse row format. + /// The data array of the compressed sparse row format. + /// The sparse matrix from the compressed sparse row format. + /// Duplicate entries will be summed together and + /// explicit zeros will be not intentionally removed. + public Matrix SparseFromCompressedSparseRowFormat(int rows, int columns, int valueCount, int[] rowPointers, int[] columnIndices, T[] values) + { + if (values == null) + throw new NullReferenceException(nameof(values)); + if (columnIndices == null) + throw new NullReferenceException(nameof(columnIndices)); + if (rowPointers == null) + throw new NullReferenceException(nameof(rowPointers)); + if (rowPointers.Length < rows) { throw new Exception($"The given array has the wrong length. Should be {rows + 1}."); } - if (nonZeroCount != csrRowPointers[rows]) + if (valueCount != rowPointers[rows]) { - throw new Exception($"{nameof(nonZeroCount)} should be same to {csrRowPointers[rows]}"); + throw new Exception($"{nameof(valueCount)} should be same to {rowPointers[rows]}"); } - var values = new T[nonZeroCount]; - Array.Copy(csrValues, values, nonZeroCount); - var columnIndices = new int[nonZeroCount]; - Array.Copy(csrColumnIndices, columnIndices, nonZeroCount); - var rowPointers = new int[rows + 1]; - Array.Copy(csrRowPointers, rowPointers, rows + 1); + // copy arrays to new memory block. + + var csrValues = new T[valueCount]; + Array.Copy(values, csrValues, valueCount); + var csrColumnIndices = new int[valueCount]; + Array.Copy(columnIndices, csrColumnIndices, valueCount); + var csrRowPointers = new int[rows + 1]; + Array.Copy(rowPointers, csrRowPointers, rows + 1); - var storage = new SparseCompressedRowMatrixStorage(rows, columns, rowPointers, columnIndices, values); - return Sparse(storage); + var storage = new SparseCompressedRowMatrixStorage(rows, columns, csrRowPointers, csrColumnIndices, csrValues); + return Sparse(storage, true); } /// @@ -1241,34 +1308,43 @@ namespace MathNet.Numerics.LinearAlgebra /// This new matrix will be independent from the given arrays. /// A new memory block will be allocated for storing the matrix. /// - public Matrix SparseFromCompressedSparseColumnFormat(int rows, int columns, int nonZeroCount, int[] cscRowIndices, int[] cscColumnPointers, T[] cscValues) - { - if (cscValues == null) - throw new NullReferenceException(nameof(cscValues)); - if (cscRowIndices == null) - throw new NullReferenceException(nameof(cscRowIndices)); - if (cscColumnPointers == null) - throw new NullReferenceException(nameof(cscColumnPointers)); - if (cscColumnPointers.Length < columns) + /// The number of rows. + /// The number of columns. + /// The number of stored values including explicit zeros. + /// The row index array of the compressed sparse column format. + /// The column pointer array of the compressed sparse column format. + /// The data array of the compressed sparse column format. + /// The sparse matrix from the compressed sparse column format. + /// Duplicate entries will be summed together and + /// explicit zeros will be not intentionally removed. + public Matrix SparseFromCompressedSparseColumnFormat(int rows, int columns, int valueCount, int[] rowIndices, int[] columnPointers, T[] values) + { + if (values == null) + throw new NullReferenceException(nameof(values)); + if (rowIndices == null) + throw new NullReferenceException(nameof(rowIndices)); + if (columnPointers == null) + throw new NullReferenceException(nameof(columnPointers)); + if (columnPointers.Length < columns) { throw new Exception($"The given array has the wrong length. Should be {columns + 1}."); } - if (nonZeroCount != cscColumnPointers[columns]) + if (valueCount != columnPointers[columns]) { - throw new Exception($"{nameof(nonZeroCount)} should be same to {cscColumnPointers[columns]}"); + throw new Exception($"{nameof(valueCount)} should be same to {columnPointers[columns]}"); } // convert from CSC to CSR - var csrValues = new T[nonZeroCount]; + var csrValues = new T[valueCount]; var csrRowPointers = new int[rows + 1]; - var csrColumnIndices = new int[nonZeroCount]; + var csrColumnIndices = new int[valueCount]; for (int i = 0; i < columns; i++) { - for (int j = cscColumnPointers[i]; j < cscColumnPointers[i + 1]; j++) + for (int j = columnPointers[i]; j < columnPointers[i + 1]; j++) { - csrRowPointers[cscRowIndices[j] + 1]++; + csrRowPointers[rowIndices[j] + 1]++; } } for (int i = 1; i < rows + 1; i++) @@ -1278,19 +1354,58 @@ namespace MathNet.Numerics.LinearAlgebra var curr = new int[rows]; for (int i = 0; i < columns; i++) { - for (int j = cscColumnPointers[i]; j < cscColumnPointers[i + 1]; j++) + for (int j = columnPointers[i]; j < columnPointers[i + 1]; j++) { - var loc = csrRowPointers[cscRowIndices[j]] + curr[cscRowIndices[j]]; - curr[cscRowIndices[j]]++; + var loc = csrRowPointers[rowIndices[j]] + curr[rowIndices[j]]; + curr[rowIndices[j]]++; csrColumnIndices[loc] = i; - csrValues[loc] = cscValues[j]; + csrValues[loc] = values[j]; } } var storage = new SparseCompressedRowMatrixStorage(rows, columns, csrRowPointers, csrColumnIndices, csrValues); - return Sparse(storage); + return Sparse(storage, true); + } + + // Eliminate duplicate entries by adding them together. + internal void SumDuplicates(SparseCompressedRowMatrixStorage storage) + { + int valueCount = 0; + for (int i = 0; i < storage.RowCount; i++) + { + int index = storage.RowPointers[i]; + int last = storage.RowPointers[i + 1]; + while (index < last) + { + var col = storage.ColumnIndices[index]; + var val = storage.Values[index]; + index++; + while (index < last) + { + if (storage.ColumnIndices[index] == col) + { + val = AddEntries(val, storage.Values[index]); + index++; + } + else + { + break; + } + } + storage.ColumnIndices[valueCount] = col; + storage.Values[valueCount] = val; + valueCount++; + } + storage.RowPointers[i + 1] = valueCount; + } + + // Remove extra space from arrays. + Array.Resize(ref storage.Values, valueCount); + Array.Resize(ref storage.ColumnIndices, valueCount); } + internal abstract T AddEntries(T x, T y); + /// /// Create a new diagonal matrix straight from an initialized matrix storage instance. /// The storage is used directly without copying. diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs index 85cea73f..41fe9fe8 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -84,6 +84,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage RowPointers = rowPointers; ColumnIndices = columnIndices; Values = values; + + // Explicit zeros are not intentionally removed. + // Sort ColumnIndices. + NormalizeOrdering(); } /// @@ -288,7 +292,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { MapInplace(x => x, Zeros.AllowSkip); } - + /// /// Fill zeros explicitly on the diagonal entries as required by the Intel MKL direct sparse solver. ///