Browse Source

LA: Matrix.FoldByRow/Column operates on arrays to allow arbitrary types (not just those supported by vectors)

provider
Christoph Ruegg 12 years ago
parent
commit
b50d4417a9
  1. 24
      src/Numerics/LinearAlgebra/Complex/Matrix.cs
  2. 24
      src/Numerics/LinearAlgebra/Complex32/Matrix.cs
  3. 24
      src/Numerics/LinearAlgebra/Double/Matrix.cs
  4. 35
      src/Numerics/LinearAlgebra/Matrix.cs
  5. 24
      src/Numerics/LinearAlgebra/Single/Matrix.cs
  6. 12
      src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
  7. 20
      src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
  8. 24
      src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
  9. 36
      src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
  10. 20
      src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs

24
src/Numerics/LinearAlgebra/Complex/Matrix.cs

@ -117,25 +117,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(RowCount); var ret = new double[RowCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>
@ -149,25 +149,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(ColumnCount); var ret = new double[ColumnCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>

24
src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@ -111,25 +111,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(RowCount); var ret = new double[RowCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>
@ -143,25 +143,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(ColumnCount); var ret = new double[ColumnCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>

24
src/Numerics/LinearAlgebra/Double/Matrix.cs

@ -109,25 +109,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(RowCount); var ret = new double[RowCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>
@ -141,25 +141,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(ColumnCount); var ret = new double[ColumnCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>

35
src/Numerics/LinearAlgebra/Matrix.cs

@ -33,6 +33,7 @@ using MathNet.Numerics.Properties;
using System; using System;
using System.Collections.Generic; using System.Collections.Generic;
using System.Runtime; using System.Runtime;
using MathNet.Numerics.Threading;
namespace MathNet.Numerics.LinearAlgebra 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. /// 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. /// Returns a vector with the resulting accumulator states for each row.
/// </summary> /// </summary>
public Vector<TU> FoldByRow<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip) public TU[] FoldByRow<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip)
where TU : struct, IEquatable<TU>, IFormattable
{ {
var result = Vector<TU>.Build.SameAs(this, RowCount); var result = new TU[RowCount];
Storage.FoldByRowUnchecked(result.Storage, f, (x, c) => x, DenseVectorStorage<TU>.OfInit(RowCount, i => state), zeros); if (!EqualityComparer<TU>.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; 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. /// 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. /// Returns a vector with the resulting accumulator states for each column.
/// </summary> /// </summary>
public Vector<TU> FoldByColumn<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip) public TU[] FoldByColumn<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip)
where TU : struct, IEquatable<TU>, IFormattable
{ {
var result = Vector<TU>.Build.SameAs(this, ColumnCount); var result = new TU[ColumnCount];
Storage.FoldByColumnUnchecked(result.Storage, f, (x, c) => x, DenseVectorStorage<TU>.OfInit(ColumnCount, i => state), zeros); if (!EqualityComparer<TU>.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; return result;
} }
} }

24
src/Numerics/LinearAlgebra/Single/Matrix.cs

@ -109,25 +109,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(RowCount); var ret = new double[RowCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>
@ -141,25 +141,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive); throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
} }
var ret = Vector<double>.Build.Dense(ColumnCount); var ret = new double[ColumnCount];
if (norm == 2.0) 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) 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)) 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 else
{ {
double invnorm = 1.0/norm; 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<double>.Build.Dense(ret);
} }
/// <summary> /// <summary>

12
src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@ -727,30 +727,30 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS: FOLD // FUNCTIONAL COMBINATORS: FOLD
internal override void FoldByRowUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByRowUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
TU s = state.At(i); TU s = state[i];
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
s = f(s, Data[j*RowCount + i]); s = f(s, Data[j*RowCount + i]);
} }
target.At(i, finalize(s, ColumnCount)); target[i] = finalize(s, ColumnCount);
} }
} }
internal override void FoldByColumnUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByColumnUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
int offset = j*RowCount; int offset = j*RowCount;
TU s = state.At(j); TU s = state[j];
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
s = f(s, Data[offset + i]); s = f(s, Data[offset + i]);
} }
target.At(j, finalize(s, RowCount)); target[j] = finalize(s, RowCount);
} }
} }
} }

20
src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@ -901,58 +901,58 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS: FOLD // FUNCTIONAL COMBINATORS: FOLD
internal override void FoldByRowUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByRowUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
if (zeros == Zeros.AllowSkip) if (zeros == Zeros.AllowSkip)
{ {
for (int k = 0; k < Data.Length; k++) 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++) for (int k = Data.Length; k < RowCount; k++)
{ {
target.At(k, finalize(state.At(k), 0)); target[k] = finalize(state[k], 0);
} }
} }
else else
{ {
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
TU s = state.At(i); TU s = state[i];
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
s = f(s, i == j ? Data[i] : Zero); s = f(s, i == j ? Data[i] : Zero);
} }
target.At(i, finalize(s, ColumnCount)); target[i] = finalize(s, ColumnCount);
} }
} }
} }
internal override void FoldByColumnUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByColumnUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
if (zeros == Zeros.AllowSkip) if (zeros == Zeros.AllowSkip)
{ {
for (int k = 0; k < Data.Length; k++) 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++) for (int k = Data.Length; k < ColumnCount; k++)
{ {
target.At(k, finalize(state.At(k), 0)); target[k] = finalize(state[k], 0);
} }
} }
else else
{ {
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
TU s = state.At(j); TU s = state[j];
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
s = f(s, i == j ? Data[i] : Zero); s = f(s, i == j ? Data[i] : Zero);
} }
target.At(j, finalize(s, RowCount)); target[j] = finalize(s, RowCount);
} }
} }
} }

24
src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@ -670,8 +670,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS: FOLD // FUNCTIONAL COMBINATORS: FOLD
public void FoldByRow<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) /// <remarks>The state array will not be modified, unless it is the same instance as the target array (which is allowed).</remarks>
where TU : struct, IEquatable<TU>, IFormattable public void FoldByRow<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
if (target == null) if (target == null)
{ {
@ -694,22 +694,22 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
FoldByRowUnchecked(target, f, finalize, state, zeros); FoldByRowUnchecked(target, f, finalize, state, zeros);
} }
internal virtual void FoldByRowUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) /// <remarks>The state array will not be modified, unless it is the same instance as the target array (which is allowed).</remarks>
where TU : struct, IEquatable<TU>, IFormattable internal virtual void FoldByRowUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
TU s = state.At(i); TU s = state[i];
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
s = f(s, At(i, j)); s = f(s, At(i, j));
} }
target.At(i, finalize(s, ColumnCount)); target[i] = finalize(s, ColumnCount);
} }
} }
public void FoldByColumn<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) /// <remarks>The state array will not be modified, unless it is the same instance as the target array (which is allowed).</remarks>
where TU : struct, IEquatable<TU>, IFormattable public void FoldByColumn<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
if (target == null) if (target == null)
{ {
@ -732,17 +732,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
FoldByColumnUnchecked(target, f, finalize, state, zeros); FoldByColumnUnchecked(target, f, finalize, state, zeros);
} }
internal virtual void FoldByColumnUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) /// <remarks>The state array will not be modified, unless it is the same instance as the target array (which is allowed).</remarks>
where TU : struct, IEquatable<TU>, IFormattable internal virtual void FoldByColumnUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
TU s = state.At(j); TU s = state[j];
for (int i = 0; i < RowCount; i++) for (int i = 0; i < RowCount; i++)
{ {
s = f(s, At(i, j)); s = f(s, At(i, j));
} }
target.At(j, finalize(s, RowCount)); target[j] = finalize(s, RowCount);
} }
} }
} }

36
src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@ -1757,7 +1757,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS: FOLD // FUNCTIONAL COMBINATORS: FOLD
internal override void FoldByRowUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByRowUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
if (zeros == Zeros.AllowSkip) if (zeros == Zeros.AllowSkip)
{ {
@ -1765,12 +1765,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{ {
var startIndex = RowPointers[row]; var startIndex = RowPointers[row];
var endIndex = RowPointers[row + 1]; var endIndex = RowPointers[row + 1];
TU s = state.At(row); TU s = state[row];
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
s = f(s, Values[j]); s = f(s, Values[j]);
} }
target.At(row, finalize(s, endIndex - startIndex)); target[row] = finalize(s, endIndex - startIndex);
} }
} }
else else
@ -1779,7 +1779,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{ {
var index = RowPointers[row]; var index = RowPointers[row];
var endIndex = RowPointers[row + 1]; var endIndex = RowPointers[row + 1];
TU s = state.At(row); TU s = state[row];
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
if (index < endIndex && j == ColumnIndices[index]) if (index < endIndex && j == ColumnIndices[index])
@ -1792,18 +1792,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
s = f(s, Zero); s = f(s, Zero);
} }
} }
target.At(row, finalize(s, ColumnCount)); target[row] = finalize(s, ColumnCount);
} }
} }
} }
internal override void FoldByColumnUnchecked<TU>(VectorStorage<TU> target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, VectorStorage<TU> state, Zeros zeros = Zeros.AllowSkip) internal override void FoldByColumnUnchecked<TU>(TU[] target, Func<TU, T, TU> f, Func<TU, int, TU> finalize, TU[] state, Zeros zeros = Zeros.AllowSkip)
{ {
var denseResult = target as DenseVectorStorage<TU> ?? new DenseVectorStorage<TU>(ColumnCount); if (!ReferenceEquals(state, target))
{
state.CopyTo(denseResult); Array.Copy(state, target, state.Length);
TU[] result = denseResult.Data; }
if (zeros == Zeros.AllowSkip) if (zeros == Zeros.AllowSkip)
{ {
int[] count = new int[ColumnCount]; int[] count = new int[ColumnCount];
@ -1814,13 +1813,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
var column = ColumnIndices[j]; var column = ColumnIndices[j];
result[column] = f(result[column], Values[j]); target[column] = f(target[column], Values[j]);
count[column]++; count[column]++;
} }
} }
for (int j = 0; j < ColumnCount; j++) for (int j = 0; j < ColumnCount; j++)
{ {
result[j] = finalize(result[j], count[j]); target[j] = finalize(target[j], count[j]);
} }
} }
else else
@ -1833,25 +1832,20 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{ {
if (index < endIndex && j == ColumnIndices[index]) 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); index = Math.Min(index + 1, endIndex);
} }
else else
{ {
result[j] = f(result[j], Zero); target[j] = f(target[j], Zero);
} }
} }
} }
for (int j = 0; j < ColumnCount; j++) 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);
}
} }
} }
} }

20
src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Functional.cs

@ -261,17 +261,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests
public void CanFoldRows(Matrix<T> matrix) public void CanFoldRows(Matrix<T> matrix)
{ {
// not forced // not forced
var rowSum = matrix.FoldByRow((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.AllowSkip); T[] rowSum = matrix.FoldByRow((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.AllowSkip);
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<T>.Add(a, b))), "not forced"); Assert.That(rowSum[i], Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator<T>.Add(a, b))), "not forced");
} }
// forced // forced
rowSum = matrix.FoldByRow((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.Include); rowSum = matrix.FoldByRow((s, x) => Operator<T>.Add(s, x), Operator<T>.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<T>.Add(a, b))), "forced"); Assert.That(rowSum[i], Is.EqualTo(matrix.Row(i).Enumerate().Aggregate((a, b) => Operator<T>.Add(a, b))), "forced");
} }
Assert.That(matrix.FoldByRow((s, x) => s + 1.0, 0.0, Zeros.Include), 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<T> matrix) public void CanFoldColumns(Matrix<T> matrix)
{ {
// not forced // not forced
var colSum = matrix.FoldByColumn((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.AllowSkip); T[] colSum = matrix.FoldByColumn((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.AllowSkip);
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<T>.Add(a, b))), "not forced"); Assert.That(colSum[i], Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator<T>.Add(a, b))), "not forced");
} }
// forced // forced
colSum = matrix.FoldByColumn((s, x) => Operator<T>.Add(s, x), Operator<T>.Zero, Zeros.Include); colSum = matrix.FoldByColumn((s, x) => Operator<T>.Add(s, x), Operator<T>.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<T>.Add(a, b))), "forced"); Assert.That(colSum[i], Is.EqualTo(matrix.Column(i).Enumerate().Aggregate((a, b) => Operator<T>.Add(a, b))), "forced");
} }
Assert.That(matrix.FoldByColumn((s, x) => s + 1.0, 0.0, Zeros.Include), Assert.That(matrix.FoldByColumn((s, x) => s + 1.0, 0.0, Zeros.Include),

Loading…
Cancel
Save