LA: optimize sparse*sparse and sparse*diagonal matrix products (via wo80); cleanup

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -435,7 +435,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 return;
             }
 
-            base.DoMultiply(other, result);
+            // TODO: Map is rather generic, we can do better
+            other.MapIndexed((i, j, x) => _data[i]*x, result, false);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -923,6 +923,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Matrix<Complex> other, Matrix<Complex> result)
         {
+            var sparseOther = other as SparseMatrix;
+            var sparseResult = result as SparseMatrix;
+            if (sparseOther != null && sparseResult != null)
+            {
+                DoMultiplySparse(sparseOther, sparseResult);
+                return;
+            }
+
+            var diagonalOther = other as DiagonalMatrix;
+            if (diagonalOther != null && sparseResult != null)
+            {
+                var diagonal = ((DiagonalMatrixStorage<Complex>)other.Storage).Data;
+                // TODO: Map is rather generic, we can do better
+                MapIndexed((i, j, x) => x*diagonal[j], result, false);
+                return;
+            }
+
             result.Clear();
             var columnVector = new DenseVector(other.RowCount);
 
@@ -956,6 +973,92 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        void DoMultiplySparse(SparseMatrix other, SparseMatrix result)
+        {
+            result.Clear();
+
+            var ax = _storage.Values;
+            var ap = _storage.RowPointers;
+            var ai = _storage.ColumnIndices;
+
+            var bx = other._storage.Values;
+            var bp = other._storage.RowPointers;
+            var bi = other._storage.ColumnIndices;
+
+            int rows = RowCount;
+            int cols = other.ColumnCount;
+
+            int[] cp = result._storage.RowPointers;
+
+            var marker = new int[cols];
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            int count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                // For each row of A
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    // Row number to be added
+                    int a = ai[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        if (marker[b] != i)
+                        {
+                            marker[b] = i;
+                            count++;
+                        }
+                    }
+                }
+
+                // Record non-zero count.
+                cp[i + 1] = count;
+            }
+
+            var ci = new int[count];
+            var cx = new Complex[count];
+
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                int rowStart = cp[i];
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    int a = ai[j];
+                    Complex aEntry = ax[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        Complex bEntry = bx[k];
+                        if (marker[b] < rowStart)
+                        {
+                            marker[b] = count;
+                            ci[marker[b]] = b;
+                            cx[marker[b]] = aEntry * bEntry;
+                            count++;
+                        }
+                        else
+                        {
+                            cx[marker[b]] += aEntry * bEntry;
+                        }
+                    }
+                }
+            }
+
+            result._storage.Values = cx;
+            result._storage.ColumnIndices = ci;
+            result._storage.Normalize();
+        }
+
         /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -429,7 +429,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 return;
             }
 
-            base.DoMultiply(other, result);
+            // TODO: Map is rather generic, we can do better
+            other.MapIndexed((i, j, x) => _data[i]*x, result, false);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -917,6 +917,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Matrix<Complex32> other, Matrix<Complex32> result)
         {
+            var sparseOther = other as SparseMatrix;
+            var sparseResult = result as SparseMatrix;
+            if (sparseOther != null && sparseResult != null)
+            {
+                DoMultiplySparse(sparseOther, sparseResult);
+                return;
+            }
+
+            var diagonalOther = other as DiagonalMatrix;
+            if (diagonalOther != null && sparseResult != null)
+            {
+                var diagonal = ((DiagonalMatrixStorage<Complex32>)other.Storage).Data;
+                // TODO: Map is rather generic, we can do better
+                MapIndexed((i, j, x) => x*diagonal[j], result, false);
+                return;
+            }
+
             result.Clear();
             var columnVector = new DenseVector(other.RowCount);
 
@@ -950,6 +967,92 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        void DoMultiplySparse(SparseMatrix other, SparseMatrix result)
+        {
+            result.Clear();
+
+            var ax = _storage.Values;
+            var ap = _storage.RowPointers;
+            var ai = _storage.ColumnIndices;
+
+            var bx = other._storage.Values;
+            var bp = other._storage.RowPointers;
+            var bi = other._storage.ColumnIndices;
+
+            int rows = RowCount;
+            int cols = other.ColumnCount;
+
+            int[] cp = result._storage.RowPointers;
+
+            var marker = new int[cols];
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            int count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                // For each row of A
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    // Row number to be added
+                    int a = ai[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        if (marker[b] != i)
+                        {
+                            marker[b] = i;
+                            count++;
+                        }
+                    }
+                }
+
+                // Record non-zero count.
+                cp[i + 1] = count;
+            }
+
+            var ci = new int[count];
+            var cx = new Complex32[count];
+
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                int rowStart = cp[i];
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    int a = ai[j];
+                    Complex32 aEntry = ax[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        Complex32 bEntry = bx[k];
+                        if (marker[b] < rowStart)
+                        {
+                            marker[b] = count;
+                            ci[marker[b]] = b;
+                            cx[marker[b]] = aEntry * bEntry;
+                            count++;
+                        }
+                        else
+                        {
+                            cx[marker[b]] += aEntry * bEntry;
+                        }
+                    }
+                }
+            }
+
+            result._storage.Values = cx;
+            result._storage.ColumnIndices = ci;
+            result._storage.Normalize();
+        }
+
         /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -409,7 +409,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            base.DoMultiply(other, result);
+            // TODO: Map is rather generic, we can do better
+            other.MapIndexed((i, j, x) => _data[i]*x, result, false);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -919,6 +919,23 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Matrix<double> other, Matrix<double> result)
         {
+            var sparseOther = other as SparseMatrix;
+            var sparseResult = result as SparseMatrix;
+            if (sparseOther != null && sparseResult != null)
+            {
+                DoMultiplySparse(sparseOther, sparseResult);
+                return;
+            }
+
+            var diagonalOther = other as DiagonalMatrix;
+            if (diagonalOther != null && sparseResult != null)
+            {
+                var diagonal = ((DiagonalMatrixStorage<double>)other.Storage).Data;
+                // TODO: Map is rather generic, we can do better
+                MapIndexed((i, j, x) => x*diagonal[j], result, false);
+                return;
+            }
+
             result.Clear();
             var columnVector = new DenseVector(other.RowCount);
 
@@ -952,6 +969,92 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        void DoMultiplySparse(SparseMatrix other, SparseMatrix result)
+        {
+            result.Clear();
+
+            var ax = _storage.Values;
+            var ap = _storage.RowPointers;
+            var ai = _storage.ColumnIndices;
+
+            var bx = other._storage.Values;
+            var bp = other._storage.RowPointers;
+            var bi = other._storage.ColumnIndices;
+
+            int rows = RowCount;
+            int cols = other.ColumnCount;
+
+            int[] cp = result._storage.RowPointers;
+
+            var marker = new int[cols];
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            int count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                // For each row of A
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    // Row number to be added
+                    int a = ai[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        if (marker[b] != i)
+                        {
+                            marker[b] = i;
+                            count++;
+                        }
+                    }
+                }
+
+                // Record non-zero count.
+                cp[i + 1] = count;
+            }
+
+            var ci = new int[count];
+            var cx = new double[count];
+
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                int rowStart = cp[i];
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    int a = ai[j];
+                    double aEntry = ax[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        double bEntry = bx[k];
+                        if (marker[b] < rowStart)
+                        {
+                            marker[b] = count;
+                            ci[marker[b]] = b;
+                            cx[marker[b]] = aEntry * bEntry;
+                            count++;
+                        }
+                        else
+                        {
+                            cx[marker[b]] += aEntry * bEntry;
+                        }
+                    }
+                }
+            }
+
+            result._storage.Values = cx;
+            result._storage.ColumnIndices = ci;
+            result._storage.Normalize();
+        }
+
         /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -409,7 +409,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 return;
             }
 
-            base.DoMultiply(other, result);
+            // TODO: Map is rather generic, we can do better
+            other.MapIndexed((i, j, x) => _data[i]*x, result, false);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -924,6 +924,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <param name="result">The result of the multiplication.</param>
         protected override void DoMultiply(Matrix<float> other, Matrix<float> result)
         {
+            var sparseOther = other as SparseMatrix;
+            var sparseResult = result as SparseMatrix;
+            if (sparseOther != null && sparseResult != null)
+            {
+                DoMultiplySparse(sparseOther, sparseResult);
+                return;
+            }
+
+            var diagonalOther = other as DiagonalMatrix;
+            if (diagonalOther != null && sparseResult != null)
+            {
+                var diagonal = ((DiagonalMatrixStorage<float>)other.Storage).Data;
+                // TODO: Map is rather generic, we can do better
+                MapIndexed((i, j, x) => x*diagonal[j], result, false);
+                return;
+            }
+
             result.Clear();
             var columnVector = new DenseVector(other.RowCount);
 
@@ -957,6 +974,92 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        void DoMultiplySparse(SparseMatrix other, SparseMatrix result)
+        {
+            result.Clear();
+
+            var ax = _storage.Values;
+            var ap = _storage.RowPointers;
+            var ai = _storage.ColumnIndices;
+
+            var bx = other._storage.Values;
+            var bp = other._storage.RowPointers;
+            var bi = other._storage.ColumnIndices;
+
+            int rows = RowCount;
+            int cols = other.ColumnCount;
+
+            int[] cp = result._storage.RowPointers;
+
+            var marker = new int[cols];
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            int count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                // For each row of A
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    // Row number to be added
+                    int a = ai[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        if (marker[b] != i)
+                        {
+                            marker[b] = i;
+                            count++;
+                        }
+                    }
+                }
+
+                // Record non-zero count.
+                cp[i + 1] = count;
+            }
+
+            var ci = new int[count];
+            var cx = new float[count];
+
+            for (int ib = 0; ib < cols; ib++)
+            {
+                marker[ib] = -1;
+            }
+
+            count = 0;
+            for (int i = 0; i < rows; i++)
+            {
+                int rowStart = cp[i];
+                for (int j = ap[i]; j < ap[i + 1]; j++)
+                {
+                    int a = ai[j];
+                    float aEntry = ax[j];
+                    for (int k = bp[a]; k < bp[a + 1]; k++)
+                    {
+                        int b = bi[k];
+                        float bEntry = bx[k];
+                        if (marker[b] < rowStart)
+                        {
+                            marker[b] = count;
+                            ci[marker[b]] = b;
+                            cx[marker[b]] = aEntry * bEntry;
+                            count++;
+                        }
+                        else
+                        {
+                            cx[marker[b]] += aEntry * bEntry;
+                        }
+                    }
+                }
+            }
+
+            result._storage.Values = cx;
+            result._storage.ColumnIndices = ci;
+            result._storage.Normalize();
+        }
+
         /// <summary>
         /// Multiplies this matrix with a vector and places the results into the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         internal DenseColumnMajorMatrixStorage(int rows, int columns)
             : base(rows, columns)
         {
-            Data = new T[rows * columns];
+            Data = new T[rows*columns];
         }
 
         internal DenseColumnMajorMatrixStorage(int rows, int columns, T[] data)
@@ -58,9 +58,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new ArgumentNullException("data");
             }
 
-            if (data.Length != rows * columns)
+            if (data.Length != rows*columns)
             {
-                throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows * columns));
+                throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows*columns));
             }
 
             Data = data;
@@ -97,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// </summary>
         public override T At(int row, int column)
         {
-            return Data[(column * RowCount) + row];
+            return Data[(column*RowCount) + row];
         }
 
         /// <summary>
@@ -105,7 +105,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// </summary>
         public override void At(int row, int column, T value)
         {
-            Data[(column * RowCount) + row] = value;
+            Data[(column*RowCount) + row] = value;
         }
 
         public override void Clear()
@@ -181,7 +181,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             int columns = data.Length;
             int rows = data[0].Length;
-            var array = new T[rows * columns];
+            var array = new T[rows*columns];
             for (int j = 0; j < data.Length; j++)
             {
                 Array.Copy(data[j], 0, array, j*rows, rows);
@@ -193,10 +193,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             int rows = data.Length;
             int columns = data[0].Length;
-            var array = new T[rows * columns];
+            var array = new T[rows*columns];
             for (int j = 0; j < columns; j++)
             {
-                int offset = j * rows;
+                int offset = j*rows;
                 for (int i = 0; i < rows; i++)
                 {
                     array[offset + i] = data[i][j];
@@ -209,19 +209,19 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             int columns = data.Length;
             int rows = data[0].Length;
-            var array = new T[rows * columns];
+            var array = new T[rows*columns];
             for (int j = 0; j < data.Length; j++)
             {
                 var column = data[j];
                 var denseColumn = column as DenseVectorStorage<T>;
                 if (denseColumn != null)
                 {
-                    Array.Copy(denseColumn.Data, 0, array, j * rows, rows);
+                    Array.Copy(denseColumn.Data, 0, array, j*rows, rows);
                 }
                 else
                 {
                     // FALL BACK
-                    int offset = j * rows;
+                    int offset = j*rows;
                     for (int i = 0; i < rows; i++)
                     {
                         array[offset + i] = column.At(i);
@@ -235,10 +235,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         {
             int rows = data.Length;
             int columns = data[0].Length;
-            var array = new T[rows * columns];
+            var array = new T[rows*columns];
             for (int j = 0; j < columns; j++)
             {
-                int offset = j * rows;
+                int offset = j*rows;
                 for (int i = 0; i < rows; i++)
                 {
                     array[offset + i] = data[i].At(j);
@@ -249,10 +249,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         public static DenseColumnMajorMatrixStorage<T> OfIndexedEnumerable(int rows, int columns, IEnumerable<Tuple<int, int, T>> data)
         {
-            var array = new T[rows * columns];
+            var array = new T[rows*columns];
             foreach (var item in data)
             {
-                array[(item.Item2 * rows) + item.Item1] = item.Item3;
+                array[(item.Item2*rows) + item.Item1] = item.Item3;
             }
             return new DenseColumnMajorMatrixStorage<T>(rows, columns, array);
         }
@@ -388,9 +388,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             var targetDense = target as DenseVectorStorage<T>;
             if (targetDense != null)
             {
-                for (int j = 0; j<columnCount; j++)
+                for (int j = 0; j < columnCount; j++)
                 {
-                    targetDense.Data[j + targetColumnIndex] = Data[(j + sourceColumnIndex) * RowCount + rowIndex];
+                    targetDense.Data[j + targetColumnIndex] = Data[(j + sourceColumnIndex)*RowCount + rowIndex];
                 }
                 return;
             }
@@ -399,7 +399,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
             {
-                target.At(jj, Data[(j * RowCount) + rowIndex]);
+                target.At(jj, Data[(j*RowCount) + rowIndex]);
             }
         }
 
@@ -416,7 +416,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             // FALL BACK
 
-            var offset = columnIndex * RowCount;
+            var offset = columnIndex*RowCount;
             for (int i = sourceRowIndex, ii = targetRowIndex; i < sourceRowIndex + rowCount; i++, ii++)
             {
                 target.At(ii, Data[offset + i]);
@@ -430,10 +430,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             var ret = new T[Data.Length];
             for (int i = 0; i < RowCount; i++)
             {
-                var offset = i * ColumnCount;
+                var offset = i*ColumnCount;
                 for (int j = 0; j < ColumnCount; j++)
                 {
-                    ret[offset + j] = Data[(j * RowCount) + i];
+                    ret[offset + j] = Data[(j*RowCount) + i];
                 }
             }
             return ret;
@@ -453,7 +453,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (int j = 0; j < ColumnCount; j++)
                 {
-                    ret[i, j] = Data[(j * RowCount) + i];
+                    ret[i, j] = Data[(j*RowCount) + i];
                 }
             }
             return ret;

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -182,7 +182,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (var i = 0; i < hashNum; i++)
                 {
-                    hash = hash * 31 + Data[i].GetHashCode();
+                    hash = hash*31 + Data[i].GetHashCode();
                 }
             }
             return hash;
@@ -253,7 +253,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
 
             var storage = new DiagonalMatrixStorage<T>(rows, columns);
-            foreach(var item in data)
+            foreach (var item in data)
             {
                 storage.Data[item.Item1] = item.Item2;
             }
@@ -400,7 +400,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 // column by column, but skip resulting zero columns at the beginning
 
                 int columnInit = sourceRowIndex - sourceColumnIndex;
-                int offset = (columnInit + targetColumnIndex) * target.RowCount + targetRowIndex;
+                int offset = (columnInit + targetColumnIndex)*target.RowCount + targetRowIndex;
                 int step = target.RowCount + 1;
                 int end = Math.Min(columnCount - columnInit, rowCount) + sourceRowIndex;
 
@@ -513,22 +513,22 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         public override T[] ToRowMajorArray()
         {
-            var ret = new T[RowCount * ColumnCount];
+            var ret = new T[RowCount*ColumnCount];
             var stride = ColumnCount + 1;
             for (int i = 0; i < Data.Length; i++)
             {
-                ret[i * stride] = Data[i];
+                ret[i*stride] = Data[i];
             }
             return ret;
         }
 
         public override T[] ToColumnMajorArray()
         {
-            var ret = new T[RowCount * ColumnCount];
+            var ret = new T[RowCount*ColumnCount];
             var stride = RowCount + 1;
             for (int i = 0; i < Data.Length; i++)
             {
-                ret[i * stride] = Data[i];
+                ret[i*stride] = Data[i];
             }
             return ret;
         }

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -30,7 +30,6 @@
 
 using System;
 using System.Collections.Generic;
-using System.Diagnostics;
 using System.Linq;
 using MathNet.Numerics.Properties;
 
@@ -160,18 +159,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 var valueCount = RowPointers[RowPointers.Length - 1];
 
                 // Check if the storage needs to be increased
-                if ((valueCount == Values.Length) && (valueCount < ((long)RowCount * ColumnCount)))
+                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);
+                    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 Values, (int)size);
-                    Array.Resize(ref ColumnIndices, (int) size);
+                    Array.Resize(ref ColumnIndices, (int)size);
                 }
 
                 // Move all values (with a position larger than index) in the value array to the next position
@@ -217,7 +216,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             // Check whether we need to shrink the arrays. 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))
+            if ((valueCount > 1024) && (valueCount < Values.Length/2))
             {
                 Array.Resize(ref Values, valueCount);
                 Array.Resize(ref ColumnIndices, valueCount);
@@ -264,6 +263,20 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return delta;
         }
 
+        public void Normalize()
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                int index = RowPointers[i];
+                int count = RowPointers[i + 1] - index;
+                if (count > 1)
+                {
+                    Sorting.Sort(ColumnIndices, Values, index, count);
+                }
+            }
+            MapInplace(x => x, forceMapZeros: false);
+        }
+
         public override void Clear()
         {
             Array.Clear(RowPointers, 0, RowPointers.Length);
@@ -316,7 +329,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             // Check whether we need to shrink the arrays. 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))
+            if ((valueCount > 1024) && (valueCount < Values.Length/2))
             {
                 Array.Resize(ref Values, valueCount);
                 Array.Resize(ref ColumnIndices, valueCount);