From c31bcaa71db50a8d4a44c3f1efa7b4c8d221bdcb Mon Sep 17 00:00:00 2001 From: Jong Hyun Kim Date: Mon, 10 Aug 2020 18:51:41 +0900 Subject: [PATCH] Fix bugs in converting from not-sorted and duplicate entries of COO format arrays to CSR format. --- .../MatrixStructureTheory.cs | 93 +++++++++++--- src/Numerics/LinearAlgebra/Builder.cs | 118 +++++++++++++++--- .../SparseCompressedRowMatrixStorage.cs | 5 +- 3 files changed, 180 insertions(+), 36 deletions(-) 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..4b550830 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. @@ -1128,7 +1166,6 @@ namespace MathNet.Numerics.LinearAlgebra return m; } - // Representation of Sparse Matrix // // Matrix A = [ 0 b 0 h 0 0 ] @@ -1181,28 +1218,37 @@ namespace MathNet.Numerics.LinearAlgebra for (int i = 0; i < nonZeroCount; i++) { - csrRowPointers[cooRowIndices[i] + 1]++; + csrRowPointers[cooRowIndices[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]; + csrRowPointers[rows] = nonZeroCount; + for (int i = 0; i < nonZeroCount; i++) { - int row = cooRowIndices[i]; - var loc = csrRowPointers[row] + curr[row]; - curr[row]++; + var row = cooRowIndices[i]; + var loc = csrRowPointers[row]; csrColumnIndices[loc] = cooColumnIndices[i]; csrValues[loc] = cooValues[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. @@ -1233,7 +1279,7 @@ namespace MathNet.Numerics.LinearAlgebra Array.Copy(csrRowPointers, rowPointers, rows + 1); var storage = new SparseCompressedRowMatrixStorage(rows, columns, rowPointers, columnIndices, values); - return Sparse(storage); + return Sparse(storage, true); } /// @@ -1288,9 +1334,47 @@ namespace MathNet.Numerics.LinearAlgebra } var storage = new SparseCompressedRowMatrixStorage(rows, columns, csrRowPointers, csrColumnIndices, csrValues); - return Sparse(storage); + return Sparse(storage, true); } + internal void SumDuplicates(SparseCompressedRowMatrixStorage storage) + { + int nnz = 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[nnz] = col; + storage.Values[nnz] = val; + nnz++; + } + storage.RowPointers[i + 1] = nnz; + } + + // Remove extra space from arrays. + Array.Resize(ref storage.Values, nnz); + Array.Resize(ref storage.ColumnIndices, nnz); + } + + 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..7cd7ec57 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -84,6 +84,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage RowPointers = rowPointers; ColumnIndices = columnIndices; Values = values; + + // columnIndices may be not ordered. + Normalize(); } /// @@ -288,7 +291,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. ///