Merge pull request #912 from jkalias/sparse-pointwise-multiplication-division

Pointwise multiplication and division on sparse matrices

Author: cdrnet

src/Numerics.Tests/LinearAlgebraTests/MatrixTests.cs

@@ -28,7 +28,6 @@
 // </copyright>
 
 using System;
-using System.Linq;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Double;
 using MathNet.Numerics.LinearAlgebra.Storage;
@@ -97,5 +96,25 @@ public void SparseCompressedRowMatrixStorage_CoordinateFormatDuplicateRemovalChe
                 Assert.True(Math.Abs(actual - expected) < tol, $"Expected {expected:E6} at ({r}, {c}), but got {actual:E6}");
             }
         }
+
+        [Test]
+        public void PointwiseMultiplication_SparseReturnsSameResultAsDenseTest()
+        {
+            var x = SparseMatrix.OfDiagonalArray(new double[] { 1, 2, 3, 4 });
+            var y = DenseMatrix.OfDiagonalArray(new double[] { 5, -6, 7, -8 });
+            x.PointwiseMultiply(y, x);
+            var result = SparseMatrix.OfDiagonalArray(new double[] { 5, -12, 21, -32 });
+            Assert.AreEqual(result, x);
+        }
+
+        [Test]
+        public void PointwiseDivision_SparseReturnsSameResultAsDenseTest()
+        {
+            var x = SparseMatrix.OfDiagonalArray(new double[] { 30, 20, 10 });
+            var y = DenseMatrix.OfDiagonalArray(new double[] { 3, 4, 5 });
+            x.PointwiseDivide(y, x);
+            var result = SparseMatrix.OfDiagonalArray(new double[] { 10, 5, 2 });
+            Assert.AreEqual(result, x);
+        }
     }
 }

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -1176,7 +1176,10 @@ protected override void DoTransposeThisAndMultiply(Vector<double> rightSide, Vec
         /// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
         protected override void DoPointwiseMultiply(Matrix<double> other, Matrix<double> result)
         {
-            result.Clear();
+            if (!ReferenceEquals(this, result))
+            {
+                result.Clear();
+            }
 
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;
@@ -1203,7 +1206,10 @@ protected override void DoPointwiseMultiply(Matrix<double> other, Matrix<double>
         /// <param name="result">The matrix to store the result of the pointwise division.</param>
         protected override void DoPointwiseDivide(Matrix<double> divisor, Matrix<double> result)
         {
-            result.Clear();
+            if (!ReferenceEquals(this, result))
+            {
+                result.Clear();
+            }
 
             var rowPointers = _storage.RowPointers;
             var columnIndices = _storage.ColumnIndices;