LA: storage-aware non-inplace map on matrices

src/FSharp/LinearAlgebra.Matrix.fs

@@ -207,22 +207,18 @@ module Matrix =
 
 
     /// In-place map of every matrix element using a function.
-    let inline mapInPlace f (A: #Matrix<_>) =
-        A.MapInplace((fun x -> f x), true)
-
-    /// In-place map of every matrix element using a position dependent function.
-    let inline mapiInPlace f (A: #Matrix<_>) =
-        A.MapIndexedInplace((fun i j x -> f i j x), true)
+    let inline mapInPlace f (A: #Matrix<_>) = A.MapInplace((fun x -> f x), true)
 
     /// In-place map of every matrix element using a function.
     /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapSkipZerosInPlace f (A: #Matrix<_>) =
-        A.MapInplace((fun x -> f x), false)
+    let inline mapSkipZerosInPlace f (A: #Matrix<_>) = A.MapInplace((fun x -> f x), false)
+
+    /// In-place map of every matrix element using a position dependent function.
+    let inline mapiInPlace f (A: #Matrix<_>) = A.MapIndexedInplace((fun i j x -> f i j x), true)
 
     /// In-place map of every matrix element using a position dependent function.
     /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapiSkipZerosInPlace f (A: #Matrix<_>) =
-        A.MapIndexedInplace((fun i j x -> f i j x), false)
+    let inline mapiSkipZerosInPlace f (A: #Matrix<_>) = A.MapIndexedInplace((fun i j x -> f i j x), false)
 
     /// In-place map every matrix column using the given position dependent function.
     let inline mapColsInPlace (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
@@ -236,30 +232,18 @@ module Matrix =
 
 
     /// Map every matrix element using the given function.
-    let inline map f (A: #Matrix<_>) =
-        let A = A.Clone()
-        A.MapInplace((fun x -> f x), true)
-        A
+    let inline map f (A: #Matrix<_>) = A.Map((fun x -> f x), true)
 
     /// Map every matrix element using the given function.
     /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapSkipZeros f (A: #Matrix<_>) =
-        let A = A.Clone()
-        A.MapInplace((fun x -> f x), false)
-        A
+    let inline mapSkipZeros f (A: #Matrix<_>) = A.Map((fun x -> f x), false)
 
     /// Map every matrix element using the given position dependent function.
-    let inline mapi f (A: #Matrix<_>) =
-        let A = A.Clone()
-        A.MapIndexedInplace((fun i j x -> f i j x), true)
-        A
+    let inline mapi f (A: #Matrix<_>) = A.MapIndexed((fun i j x -> f i j x), true)
 
     /// Map every matrix element using the given position dependent function.
     /// Zero-values may be skipped (relevant mostly for sparse matrices).
-    let inline mapiSkipZeros f (A: #Matrix<_>) =
-        let A = A.Clone()
-        A.MapIndexedInplace((fun i j x -> f i j x), false)
-        A
+    let inline mapiSkipZeros f (A: #Matrix<_>) = A.MapIndexed((fun i j x -> f i j x), false)
 
     /// Map every matrix column using the given position dependent function.
     let inline mapCols (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =

src/Numerics/LinearAlgebra/Matrix.cs

@@ -1455,7 +1455,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public void MapInplace(Func<T, T> f, bool forceMapZeros = false)
         {
-            Storage.MapInplace(f, forceMapZeros);
+            Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
         }
 
         /// <summary>
@@ -1466,7 +1466,57 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
         {
-            Storage.MapIndexedInplace(f, forceMapZeros);
+            Storage.MapIndexedToUnchecked(Storage, f, forceMapZeros, skipClearing: true);
+        }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and replaces the value in the result matrix.
+        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
+        /// on the actual data storage implementation (relevant mostly for sparse matrices).
+        /// </summary>
+        public void Map<TU>(Func<T, TU> f, Matrix<TU> result, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            Storage.MapTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros);
+        }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and replaces the value in the result matrix.
+        /// The index of each value (zero-based) is passed as first argument to the function.
+        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
+        /// on the actual data storage implementation (relevant mostly for sparse matrices).
+        /// </summary>
+        public void MapIndexed<TU>(Func<int, int, T, TU> f, Matrix<TU> result, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            Storage.MapIndexedTo(result.Storage, f, forceMapZeros, skipClearing: forceMapZeros);
+        }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and returns the results as a new matrix.
+        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
+        /// on the actual data storage implementation (relevant mostly for sparse matrices).
+        /// </summary>
+        public Matrix<TU> Map<TU>(Func<T, TU> f, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            var result = Matrix<TU>.Build.SameAs(this);
+            Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
+            return result;
+        }
+
+        /// <summary>
+        /// Applies a function to each value of this matrix and returns the results as a new matrix.
+        /// The index of each value (zero-based) is passed as first argument to the function.
+        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
+        /// on the actual data storage implementation (relevant mostly for sparse matrices).
+        /// </summary>
+        public Matrix<TU> MapIndexed<TU>(Func<int, int, T, TU> f, bool forceMapZeros = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            var result = Matrix<TU>.Build.SameAs(this);
+            Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
+            return result;
         }
     }
 }

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -503,26 +503,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         // FUNCTIONAL COMBINATORS
 
-        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
+            if (denseTarget != null)
+            {
+                CommonParallel.For(0, Data.Length, 4096, (a, b) =>
                 {
                     for (int i = a; i < b; i++)
                     {
-                        Data[i] = f(Data[i]);
+                        denseTarget.Data[i] = f(Data[i]);
                     }
                 });
+                return;
+            }
+
+            // FALL BACK
+
+            int index = 0;
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    target.At(i, j, f(Data[index++]));
+                }
+            }
         }
 
-        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            int index = 0;
+            var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
+            if (denseTarget != null)
+            {
+                int index = 0;
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        denseTarget.Data[index] = f(i, j, Data[index]);
+                        index++;
+                    }
+                }
+                return;
+            }
+
+            // FALL BACK
+
+            int index2 = 0;
             for (int j = 0; j < ColumnCount; j++)
             {
                 for (int i = 0; i < RowCount; i++)
                 {
-                    Data[index] = f(i, j, Data[index]);
-                    index++;
+                    target.At(i, j, f(i, j, Data[index2++]));
                 }
             }
         }

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -589,30 +589,100 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         // FUNCTIONAL COMBINATORS
 
-        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            // we deliberately ignore forceMapZeros since we would not actually
-            // support any non-zero results outside of the diagonal anyway
-            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            var processZeros = forceMapZeros || !Zero.Equals(f(Zero));
+
+            var diagonalTarget = target as DiagonalMatrixStorage<TU>;
+            if (diagonalTarget != null)
+            {
+                if (processZeros)
+                {
+                    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]);
+                        diagonalTarget.Data[i] = f(Data[i]);
                     }
                 });
+                return;
+            }
+
+            // FALL BACK
+
+            if (!skipClearing && !processZeros)
+            {
+                target.Clear();
+            }
+
+            if (processZeros)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        target.At(i, j, f(i == j ? Data[i] : Zero));
+                    }
+                }
+            }
+            else
+            {
+                for (int i = 0; i < Data.Length; i++)
+                {
+                    target.At(i, i, f(Data[i]));
+                }
+            }
         }
 
-        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            // we deliberately ignore forceMapZeros since we would not actually
-            // support any non-zero results outside of the diagonal anyway
-            CommonParallel.For(0, Data.Length, 4096, (a, b) =>
+            var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero));
+
+            var diagonalTarget = target as DiagonalMatrixStorage<TU>;
+            if (diagonalTarget != null)
+            {
+                if (processZeros)
+                {
+                    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]);
+                        diagonalTarget.Data[i] = f(i, i, Data[i]);
                     }
                 });
+                return;
+            }
+
+            // FALL BACK
+
+            if (!skipClearing && !processZeros)
+            {
+                target.Clear();
+            }
+
+            if (processZeros)
+            {
+                for (int j = 0; j < ColumnCount; j++)
+                {
+                    for (int i = 0; i < RowCount; i++)
+                    {
+                        target.At(i, j, f(i, j, i == j ? Data[i] : Zero));
+                    }
+                }
+            }
+            else
+            {
+                for (int i = 0; i < Data.Length; i++)
+                {
+                    target.At(i, i, f(i, i, Data[i]));
+                }
+            }
         }
     }
 }

src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@@ -468,24 +468,60 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         // FUNCTIONAL COMBINATORS
 
-        public virtual void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        public void MapTo<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, IFormattable
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
+            {
+                var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
+                throw new ArgumentException(message, "target");
+            }
+
+            MapToUnchecked(target, f, forceMapZeros, skipClearing);
+        }
+
+        internal virtual void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             for (int i = 0; i < RowCount; i++)
             {
                 for (int j = 0; j < ColumnCount; j++)
                 {
-                    At(i, j, f(At(i, j)));
+                    target.At(i, j, f(At(i, j)));
                 }
             }
         }
 
-        public virtual void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
+        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
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
+            {
+                var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
+                throw new ArgumentException(message, "target");
+            }
+
+            MapIndexedToUnchecked(target, f, forceMapZeros, skipClearing);
+        }
+
+        internal virtual void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
+            where TU : struct, IEquatable<TU>, IFormattable
         {
             for (int i = 0; i < RowCount; i++)
             {
                 for (int j = 0; j < ColumnCount; j++)
                 {
-                    At(i, j, f(i, j, At(i, j)));
+                    target.At(i, j, f(i, j, At(i, j)));
                 }
             }
         }

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -30,6 +30,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using MathNet.Numerics.Properties;
 
@@ -42,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         // [ruegg] public fields are OK here
 
         /// <summary>
-        /// The array containing the row indices of the existing rows. Element "i" of the array gives the index of the 
+        /// The array containing the row indices of the existing rows. Element "i" of the array gives the index of the
         /// element in the <see cref="Values"/> array that is first non-zero element in a row "i".
         /// The last value is equal to ValueCount, so that the number of non-zero entries in row "i" is always
         /// given by RowPointers[i+i] - RowPointers[i]. This array thus has length RowCount+1.
@@ -50,13 +51,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         public readonly int[] RowPointers;
 
         /// <summary>
-        /// An array containing the column indices of the non-zero values. Element "j" of the array 
+        /// An array containing the column indices of the non-zero values. Element "j" of the array
         /// is the number of the column in matrix that contains the j-th value in the <see cref="Values"/> 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 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;
@@ -214,7 +215,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
             valueCount -= 1;
 
-            // Check whether we need to shrink the arrays. This is reasonable to do if 
+            // Check whether we need to shrink the arrays. 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))
             {
@@ -374,7 +375,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         /// Returns a hash code for this instance.
         /// </summary>
         /// <returns>
-        /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table. 
+        /// A hash code for this instance, suitable for use in hashing algorithms and data structures like a hash table.
         /// </returns>
         public override int GetHashCode()
         {
@@ -1104,25 +1105,82 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
 
         // FUNCTIONAL COMBINATORS
 
-        public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
+        internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
         {
-            var newRowPointers = new int[RowCount+1];
-            var newColumnIndices = new List<int>();
-            var newValues = new List<T>();
+            var processZeros = forceMapZeros || !Zero.Equals(f(Zero));
+
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<TU>;
+            if (sparseTarget != null)
+            {
+                var newRowPointers = sparseTarget.RowPointers;
+                var newColumnIndices = new List<int>();
+                var newValues = new List<TU>();
+
+                if (processZeros)
+                {
+                    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);
+                            }
+                        }
+                    }
+                }
+                else
+                {
+                    for (int row = 0; row < RowCount; row++)
+                    {
+                        newRowPointers[row] = newValues.Count;
+                        var startIndex = RowPointers[row];
+                        var endIndex = RowPointers[row + 1];
+                        for (var j = startIndex; j < endIndex; j++)
+                        {
+                            var item = f(Values[j]);
+                            if (!Zero.Equals(item))
+                            {
+                                newValues.Add(item);
+                                newColumnIndices.Add(ColumnIndices[j]);
+                            }
+                        }
+                    }
+                }
+
+                sparseTarget.ColumnIndices = newColumnIndices.ToArray();
+                sparseTarget.Values = newValues.ToArray();
+                newRowPointers[RowCount] = newValues.Count;
+                return;
+            }
+
+            // FALL BACK
 
-            if (forceMapZeros || !Zero.Equals(f(Zero)))
+            if (!skipClearing && !processZeros)
+            {
+                target.Clear();
+            }
+
+            if (processZeros)
             {
-                int k = 0;
                 for (int row = 0; row < RowCount; row++)
                 {
-                    newRowPointers[row] = newValues.Count;
-                    for (int col = 0; col < ColumnCount; col++)
+                    var index = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    for (int j = 0; j < ColumnCount; j++)
                     {
-                        var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero);
-                        if (!Zero.Equals(item))
+                        if (j == ColumnIndices[index])
                         {
-                            newValues.Add(item);
-                            newColumnIndices.Add(col);
+                            target.At(row, j, f(Values[j]));
+                            index = Math.Min(index + 1, endIndex);
+                        }
+                        else
+                        {
+                            target.At(row, j, f(Zero));
                         }
                     }
                 }
@@ -1131,46 +1189,92 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (int row = 0; row < RowCount; row++)
                 {
-                    newRowPointers[row] = newValues.Count;
                     var startIndex = RowPointers[row];
                     var endIndex = RowPointers[row + 1];
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var item = f(Values[j]);
-                        if (!Zero.Equals(item))
+                        target.At(row, ColumnIndices[j], f(Values[j]));
+                    }
+                }
+            }
+        }
+
+        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));
+
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<TU>;
+            if (sparseTarget != null)
+            {
+                var newRowPointers = sparseTarget.RowPointers;
+                var newColumnIndices = new List<int>();
+                var newValues = new List<TU>();
+
+                if (processZeros)
+                {
+                    int k = 0;
+                    for (int row = 0; row < RowCount; row++)
+                    {
+                        newRowPointers[row] = newValues.Count;
+                        for (int col = 0; col < ColumnCount; col++)
                         {
-                            newValues.Add(item);
-                            newColumnIndices.Add(ColumnIndices[j]);
+                            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);
+                            }
                         }
                     }
                 }
+                else
+                {
+                    for (int row = 0; row < RowCount; row++)
+                    {
+                        newRowPointers[row] = newValues.Count;
+                        var startIndex = RowPointers[row];
+                        var endIndex = RowPointers[row + 1];
+                        for (var j = startIndex; j < endIndex; j++)
+                        {
+                            var item = f(row, ColumnIndices[j], Values[j]);
+                            if (!Zero.Equals(item))
+                            {
+                                newValues.Add(item);
+                                newColumnIndices.Add(ColumnIndices[j]);
+                            }
+                        }
+                    }
+                }
+
+                sparseTarget.ColumnIndices = newColumnIndices.ToArray();
+                sparseTarget.Values = newValues.ToArray();
+                newRowPointers[RowCount] = newValues.Count;
+                return;
             }
 
-            ColumnIndices = newColumnIndices.ToArray();
-            Values = newValues.ToArray();
-            newRowPointers[RowCount] = newValues.Count;
-            Array.Copy(newRowPointers, RowPointers, newRowPointers.Length);
-        }
+            // FALL BACK
 
-        public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
-        {
-            var newRowPointers = new int[RowCount+1];
-            var newColumnIndices = new List<int>();
-            var newValues = new List<T>();
+            if (!skipClearing && !processZeros)
+            {
+                target.Clear();
+            }
 
-            if (forceMapZeros || !Zero.Equals(f(0, 0, Zero)))
+            if (processZeros)
             {
-                int k = 0;
                 for (int row = 0; row < RowCount; row++)
                 {
-                    newRowPointers[row] = newValues.Count;
-                    for (int col = 0; col < ColumnCount; col++)
+                    var index = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    for (int j = 0; j < ColumnCount; j++)
                     {
-                        var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero);
-                        if (!Zero.Equals(item))
+                        if (j == ColumnIndices[index])
+                        {
+                            target.At(row, j, f(row, j, Values[j]));
+                            index = Math.Min(index + 1, endIndex);
+                        }
+                        else
                         {
-                            newValues.Add(item);
-                            newColumnIndices.Add(col);
+                            target.At(row, j, f(row, j, Zero));
                         }
                     }
                 }
@@ -1179,25 +1283,14 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             {
                 for (int row = 0; row < RowCount; row++)
                 {
-                    newRowPointers[row] = newValues.Count;
                     var startIndex = RowPointers[row];
                     var endIndex = RowPointers[row + 1];
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var item = f(row, ColumnIndices[j], Values[j]);
-                        if (!Zero.Equals(item))
-                        {
-                            newValues.Add(item);
-                            newColumnIndices.Add(ColumnIndices[j]);
-                        }
+                        target.At(row, ColumnIndices[j], f(row, ColumnIndices[j], Values[j]));
                     }
                 }
             }
-
-            ColumnIndices = newColumnIndices.ToArray();
-            Values = newValues.ToArray();
-            newRowPointers[RowCount] = newValues.Count;
-            Array.Copy(newRowPointers, RowPointers, newRowPointers.Length);
         }
     }
 }