From ebc5471d5c873a004f7cd45acd5dea2e9d7706d4 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 25 Apr 2014 01:37:11 +0200 Subject: [PATCH] LA: special case for matrix inplace map --- src/Numerics/LinearAlgebra/Matrix.cs | 4 +- .../Storage/DenseColumnMajorMatrixStorage.cs | 42 ++++++- .../Storage/DiagonalMatrixStorage.cs | 32 +++++ .../LinearAlgebra/Storage/MatrixStorage.cs | 22 ++++ .../SparseCompressedRowMatrixStorage.cs | 114 +++++++++++++++++- 5 files changed, 202 insertions(+), 12 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Matrix.cs b/src/Numerics/LinearAlgebra/Matrix.cs index 245973f0..991b4664 100644 --- a/src/Numerics/LinearAlgebra/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Matrix.cs @@ -1473,7 +1473,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapInplace(Func f, bool forceMapZeros = false) { - Storage.MapToUnchecked(Storage, f, forceMapZeros, skipClearing: true); + Storage.MapInplace(f, forceMapZeros); } /// @@ -1484,7 +1484,7 @@ namespace MathNet.Numerics.LinearAlgebra /// public void MapIndexedInplace(Func f, bool forceMapZeros = false) { - Storage.MapIndexedToUnchecked(Storage, f, forceMapZeros, skipClearing: true); + Storage.MapIndexedInplace(f, forceMapZeros); } /// diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs index 35947b20..ec8cce0c 100644 --- 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 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(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) { var denseTarget = target as DenseColumnMajorMatrixStorage; 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 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++; + } + } + }); + } } } diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs index b863ec8e..d42a0f32 100644 --- 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 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(MatrixStorage target, Func 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 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]); + } + }); + } } } diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs index 754e7514..98a92bf8 100644 --- 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 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(MatrixStorage target, Func f, bool forceMapZeros = false, bool skipClearing = false) where TU : struct, IEquatable, IFormattable { @@ -540,5 +551,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } } + + public virtual void MapIndexedInplace(Func 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))); + } + } + } } } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs index 3d8c7e5a..ec85d0e1 100644 --- 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(); - var newValues = new List(); + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(Values.Length); if (processZeros) { @@ -1199,6 +1199,59 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } + public override void MapInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros || !Zero.Equals(f(Zero))) + { + var newRowPointers = RowPointers; + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(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(MatrixStorage target, Func 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(); - var newValues = new List(); + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(Values.Length); if (processZeros) { @@ -1292,5 +1345,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage } } } + + public override void MapIndexedInplace(Func f, bool forceMapZeros = false) + { + if (forceMapZeros || !Zero.Equals(f(0, 0, Zero))) + { + var newRowPointers = RowPointers; + var newColumnIndices = new List(ColumnIndices.Length); + var newValues = new List(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; + } + } } }