Pointwise multiplication and division on sparse matrices

The result of pointwise multiplication and division on sparse matrices is passed as an argument. The result is cleared without checking if it is equal to the current instance, thus causing a full deletion of the sparse matrix contents.

Author: jkalias

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;