LA Storage: sparse matrix storage

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -6,14 +6,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     internal class DenseColumnMajorMatrixStorage<T>
         where T : struct, IEquatable<T>, 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;
         }
 

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -0,0 +1,229 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Storage
+{
+    internal class SparseCompressedRowMatrixStorage<T>
+        where T : struct, IEquatable<T>, IFormattable
+    {
+        // [ruegg] public fields are OK here
+
+        public readonly int RowCount;
+        public readonly int ColumnCount;
+        readonly T _zero;
+
+        /// <summary>
+        /// The array containing the row indices of the existing rows. Element "j" of the array gives the index of the 
+        /// element in the <see cref="Values"/> array that is first non-zero element in a row "j"
+        /// </summary>
+        public readonly int[] RowPointers;
+
+        /// <summary>
+        /// 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 <see cref="_nonZeroValues"/> array.
+        /// </summary>
+        public int[] ColumnIndices;
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        public T[] Values;
+
+        /// <summary>
+        /// Gets the number of non zero elements in the matrix.
+        /// </summary>
+        /// <value>The number of non zero elements.</value>
+        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);
+        }
+
+        /// <summary>
+        /// Retrieves the requested element without range checking.
+        /// </summary>
+        /// <param name="row">
+        /// The row of the element.
+        /// </param>
+        /// <param name="column">
+        /// The column of the element.
+        /// </param>
+        /// <returns>
+        /// The requested element.
+        /// </returns>
+        public T GetValueAt(int row, int column)
+        {
+            var index = FindItem(row, column);
+            return index >= 0 ? Values[index] : _zero;
+        }
+
+        /// <summary>
+        /// Created this method because we cannot call "virtual At" in constructor of the class, but we need to do it
+        /// </summary>
+        /// <param name="row"> The row of the element. </param>
+        /// <param name="column"> The column of the element. </param>
+        /// <param name="value"> The value to set the element to. </param>
+        /// <remarks>WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks</remarks>
+        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;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Delete value from internal storage
+        /// </summary>
+        /// <param name="itemIndex">Index of value in nonZeroValues array</param>
+        /// <param name="row">Row number of matrix</param>
+        /// <remarks>WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks</remarks>
+        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);
+            }
+        }
+
+        /// <summary>
+        /// Find item Index in nonZeroValues array
+        /// </summary>
+        /// <param name="row">Matrix row index</param>
+        /// <param name="column">Matrix column index</param>
+        /// <returns>Item index</returns>
+        /// <remarks>WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks</remarks>
+        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);
+        }
+
+        /// <summary>
+        /// Calculates the amount with which to grow the storage array's if they need to be
+        /// increased in size.
+        /// </summary>
+        /// <returns>The amount grown.</returns>
+        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;
+        }
+    }
+}

src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs

@@ -6,25 +6,30 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
     internal class SparseDiagonalMatrixStorage<T>
         where T : struct, IEquatable<T>, 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]

src/Numerics/Numerics.csproj

@@ -331,6 +331,7 @@
     <Compile Include="LinearAlgebra\Generic\Vector.cs" />
     <Compile Include="LinearAlgebra\Single\Vector.cs" />
     <Compile Include="LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs" />
+    <Compile Include="LinearAlgebra\Storage\SparseCompressedRowMatrixStorage.cs" />
     <Compile Include="LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs" />
     <Compile Include="Permutation.cs" />
     <Compile Include="Distributions\Continuous\Beta.cs" />

src/Portable/Portable.csproj

@@ -891,6 +891,9 @@
     <Compile Include="..\Numerics\LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs">
       <Link>LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Storage\SparseCompressedRowMatrixStorage.cs">
+      <Link>LinearAlgebra\Storage\SparseCompressedRowMatrixStorage.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs">
       <Link>LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs</Link>
     </Compile>