Merge pull request #409 from arthurvb/408-sparsedensemultiplication

Multiplying sparse with dense matrices with DenseColumnMajorMatrixStorage

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -911,12 +911,41 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
 
             result.Clear();
-            var columnVector = new DenseVector(other.RowCount);
+            
 
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;
             var values = _storage.Values;
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex>;
+            if (denseOther != null)
+            {
+                // in this case we can directly address the underlying data-array
+                for (var row = 0; row < RowCount; row++)
+                {
+                    var startIndex = rowPointers[row];
+                    var endIndex = rowPointers[row + 1];
+
+                    if (startIndex == endIndex)
+                    {
+                        continue;
+                    }
+
+                    for (var column = 0; column < other.ColumnCount; column++)
+                    {
+                        int otherColumnStartPosition = column * other.RowCount;
+                        var sum = Complex.Zero;
+                        for (var index = startIndex; index < endIndex; index++)
+                        {
+                            sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
+                        }
 
+                        result.At(row, column, sum);
+                    }
+                }
+                return;
+            }
+
+            var columnVector = new DenseVector(other.RowCount);
             for (var row = 0; row < RowCount; row++)
             {
                 var startIndex = rowPointers[row];

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -905,12 +905,41 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
 
             result.Clear();
-            var columnVector = new DenseVector(other.RowCount);
+            
 
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;
             var values = _storage.Values;
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex32>;
+            if (denseOther != null)
+            {
+                // in this case we can directly address the underlying data-array
+                for (var row = 0; row < RowCount; row++)
+                {
+                    var startIndex = rowPointers[row];
+                    var endIndex = rowPointers[row + 1];
+
+                    if (startIndex == endIndex)
+                    {
+                        continue;
+                    }
+
+                    for (var column = 0; column < other.ColumnCount; column++)
+                    {
+                        int otherColumnStartPosition = column * other.RowCount;
+                        var sum = Complex32.Zero;
+                        for (var index = startIndex; index < endIndex; index++)
+                        {
+                            sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
+                        }
 
+                        result.At(row, column, sum);
+                    }
+                }
+                return;
+            }
+
+            var columnVector = new DenseVector(other.RowCount);
             for (var row = 0; row < RowCount; row++)
             {
                 var startIndex = rowPointers[row];

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -906,12 +906,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             result.Clear();
-            var columnVector = new DenseVector(other.RowCount);
-
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;
             var values = _storage.Values;
 
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<double>;
+            if (denseOther != null)
+            {
+                // in this case we can directly address the underlying data-array
+                for (var row = 0; row < RowCount; row++)
+                {
+                    var startIndex = rowPointers[row];
+                    var endIndex = rowPointers[row + 1];
+
+                    if (startIndex == endIndex)
+                    {
+                        continue;
+                    }
+
+                    for (var column = 0; column < other.ColumnCount; column++)
+                    {
+                        int otherColumnStartPosition = column * other.RowCount;
+                        var sum = 0d;
+                        for (var index = startIndex; index < endIndex; index++)
+                        {
+                            sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
+                        }
+
+                        result.At(row, column, sum);
+                    }
+                }
+                return;
+            }
+
+            var columnVector = new DenseVector(other.RowCount);
             for (var row = 0; row < RowCount; row++)
             {
                 var startIndex = rowPointers[row];

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -909,12 +909,39 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             result.Clear();
-            var columnVector = new DenseVector(other.RowCount);
-
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;
             var values = _storage.Values;
+            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<float>;
+            if (denseOther != null)
+            {
+                // in this case we can directly address the underlying data-array
+                for (var row = 0; row < RowCount; row++)
+                {
+                    var startIndex = rowPointers[row];
+                    var endIndex = rowPointers[row + 1];
 
+                    if (startIndex == endIndex)
+                    {
+                        continue;
+                    }
+
+                    for (var column = 0; column < other.ColumnCount; column++)
+                    {
+                        int otherColumnStartPosition = column * other.RowCount;
+                        var sum = 0f;
+                        for (var index = startIndex; index < endIndex; index++)
+                        {
+                            sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
+                        }
+
+                        result.At(row, column, sum);
+                    }
+                }
+                return;
+            }
+
+            var columnVector = new DenseVector(other.RowCount);
             for (var row = 0; row < RowCount; row++)
             {
                 var startIndex = rowPointers[row];