Fix bugs in converting from not-sorted and duplicate entries of COO format arrays to CSR format.

src/Numerics.Tests/LinearAlgebraTests/MatrixStructureTheory.cs

@@ -27,11 +27,10 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-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<T>.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<T>.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<T>.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<T>.Build.Random(5, 0).ToArray();
+
+            var A = Matrix<T>.Build.SparseFromCoordinateFormat(rowCount, columnCount, valueCount, cooRowIndices, cooColumnIndices, cooValues);
+
+            cooRowIndices.Reverse();
+            cooColumnIndices.Reverse();
+            cooValues.Reverse();
+
+            var B = Matrix<T>.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<T>.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<T>.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<T>.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<T>.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<T>.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<T>.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]

src/Numerics/LinearAlgebra/Builder.cs

@@ -48,8 +48,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DenseMatrix(storage);
         }
 
-        public override Matrix<double> Sparse(SparseCompressedRowMatrixStorage<double> storage)
+        public override Matrix<double> Sparse(SparseCompressedRowMatrixStorage<double> storage, bool cleanup = false)
         {
+            if (cleanup)
+            {
+                SumDuplicates(storage);
+            }
             return new SparseMatrix(storage);
         }
 
@@ -73,6 +77,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 new ResidualStopCriterion<double>(1e-12)
             };
         }
+
+        internal override double AddEntries(double x, double y)
+        {
+            return x + y;
+        }
     }
 
     internal class VectorBuilder : VectorBuilder<double>
@@ -111,8 +120,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DenseMatrix(storage);
         }
 
-        public override Matrix<float> Sparse(SparseCompressedRowMatrixStorage<float> storage)
+        public override Matrix<float> Sparse(SparseCompressedRowMatrixStorage<float> storage, bool cleanup = false)
         {
+            if (cleanup)
+            {
+                SumDuplicates(storage);
+            }
             return new SparseMatrix(storage);
         }
 
@@ -136,6 +149,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 new ResidualStopCriterion<float>(1e-6)
             };
         }
+
+        internal override float AddEntries(float x, float y)
+        {
+            return x + y;
+        }
     }
 
     internal class VectorBuilder : VectorBuilder<float>
@@ -176,8 +194,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DenseMatrix(storage);
         }
 
-        public override Matrix<Complex> Sparse(SparseCompressedRowMatrixStorage<Complex> storage)
+        public override Matrix<Complex> Sparse(SparseCompressedRowMatrixStorage<Complex> storage, bool cleanup = false)
         {
+            if (cleanup)
+            {
+                SumDuplicates(storage);
+            }
             return new SparseMatrix(storage);
         }
 
@@ -201,6 +223,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 new ResidualStopCriterion<Complex>(1e-12)
             };
         }
+
+        internal override Complex AddEntries(Complex x, Complex y)
+        {
+            return x + y;
+        }
     }
 
     internal class VectorBuilder : VectorBuilder<Complex>
@@ -239,8 +266,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DenseMatrix(storage);
         }
 
-        public override Matrix<Numerics.Complex32> Sparse(SparseCompressedRowMatrixStorage<Numerics.Complex32> storage)
+        public override Matrix<Numerics.Complex32> Sparse(SparseCompressedRowMatrixStorage<Numerics.Complex32> storage, bool cleanup = false)
         {
+            if (cleanup)
+            {
+                SumDuplicates(storage);
+            }
             return new SparseMatrix(storage);
         }
 
@@ -264,6 +295,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 new ResidualStopCriterion<Numerics.Complex32>(1e-6)
             };
         }
+
+        internal override Numerics.Complex32 AddEntries(Numerics.Complex32 x, Numerics.Complex32 y)
+        {
+            return x + y;
+        }
     }
 
     internal class VectorBuilder : VectorBuilder<Numerics.Complex32>
@@ -787,7 +823,9 @@ namespace MathNet.Numerics.LinearAlgebra
         /// Intended for advanced scenarios where you're working directly with
         /// storage for performance or interop reasons.
         /// </summary>
-        public abstract Matrix<T> Sparse(SparseCompressedRowMatrixStorage<T> storage);
+        /// <param name="storage">The SparseCompressedRowMatrixStorage</param>
+        /// <param name="cleanup">Remove and sum duplicate entries.</param>
+        public abstract Matrix<T> Sparse(SparseCompressedRowMatrixStorage<T> storage, bool cleanup = false);
 
         /// <summary>
         /// 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<T>(rows, columns, csrRowPointers, csrColumnIndices, csrValues);
-            return Sparse(storage);
+            return Sparse(storage, true);
         }
 
-
         /// <summary>
         /// 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<T>(rows, columns, rowPointers, columnIndices, values);
-            return Sparse(storage);
+            return Sparse(storage, true);
         }
 
         /// <summary>
@@ -1288,9 +1334,47 @@ namespace MathNet.Numerics.LinearAlgebra
             }
 
             var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns, csrRowPointers, csrColumnIndices, csrValues);
-            return Sparse(storage);
+            return Sparse(storage, true);
         }
 
+        internal void SumDuplicates(SparseCompressedRowMatrixStorage<T> 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);
+
         /// <summary>
         /// Create a new diagonal matrix straight from an initialized matrix storage instance.
         /// The storage is used directly without copying.

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();
         }
 
         /// <summary>
@@ -288,7 +291,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             MapInplace(x => x, Zeros.AllowSkip);
         }
-
+ 
         /// <summary>
         /// Fill zeros explicitly on the diagonal entries as required by the Intel MKL direct sparse solver.
         /// </summary>