From b50d4417a9e7567fe370249caf8450c8e2743b45 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 15 Jun 2014 23:35:00 +0200 Subject: [PATCH] LA: Matrix.FoldByRow/Column operates on arrays to allow arbitrary types (not just those supported by vectors) --- src/Numerics/LinearAlgebra/Complex/Matrix.cs | 24 ++++++------- .../LinearAlgebra/Complex32/Matrix.cs | 24 ++++++------- src/Numerics/LinearAlgebra/Double/Matrix.cs | 24 ++++++------- src/Numerics/LinearAlgebra/Matrix.cs | 35 +++++++++++++----- src/Numerics/LinearAlgebra/Single/Matrix.cs | 24 ++++++------- .../Storage/DenseColumnMajorMatrixStorage.cs | 12 +++---- .../Storage/DiagonalMatrixStorage.cs | 20 +++++------ .../LinearAlgebra/Storage/MatrixStorage.cs | 24 ++++++------- .../SparseCompressedRowMatrixStorage.cs | 36 ++++++++----------- .../MatrixStructureTheory.Functional.cs | 20 +++++------ 10 files changed, 128 insertions(+), 115 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index c62b9c01..9be6cdbc 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -117,25 +117,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(RowCount); + var ret = new double[RowCount]; if (norm == 2.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// @@ -149,25 +149,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(ColumnCount); + var ret = new double[ColumnCount]; if (norm == 2.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index 1f1d0042..764b14a5 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -111,25 +111,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(RowCount); + var ret = new double[RowCount]; if (norm == 2.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// @@ -143,25 +143,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(ColumnCount); + var ret = new double[ColumnCount]; if (norm == 2.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x.Magnitude, (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index aac56517..07086314 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -109,25 +109,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(RowCount); + var ret = new double[RowCount]; if (norm == 2.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// @@ -141,25 +141,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(ColumnCount); + var ret = new double[ColumnCount]; if (norm == 2.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// diff --git a/src/Numerics/LinearAlgebra/Matrix.cs b/src/Numerics/LinearAlgebra/Matrix.cs index 7ca21001..86e6ca4d 100644 --- a/src/Numerics/LinearAlgebra/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Matrix.cs @@ -33,6 +33,7 @@ using MathNet.Numerics.Properties; using System; using System.Collections.Generic; using System.Runtime; +using MathNet.Numerics.Threading; namespace MathNet.Numerics.LinearAlgebra { @@ -1614,11 +1615,20 @@ namespace MathNet.Numerics.LinearAlgebra /// For each row, applies a function f to each element of the row, threading an accumulator argument through the computation. /// Returns a vector with the resulting accumulator states for each row. /// - public Vector FoldByRow(Func f, TU state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + public TU[] FoldByRow(Func f, TU state, Zeros zeros = Zeros.AllowSkip) { - var result = Vector.Build.SameAs(this, RowCount); - Storage.FoldByRowUnchecked(result.Storage, f, (x, c) => x, DenseVectorStorage.OfInit(RowCount, i => state), zeros); + var result = new TU[RowCount]; + if (!EqualityComparer.Default.Equals(state, default(TU))) + { + CommonParallel.For(0, result.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = state; + } + }); + } + Storage.FoldByRowUnchecked(result, f, (x, c) => x, result, zeros); return result; } @@ -1626,11 +1636,20 @@ namespace MathNet.Numerics.LinearAlgebra /// For each column, applies a function f to each element of the column, threading an accumulator argument through the computation. /// Returns a vector with the resulting accumulator states for each column. /// - public Vector FoldByColumn(Func f, TU state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + public TU[] FoldByColumn(Func f, TU state, Zeros zeros = Zeros.AllowSkip) { - var result = Vector.Build.SameAs(this, ColumnCount); - Storage.FoldByColumnUnchecked(result.Storage, f, (x, c) => x, DenseVectorStorage.OfInit(ColumnCount, i => state), zeros); + var result = new TU[ColumnCount]; + if (!EqualityComparer.Default.Equals(state, default(TU))) + { + CommonParallel.For(0, result.Length, 4096, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = state; + } + }); + } + Storage.FoldByColumnUnchecked(result, f, (x, c) => x, result, zeros); return result; } } diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 20a6dd69..098111e1 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -109,25 +109,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(RowCount); + var ret = new double[RowCount]; if (norm == 2.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByRowUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByRowUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// @@ -141,25 +141,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); } - var ret = Vector.Build.Dense(ColumnCount); + var ret = new double[ColumnCount]; if (norm == 2.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret, Zeros.AllowSkip); } else if (norm == 1.0) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Abs(x), (x, c) => x, ret, Zeros.AllowSkip); } else if (double.IsPositiveInfinity(norm)) { - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret, Zeros.AllowSkip); } else { double invnorm = 1.0/norm; - Storage.FoldByColumnUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip); + Storage.FoldByColumnUnchecked(ret, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret, Zeros.AllowSkip); } - return ret; + return Vector.Build.Dense(ret); } /// diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs index ceef2500..29a3e660 100644 --- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs @@ -727,30 +727,30 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS: FOLD - internal override void FoldByRowUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByRowUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { for (int i = 0; i < RowCount; i++) { - TU s = state.At(i); + TU s = state[i]; for (int j = 0; j < ColumnCount; j++) { s = f(s, Data[j*RowCount + i]); } - target.At(i, finalize(s, ColumnCount)); + target[i] = finalize(s, ColumnCount); } } - internal override void FoldByColumnUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByColumnUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { for (int j = 0; j < ColumnCount; j++) { int offset = j*RowCount; - TU s = state.At(j); + TU s = state[j]; for (int i = 0; i < RowCount; i++) { s = f(s, Data[offset + i]); } - target.At(j, finalize(s, RowCount)); + target[j] = finalize(s, RowCount); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs index 810b1d7d..3bd8c9ec 100644 --- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs @@ -901,58 +901,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS: FOLD - internal override void FoldByRowUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByRowUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { if (zeros == Zeros.AllowSkip) { for (int k = 0; k < Data.Length; k++) { - target.At(k, finalize(f(state.At(k), Data[k]), 1)); + target[k] = finalize(f(state[k], Data[k]), 1); } for (int k = Data.Length; k < RowCount; k++) { - target.At(k, finalize(state.At(k), 0)); + target[k] = finalize(state[k], 0); } } else { for (int i = 0; i < RowCount; i++) { - TU s = state.At(i); + TU s = state[i]; for (int j = 0; j < ColumnCount; j++) { s = f(s, i == j ? Data[i] : Zero); } - target.At(i, finalize(s, ColumnCount)); + target[i] = finalize(s, ColumnCount); } } } - internal override void FoldByColumnUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByColumnUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { if (zeros == Zeros.AllowSkip) { for (int k = 0; k < Data.Length; k++) { - target.At(k, finalize(f(state.At(k), Data[k]), 1)); + target[k] = finalize(f(state[k], Data[k]), 1); } for (int k = Data.Length; k < ColumnCount; k++) { - target.At(k, finalize(state.At(k), 0)); + target[k] = finalize(state[k], 0); } } else { for (int j = 0; j < ColumnCount; j++) { - TU s = state.At(j); + TU s = state[j]; for (int i = 0; i < RowCount; i++) { s = f(s, i == j ? Data[i] : Zero); } - target.At(j, finalize(s, RowCount)); + target[j] = finalize(s, RowCount); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs index ab55d2d9..83d6f0ab 100644 --- a/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs @@ -670,8 +670,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS: FOLD - public void FoldByRow(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + /// The state array will not be modified, unless it is the same instance as the target array (which is allowed). + public void FoldByRow(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { if (target == null) { @@ -694,22 +694,22 @@ namespace MathNet.Numerics.LinearAlgebra.Storage FoldByRowUnchecked(target, f, finalize, state, zeros); } - internal virtual void FoldByRowUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + /// The state array will not be modified, unless it is the same instance as the target array (which is allowed). + internal virtual void FoldByRowUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { for (int i = 0; i < RowCount; i++) { - TU s = state.At(i); + TU s = state[i]; for (int j = 0; j < ColumnCount; j++) { s = f(s, At(i, j)); } - target.At(i, finalize(s, ColumnCount)); + target[i] = finalize(s, ColumnCount); } } - public void FoldByColumn(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + /// The state array will not be modified, unless it is the same instance as the target array (which is allowed). + public void FoldByColumn(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { if (target == null) { @@ -732,17 +732,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage FoldByColumnUnchecked(target, f, finalize, state, zeros); } - internal virtual void FoldByColumnUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) - where TU : struct, IEquatable, IFormattable + /// The state array will not be modified, unless it is the same instance as the target array (which is allowed). + internal virtual void FoldByColumnUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { for (int j = 0; j < ColumnCount; j++) { - TU s = state.At(j); + TU s = state[j]; for (int i = 0; i < RowCount; i++) { s = f(s, At(i, j)); } - target.At(j, finalize(s, RowCount)); + target[j] = finalize(s, RowCount); } } } diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs index 084664fc..25417f61 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs @@ -1757,7 +1757,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage // FUNCTIONAL COMBINATORS: FOLD - internal override void FoldByRowUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByRowUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { if (zeros == Zeros.AllowSkip) { @@ -1765,12 +1765,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { var startIndex = RowPointers[row]; var endIndex = RowPointers[row + 1]; - TU s = state.At(row); + TU s = state[row]; for (var j = startIndex; j < endIndex; j++) { s = f(s, Values[j]); } - target.At(row, finalize(s, endIndex - startIndex)); + target[row] = finalize(s, endIndex - startIndex); } } else @@ -1779,7 +1779,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { var index = RowPointers[row]; var endIndex = RowPointers[row + 1]; - TU s = state.At(row); + TU s = state[row]; for (int j = 0; j < ColumnCount; j++) { if (index < endIndex && j == ColumnIndices[index]) @@ -1792,18 +1792,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage s = f(s, Zero); } } - target.At(row, finalize(s, ColumnCount)); + target[row] = finalize(s, ColumnCount); } } } - internal override void FoldByColumnUnchecked(VectorStorage target, Func f, Func finalize, VectorStorage state, Zeros zeros = Zeros.AllowSkip) + internal override void FoldByColumnUnchecked(TU[] target, Func f, Func finalize, TU[] state, Zeros zeros = Zeros.AllowSkip) { - var denseResult = target as DenseVectorStorage ?? new DenseVectorStorage(ColumnCount); - - state.CopyTo(denseResult); - TU[] result = denseResult.Data; - + if (!ReferenceEquals(state, target)) + { + Array.Copy(state, target, state.Length); + } if (zeros == Zeros.AllowSkip) { int[] count = new int[ColumnCount]; @@ -1814,13 +1813,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage for (var j = startIndex; j < endIndex; j++) { var column = ColumnIndices[j]; - result[column] = f(result[column], Values[j]); + target[column] = f(target[column], Values[j]); count[column]++; } } for (int j = 0; j < ColumnCount; j++) { - result[j] = finalize(result[j], count[j]); + target[j] = finalize(target[j], count[j]); } } else @@ -1833,25 +1832,20 @@ namespace MathNet.Numerics.LinearAlgebra.Storage { if (index < endIndex && j == ColumnIndices[index]) { - result[j] = f(result[j], Values[index]); + target[j] = f(target[j], Values[index]); index = Math.Min(index + 1, endIndex); } else { - result[j] = f(result[j], Zero); + target[j] = f(target[j], Zero); } } } for (int j = 0; j < ColumnCount; j++) { - result[j] = finalize(result[j], RowCount); + target[j] = finalize(target[j], RowCount); } } - - if (!ReferenceEquals(denseResult, target)) - { - denseResult.CopyTo(target); - } } } } diff --git a/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs index b1fb70eb..44ce7271 100644 --- a/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs +++ b/src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs @@ -261,17 +261,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests public void CanFoldRows(Matrix matrix) { // not forced - var rowSum = matrix.FoldByRow((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.AllowSkip); - for (int i = 0; i < rowSum.Count; i++) + T[] rowSum = matrix.FoldByRow((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.AllowSkip); + for (int i = 0; i < rowSum.Length; i++) { - Assert.That(rowSum.At(i), Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "not forced"); + Assert.That(rowSum[i], Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "not forced"); } // forced rowSum = matrix.FoldByRow((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.Include); - for (int i = 0; i < rowSum.Count; i++) + for (int i = 0; i < rowSum.Length; i++) { - Assert.That(rowSum.At(i), Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "forced"); + Assert.That(rowSum[i], Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "forced"); } Assert.That(matrix.FoldByRow((s, x) => s + 1.0, 0.0, Zeros.Include), @@ -282,17 +282,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests public void CanFoldColumns(Matrix matrix) { // not forced - var colSum = matrix.FoldByColumn((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.AllowSkip); - for (int i = 0; i < colSum.Count; i++) + T[] colSum = matrix.FoldByColumn((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.AllowSkip); + for (int i = 0; i < colSum.Length; i++) { - Assert.That(colSum.At(i), Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "not forced"); + Assert.That(colSum[i], Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "not forced"); } // forced colSum = matrix.FoldByColumn((s, x) => Operator.Add(s, x), Operator.Zero, Zeros.Include); - for (int i = 0; i < colSum.Count; i++) + for (int i = 0; i < colSum.Length; i++) { - Assert.That(colSum.At(i), Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "forced"); + Assert.That(colSum[i], Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator.Add(a, b))), "forced"); } Assert.That(matrix.FoldByColumn((s, x) => s + 1.0, 0.0, Zeros.Include),