LA: special case for matrix inplace map

12 years ago · ebc5471d5c
5 changed files with 202 additions and 12 deletions
--- a/src/Numerics/LinearAlgebra/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Matrix.cs
@ -1473,7 +1473,7 @@ namespace MathNet.Numerics.LinearAlgebra
        /// </summary>
        public void MapInplace(Func<T, T> f, bool forceMapZeros = false)
        {
-            Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
+            Storage.MapInplace(f, forceMapZeros);
        }

        /// <summary>
@ -1484,7 +1484,7 @@ namespace MathNet.Numerics.LinearAlgebra
        /// </summary>
        public void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
        {
-            Storage.MapIndexedToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
+            Storage.MapIndexedInplace(f, forceMapZeros);
        }

        /// <summary>
--- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
@ -530,20 +530,34 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            }
        }

+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    Data[i] = f(Data[i]);
+                }
+            });
+        }
+
        internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
        {
            var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
            if (denseTarget != null)
            {
-                int index = 0;
-                for (int j = 0; j < ColumnCount; j++)
+                CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) =>
                {
-                    for (int i = 0; i < RowCount; i++)
+                    int index = a*RowCount;
+                    for (int j = a; j < b; j++)
                    {
-                        denseTarget.Data[index] = f(i, j, Data[index]);
-                        index++;
+                        for (int i = 0; i < RowCount; i++)
+                        {
+                            denseTarget.Data[index] = f(i, j, Data[index]);
+                            index++;
+                        }
                    }
-                }
+                });
                return;
            }

@ -558,5 +572,21 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                }
            }
        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) =>
+            {
+                int index = a*RowCount;
+                for (int j = a; j < b; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        Data[index] = f(i, j, Data[index]);
+                        index++;
+                    }
+                }
+            });
+        }
    }
 }
--- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
@ -637,6 +637,22 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            }
        }

+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            if (forceMapZeros)
+            {
+                throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
+            }
+
+            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    Data[i] = f(Data[i]);
+                }
+            });
+        }
+
        internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
        {
            var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero));
@ -684,5 +700,21 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                }
            }
        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            if (forceMapZeros)
+            {
+                throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
+            }
+
+            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    Data[i] = f(i, i, Data[i]);
+                }
+            });
+        }
    }
 }
--- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
@ -512,6 +512,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            }
        }

+        public virtual void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    At(i, j, f(At(i, j)));
+                }
+            }
+        }
+
        public void MapIndexedTo<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
            where TU : struct, IEquatable<TU>, IFormattable
        {
@ -540,5 +551,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                }
            }
        }
+
+        public virtual void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    At(i, j, f(i, j, At(i, j)));
+                }
+            }
+        }
    }
 }
--- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
@ -1113,8 +1113,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            if (sparseTarget != null)
            {
                var newRowPointers = sparseTarget.RowPointers;
-                var newColumnIndices = new List<int>();
-                var newValues = new List<TU>();
+                var newColumnIndices = new List<int>(ColumnIndices.Length);
+                var newValues = new List<TU>(Values.Length);

                if (processZeros)
                {
@ -1199,6 +1199,59 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            }
        }

+        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        {
+            if (forceMapZeros || !Zero.Equals(f(Zero)))
+            {
+                var newRowPointers = RowPointers;
+                var newColumnIndices = new List<int>(ColumnIndices.Length);
+                var newValues = new List<T>(Values.Length);
+
+                int k = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    newRowPointers[row] = newValues.Count;
+                    for (int col = 0; col < ColumnCount; col++)
+                    {
+                        var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero);
+                        if (!Zero.Equals(item))
+                        {
+                            newValues.Add(item);
+                            newColumnIndices.Add(col);
+                        }
+                    }
+                }
+
+                ColumnIndices = newColumnIndices.ToArray();
+                Values = newValues.ToArray();
+                newRowPointers[RowCount] = newValues.Count;
+            }
+            else
+            {
+                // we can safely do this in-place:
+                int nonZero = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    var startIndex = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    RowPointers[row] = nonZero;
+                    for (var j = startIndex; j < endIndex; j++)
+                    {
+                        var item = f(Values[j]);
+                        if (!Zero.Equals(item))
+                        {
+                            Values[nonZero] = item;
+                            ColumnIndices[nonZero] = ColumnIndices[j];
+                            nonZero++;
+                        }
+                    }
+                }
+                Array.Resize(ref ColumnIndices, nonZero);
+                Array.Resize(ref Values, nonZero);
+                RowPointers[RowCount] = nonZero;
+            }
+        }
+
        internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
        {
            var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero));
@ -1207,8 +1260,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
            if (sparseTarget != null)
            {
                var newRowPointers = sparseTarget.RowPointers;
-                var newColumnIndices = new List<int>();
-                var newValues = new List<TU>();
+                var newColumnIndices = new List<int>(ColumnIndices.Length);
+                var newValues = new List<TU>(Values.Length);

                if (processZeros)
                {
@ -1292,5 +1345,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                }
            }
        }
+
+        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        {
+            if (forceMapZeros || !Zero.Equals(f(0, 0, Zero)))
+            {
+                var newRowPointers = RowPointers;
+                var newColumnIndices = new List<int>(ColumnIndices.Length);
+                var newValues = new List<T>(Values.Length);
+
+                int k = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    newRowPointers[row] = newValues.Count;
+                    for (int col = 0; col < ColumnCount; col++)
+                    {
+                        var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero);
+                        if (!Zero.Equals(item))
+                        {
+                            newValues.Add(item);
+                            newColumnIndices.Add(col);
+                        }
+                    }
+                }
+
+                ColumnIndices = newColumnIndices.ToArray();
+                Values = newValues.ToArray();
+                newRowPointers[RowCount] = newValues.Count;
+            }
+            else
+            {
+                // we can safely do this in-place:
+                int nonZero = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    var startIndex = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    RowPointers[row] = nonZero;
+                    for (var j = startIndex; j < endIndex; j++)
+                    {
+                        var item = f(row, ColumnIndices[j], Values[j]);
+                        if (!Zero.Equals(item))
+                        {
+                            Values[nonZero] = item;
+                            ColumnIndices[nonZero] = ColumnIndices[j];
+                            nonZero++;
+                        }
+                    }
+                }
+                Array.Resize(ref ColumnIndices, nonZero);
+                Array.Resize(ref Values, nonZero);
+                RowPointers[RowCount] = nonZero;
+            }
+        }
    }
 }