From 1ca9cfc1626005230dc6894918b426e7c59f2cf4 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 15 Aug 2014 14:09:30 +0200 Subject: [PATCH] Perf: matrix product bench with experimental providers --- .../ManagedLinearAlgebraProvider.Complex.cs | 2 +- .../ManagedLinearAlgebraProvider.Complex32.cs | 2 +- .../ManagedLinearAlgebraProvider.Double.cs | 2 +- .../ManagedLinearAlgebraProvider.Single.cs | 2 +- .../LinearAlgebra/DenseMatrixProduct.cs | 557 ++++++++++++++++++ .../LinearAlgebra/DenseVectorAdd.cs | 132 +++-- src/Performance/Performance.csproj | 4 + src/Performance/Program.cs | 47 +- .../StatisticsTests/StatisticsTests.cs | 10 + 9 files changed, 693 insertions(+), 65 deletions(-) create mode 100644 src/Performance/LinearAlgebra/DenseMatrixProduct.cs diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index e4a0812d..4dfb5cf2 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -632,7 +632,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } else if (!beta.IsOne()) { - Control.LinearAlgebraProvider.ScaleArray(beta, c, c); + ScaleArray(beta, c, c); } if (alpha.IsZero()) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 02c950c4..87bc50ff 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -629,7 +629,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } else if (!beta.IsOne()) { - Control.LinearAlgebraProvider.ScaleArray(beta, c, c); + ScaleArray(beta, c, c); } if (alpha.IsZero()) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index d1e3318a..77f49d6f 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -624,7 +624,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } else if (beta != 1.0) { - Control.LinearAlgebraProvider.ScaleArray(beta, c, c); + ScaleArray(beta, c, c); } if (alpha == 0.0) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 98dafcca..cc957a7d 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -624,7 +624,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } else if (beta != 1.0f) { - Control.LinearAlgebraProvider.ScaleArray(beta, c, c); + ScaleArray(beta, c, c); } if (alpha == 0.0f) diff --git a/src/Performance/LinearAlgebra/DenseMatrixProduct.cs b/src/Performance/LinearAlgebra/DenseMatrixProduct.cs new file mode 100644 index 00000000..e62c782d --- /dev/null +++ b/src/Performance/LinearAlgebra/DenseMatrixProduct.cs @@ -0,0 +1,557 @@ +using System; +using Binarysharp.Benchmark; +using MathNet.Numerics; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Providers.LinearAlgebra; +using MathNet.Numerics.Providers.LinearAlgebra.Mkl; +using MathNet.Numerics.Threading; + +namespace Performance.LinearAlgebra +{ + public class DenseMatrixProduct + { + readonly int _rounds; + readonly Matrix _a; + readonly Matrix _b; + + readonly ILinearAlgebraProvider _managed = new ManagedLinearAlgebraProvider(); + readonly ILinearAlgebraProvider _mkl = new MklLinearAlgebraProvider(); + readonly ILinearAlgebraProvider _safeProvider = new SafeProvider(); + readonly ILinearAlgebraProvider _unsafeProvider = new UnsafeProvider(); + readonly ILinearAlgebraProvider _experimentalProvider = new ExperimentalProvider(); + + public DenseMatrixProduct(int size, int rounds) + { + _rounds = rounds; + + _b = Matrix.Build.Random(size, size); + _a = Matrix.Build.Random(size, size); + + _managed.InitializeVerify(); + _safeProvider.InitializeVerify(); + _unsafeProvider.InitializeVerify(); + _experimentalProvider.InitializeVerify(); + +#if NATIVEMKL + _mkl.InitializeVerify(); +#endif + } + + public static void Verify(int size) + { + var x = new DenseMatrixProduct(size, 1); + var managedResult = x.ManagedProvider(); + var mklResult = x.MklProvider(); + var safeResult = x.SafeProvider(); + var unsafeResult = x.UnsafeProvider(); + var experimentalResult = x.ExperimentalProvider(); + + Console.WriteLine(managedResult.ToString()); + //Console.WriteLine(mklResult.ToString()); + //Console.WriteLine(safeResult.ToString()); + //Console.WriteLine(unsafeResult.ToString()); + //Console.WriteLine(experimentalResult.ToString()); + + if (!managedResult.AlmostEqual(mklResult, 1e-12)) + { + throw new Exception("MklProvider"); + } + if (!managedResult.AlmostEqual(safeResult, 1e-12)) + { + throw new Exception("SafeProvider"); + } + if (!managedResult.AlmostEqual(unsafeResult, 1e-12)) + { + throw new Exception("UnsafeProvider"); + } + if (!managedResult.AlmostEqual(experimentalResult, 1e-12)) + { + throw new Exception("ExperimentalProvider"); + } + } + + [BenchSharkTask("ManagedProvider")] + public Matrix ManagedProvider() + { + Control.LinearAlgebraProvider = _managed; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a*z; + } + return z; + } + + [BenchSharkTask("MklProvider")] + public Matrix MklProvider() + { + Control.LinearAlgebraProvider = _mkl; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a*z; + } + return z; + } + + [BenchSharkTask("SafeProvider")] + public Matrix SafeProvider() + { + Control.LinearAlgebraProvider = _safeProvider; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a*z; + } + return z; + } + + [BenchSharkTask("UnsafeProvider")] + public Matrix UnsafeProvider() + { + Control.LinearAlgebraProvider = _unsafeProvider; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a*z; + } + return z; + } + + [BenchSharkTask("ExperimentalProvider")] + public Matrix ExperimentalProvider() + { + Control.LinearAlgebraProvider = _experimentalProvider; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a*z; + } + return z; + } + } + + public class SafeProvider : ManagedLinearAlgebraProvider + { + public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + { + if (rowsX + columnsY <= Control.ParallelizeOrder) + { + for (int i = 0; i < rowsX; ++i) + { + for (int j = 0; j < columnsY; ++j) + { + var jrowsY = j*rowsY; + double sum = 0.0; + for (int k = 0; k < columnsX; ++k) + { + sum += x[k*rowsX + i]*y[jrowsY + k]; + } + result[j*rowsX + i] = sum; + } + } + + return; + } + + double[] xdata; + if (ReferenceEquals(x, result)) + { + xdata = (double[])x.Clone(); + } + else + { + xdata = x; + } + + double[] ydata; + if (ReferenceEquals(y, result)) + { + ydata = (double[])y.Clone(); + } + else + { + ydata = y; + } + + Array.Clear(result, 0, result.Length); + + CacheObliviousMatrixMultiply(xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, 0); + } + + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + { + if (transposeA == Transpose.DontTranspose && transposeB == Transpose.DontTranspose && alpha == 1.0 && beta == 0.0) + { + MatrixMultiply(a, rowsA, columnsA, b, rowsB, columnsB, c); + return; + } + + base.MatrixMultiplyWithUpdate(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c); + } + + static void CacheObliviousMatrixMultiply(double[] matrixA, int shiftArow, int shiftAcol, double[] matrixB, int shiftBrow, int shiftBcol, double[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, int level) + { + if (m + n <= Control.ParallelizeOrder) + { + for (var m1 = 0; m1 < m; m1++) + { + var matArowPos = m1 + shiftArow; + var matCrowPos = m1 + shiftCrow; + for (var n1 = 0; n1 < n; ++n1) + { + var boffset = ((n1 + shiftBcol)*constK) + shiftBrow; + double sum = 0; + for (var k1 = 0; k1 < k; ++k1) + { + sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*matrixB[boffset + k1]; + } + + result[((n1 + shiftCcol)*constM) + matCrowPos] += sum; + } + } + + return; + } + + // divide and conquer + int m2 = m/2, n2 = n/2, k2 = k/2; + + level++; + if (level <= 2) + { + CommonParallel.Invoke( + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, level)); + + CommonParallel.Invoke( + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, level)); + } + else + { + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, level); + } + } + } + + public unsafe class UnsafeProvider : ManagedLinearAlgebraProvider + { + public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + { + if (rowsX + columnsY <= Control.ParallelizeOrder) + { + fixed (double* resultPtr = &result[0]) + fixed (double* xPtr = &x[0]) + fixed (double* yPtr = &y[0]) + { + double* a = xPtr; + double* c = resultPtr; + for (int i = 0; i < rowsX; ++i) + { + double* b = yPtr; + double* cj = c; + for (int j = 0; j < columnsY; ++j) + { + double sum = 0.0; + for (int k = 0; k < columnsX; ++k) + { + sum += a[k*rowsX]*b[k]; + } + *cj = sum; + cj += rowsX; + b += rowsY; + } + a++; + c++; + } + } + + return; + } + + double[] xdata; + if (ReferenceEquals(x, result)) + { + xdata = (double[])x.Clone(); + } + else + { + xdata = x; + } + + double[] ydata; + if (ReferenceEquals(y, result)) + { + ydata = (double[])y.Clone(); + } + else + { + ydata = y; + } + + Array.Clear(result, 0, result.Length); + + CacheObliviousMatrixMultiply(xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, 0); + } + + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + { + if (transposeA == Transpose.DontTranspose && transposeB == Transpose.DontTranspose && alpha == 1.0 && beta == 0.0) + { + MatrixMultiply(a, rowsA, columnsA, b, rowsB, columnsB, c); + return; + } + + base.MatrixMultiplyWithUpdate(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c); + } + + static void CacheObliviousMatrixMultiply(double[] matrixA, int shiftArow, int shiftAcol, double[] matrixB, int shiftBrow, int shiftBcol, double[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, int level) + { + if (m + n <= Control.ParallelizeOrder) + { + fixed (double* resultPtr = &result[0]) + fixed (double* aPtr = &matrixA[0]) + fixed (double* bPtr = &matrixB[0]) + { + double* a = aPtr + shiftArow; + double* c = resultPtr + shiftCrow; + for (var m1 = 0; m1 < m; m1++) + { + for (var n1 = 0; n1 < n; ++n1) + { + double* b = bPtr + (n1 + shiftBcol)*constK + shiftBrow; + double sum = 0; + for (var k1 = 0; k1 < k; ++k1) + { + sum += a[((k1 + shiftAcol)*constM)]*b[k1]; + } + + c[((n1 + shiftCcol)*constM)] += sum; + } + a++; + c++; + } + } + + return; + } + + // divide and conquer + int m2 = m/2, n2 = n/2, k2 = k/2; + + level++; + if (level <= 2) + { + CommonParallel.Invoke( + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, level)); + + CommonParallel.Invoke( + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, level), + () => CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, level)); + } + else + { + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, level); + + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, level); + CacheObliviousMatrixMultiply(matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, level); + } + } + } + + public class ExperimentalProvider : ManagedLinearAlgebraProvider + { + public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); + } + + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + if (transposeA != Transpose.DontTranspose) + { + Swap(ref rowsA, ref columnsA); + } + + if (transposeB != Transpose.DontTranspose) + { + Swap(ref rowsB, ref columnsB); + } + + if (columnsA != rowsB) + { + throw new ArgumentOutOfRangeException(String.Format("columnsA ({0}) != rowsB ({1})", columnsA, rowsB)); + } + + if (rowsA*columnsA != a.Length) + { + throw new ArgumentOutOfRangeException(String.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsA, columnsA, a.Length)); + } + + if (rowsB*columnsB != b.Length) + { + throw new ArgumentOutOfRangeException(String.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsB, columnsB, b.Length)); + } + + if (rowsA*columnsB != c.Length) + { + throw new ArgumentOutOfRangeException(String.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsA, columnsB, c.Length)); + } + + // handle the degenerate cases + if (beta == 0.0) + { + Array.Clear(c, 0, c.Length); + } + else if (beta != 1.0) + { + ScaleArray(beta, c, c); + } + + if (alpha == 0.0) + { + return; + } + + // Extract column arrays + var columnDataB = new double[columnsB][]; + for (int i = 0; i < columnDataB.Length; i++) + { + columnDataB[i] = GetColumn(transposeB, i, rowsB, columnsB, b); + } + + var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; + if (shouldNotParallelize) + { + for (int i = 0; i < rowsA; i++) + { + var row = GetRow(transposeA, i, rowsA, columnsA, a); + for (int j = 0; j < columnsB; j++) + { + var col = columnDataB[j]; + double sum = 0; + for (int ii = 0; ii < row.Length; ii++) + { + sum += row[ii]*col[ii]; + } + + c[j*rowsA + i] += alpha*sum; + } + } + } + else + { + CommonParallel.For(0, rowsA, 1, (u, v) => + { + for (int i = u; i < v; i++) + { + // for each row in a + var row = GetRow(transposeA, i, rowsA, columnsA, a); + for (int j = 0; j < columnsB; j++) + { + var column = columnDataB[j]; + double sum = 0; + for (int ii = 0; ii < row.Length; ii++) + { + sum += row[ii]*column[ii]; + } + + c[j*rowsA + i] += alpha*sum; + } + } + }); + } + } + + static void Swap(ref int first, ref int second) + { + var prior = first; + first = second; + second = prior; + } + + /// + /// Assumes that and have already been transposed. + /// + static double[] GetRow(Transpose transpose, int rowindx, int numRows, int numCols, double[] matrix) + { + var ret = new double[numCols]; + if (transpose == Transpose.DontTranspose) + { + for (int i = 0; i < numCols; i++) + { + ret[i] = matrix[(i*numRows) + rowindx]; + } + } + else + { + Array.Copy(matrix, rowindx*numCols, ret, 0, numCols); + } + + return ret; + } + + /// + /// Assumes that and have already been transposed. + /// + static double[] GetColumn(Transpose transpose, int colindx, int numRows, int numCols, double[] matrix) + { + var ret = new double[numRows]; + if (transpose == Transpose.DontTranspose) + { + Array.Copy(matrix, colindx*numRows, ret, 0, numRows); + } + else + { + for (int i = 0; i < numRows; i++) + { + ret[i] = matrix[(i*numCols) + colindx]; + } + } + + return ret; + } + } +} diff --git a/src/Performance/LinearAlgebra/DenseVectorAdd.cs b/src/Performance/LinearAlgebra/DenseVectorAdd.cs index 66476854..c4ce38f0 100644 --- a/src/Performance/LinearAlgebra/DenseVectorAdd.cs +++ b/src/Performance/LinearAlgebra/DenseVectorAdd.cs @@ -11,103 +11,139 @@ namespace Performance.LinearAlgebra { public class DenseVectorAdd { - readonly Vector a; - readonly Vector b; + readonly int _rounds; + readonly Vector _a; + readonly Vector _b; - readonly ILinearAlgebraProvider managed = new ManagedLinearAlgebraProvider(); - readonly ILinearAlgebraProvider mkl = new MklLinearAlgebraProvider(); + readonly ILinearAlgebraProvider _managed = new ManagedLinearAlgebraProvider(); + readonly ILinearAlgebraProvider _mkl = new MklLinearAlgebraProvider(); - public DenseVectorAdd(int size) + public DenseVectorAdd(int size, int rounds) { - b = Vector.Build.Random(size); - a = Vector.Build.Random(size); + _rounds = rounds; - managed.InitializeVerify(); - Control.LinearAlgebraProvider = managed; + _b = Vector.Build.Random(size); + _a = Vector.Build.Random(size); + + _managed.InitializeVerify(); + Control.LinearAlgebraProvider = _managed; #if NATIVEMKL - mkl.InitializeVerify(); - Console.WriteLine("MklProvider: {0}", mkl); - //Control.LinearAlgebraProvider = mkl; + _mkl.InitializeVerify(); #endif } [BenchSharkTask("AddOperator")] public Vector AddOperator() { - return a + b; + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a + z; + } + return z; } [BenchSharkTask("Map2")] public Vector Map2() { - return a.Map2((u, v) => u + v, b); + var z = _b; + for (int i = 0; i < _rounds; i++) + { + z = _a.Map2((u, v) => u + v, z); + } + return z; } [BenchSharkTask("Loop")] public Vector Loop() { - var aa = ((DenseVectorStorage)a.Storage).Data; - var ab = ((DenseVectorStorage)b.Storage).Data; - var ar = new Double[aa.Length]; - for (int i = 0; i < ar.Length; i++) + var z = _b; + for (int i = 0; i < _rounds; i++) { - ar[i] = aa[i] + ab[i]; + var aa = ((DenseVectorStorage)_a.Storage).Data; + var az = ((DenseVectorStorage)z.Storage).Data; + var ar = new Double[aa.Length]; + for (int k = 0; k < ar.Length; k++) + { + ar[k] = aa[k] + az[k]; + } + z = Vector.Build.Dense(ar); } - return Vector.Build.Dense(ar); + return z; } [BenchSharkTask("ParallelLoop4096")] public Vector ParallelLoop4096() { - var aa = ((DenseVectorStorage)a.Storage).Data; - var ab = ((DenseVectorStorage)b.Storage).Data; - var ar = new Double[aa.Length]; - CommonParallel.For(0, ar.Length, 4096, (u, v) => + var z = _b; + for (int i = 0; i < _rounds; i++) { - for (int i = u; i < v; i++) + var aa = ((DenseVectorStorage)_a.Storage).Data; + var az = ((DenseVectorStorage)z.Storage).Data; + var ar = new Double[aa.Length]; + CommonParallel.For(0, ar.Length, 4096, (u, v) => { - ar[i] = aa[i] + ab[i]; - } - }); - return Vector.Build.Dense(ar); + for (int k = u; k < v; k++) + { + ar[k] = aa[k] + az[k]; + } + }); + z = Vector.Build.Dense(ar); + } + return z; } [BenchSharkTask("ParallelLoop32768")] public Vector ParallelLoop32768() { - var aa = ((DenseVectorStorage)a.Storage).Data; - var ab = ((DenseVectorStorage)b.Storage).Data; - var ar = new Double[aa.Length]; - CommonParallel.For(0, ar.Length, 32768, (u, v) => + var z = _b; + for (int i = 0; i < _rounds; i++) { - for (int i = u; i < v; i++) + var aa = ((DenseVectorStorage)_a.Storage).Data; + var az = ((DenseVectorStorage)z.Storage).Data; + var ar = new Double[aa.Length]; + CommonParallel.For(0, ar.Length, 32768, (u, v) => { - ar[i] = aa[i] + ab[i]; - } - }); - return Vector.Build.Dense(ar); + for (int k = u; k < v; k++) + { + ar[k] = aa[k] + az[k]; + } + }); + z = Vector.Build.Dense(ar); + } + return z; } [BenchSharkTask("ManagedProvider")] public Vector ManagedProvider() { - var aa = ((DenseVectorStorage)a.Storage).Data; - var ab = ((DenseVectorStorage)b.Storage).Data; - var ar = new Double[aa.Length]; - managed.AddArrays(aa, ab, ar); - return Vector.Build.Dense(ar); + var z = _b; + for (int i = 0; i < _rounds; i++) + { + var aa = ((DenseVectorStorage)_a.Storage).Data; + var az = ((DenseVectorStorage)z.Storage).Data; + var ar = new Double[aa.Length]; + _managed.AddArrays(aa, az, ar); + z = Vector.Build.Dense(ar); + } + return z; } #if NATIVEMKL [BenchSharkTask("MklProvider")] public Vector MklProvider() { - var aa = ((DenseVectorStorage)a.Storage).Data; - var ab = ((DenseVectorStorage)b.Storage).Data; - var ar = new Double[aa.Length]; - mkl.AddArrays(aa, ab, ar); - return Vector.Build.Dense(ar); + var z = _b; + for (int i = 0; i < _rounds; i++) + { + var aa = ((DenseVectorStorage)_a.Storage).Data; + var az = ((DenseVectorStorage)z.Storage).Data; + var ar = new Double[aa.Length]; + _mkl.AddArrays(aa, az, ar); + z = Vector.Build.Dense(ar); + } + return z; } #endif } diff --git a/src/Performance/Performance.csproj b/src/Performance/Performance.csproj index 9fb79646..346f3c02 100644 --- a/src/Performance/Performance.csproj +++ b/src/Performance/Performance.csproj @@ -20,6 +20,8 @@ DEBUG;TRACE prompt 4 + true + x64 pdbonly @@ -29,6 +31,7 @@ prompt 4 x64 + true @@ -60,6 +63,7 @@ + diff --git a/src/Performance/Program.cs b/src/Performance/Program.cs index cb9717da..78ae3239 100644 --- a/src/Performance/Program.cs +++ b/src/Performance/Program.cs @@ -1,6 +1,8 @@ -using System.Linq; +using System; +using System.Linq; using Binarysharp.Benchmark; using ConsoleDump; +using MathNet.Numerics.Statistics; namespace Performance { @@ -8,24 +10,43 @@ namespace Performance { public static void Main() { - Run(new LinearAlgebra.DenseVectorAdd(10000000), 10, "Large (10'000'000)"); - Run(new LinearAlgebra.DenseVectorAdd(100), 10000, "Small (100)"); - } + //Benchmark(new LinearAlgebra.DenseVectorAdd(10000000,1), 10, "Large (10'000'000) - 10x1 iterations"); + //Benchmark(new LinearAlgebra.DenseVectorAdd(100,1000), 100, "Small (100) - 100x1000 iterations"); - static void Run(uint iterations, string suffix = null) where T:new() - { - var bench = new BenchShark(); - var result = bench.EvaluateDecoratedTasks(iterations); - var label = string.IsNullOrEmpty(suffix) ? typeof (T).FullName : string.Concat(typeof (T).FullName, ": ", suffix); - result.FastestEvaluations.Select(x => new { x.Name, x.BestExecutionTime, x.AverageExecutionTime, x.WorstExecutionTime }).Dump(label); + LinearAlgebra.DenseMatrixProduct.Verify(5); + LinearAlgebra.DenseMatrixProduct.Verify(100); + Benchmark(new LinearAlgebra.DenseMatrixProduct(10,100), 100, "10 - 100x100 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(25, 100), 100, "25 - 100x100 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(50, 10), 100, "50 - 100x10 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(100, 10), 100, "100 - 100x10 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(250, 1), 10, "250 - 10x1 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(500,1), 10, "500 - 10x1 iterations"); + Benchmark(new LinearAlgebra.DenseMatrixProduct(1000,1), 2, "1000 - 2x1 iterations"); } - static void Run(object obj, uint iterations, string suffix = null) + static void Benchmark(object obj, uint iterations, string suffix = null) { - var bench = new BenchShark(); + var bench = new BenchShark(true); var result = bench.EvaluateDecoratedTasks(obj, iterations); + var results = result.FastestEvaluations.Select(x => + { + var series = x.Iterations.Select(it => (double)it.ElapsedTicks).ToArray(); + Array.Sort(series); + var summary = SortedArrayStatistics.FiveNumberSummary(series); + var ms = ArrayStatistics.MeanStandardDeviation(series); + return new { x.Name, Mean = ms.Item1, StdDev = ms.Item2, Min = summary[0], Q1 = summary[1], Median = summary[2], Q3 = summary[3], Max = summary[4] }; + }).ToArray(); + var top = results[0]; + var managed = results.Single(x => x.Name.StartsWith("Managed")); var label = string.IsNullOrEmpty(suffix) ? obj.GetType().FullName : string.Concat(obj.GetType().FullName, ": ", suffix); - result.FastestEvaluations.Select(x => new { x.Name, x.BestExecutionTime, x.AverageExecutionTime, x.WorstExecutionTime }).Dump(label); + results.Select(x => new + { + x.Name, + Mean = Math.Round(x.Mean), StdDev = Math.Round(x.StdDev), + Min = Math.Round(x.Min), Q1 = Math.Round(x.Q1), Median = Math.Round(x.Median), Q3 = Math.Round(x.Q3), Max = Math.Round(x.Max), + TopSlowdown = Math.Round(x.Median/top.Median, 2), + ManagedSpeedup = Math.Round(managed.Median/x.Median, 2) + }).Dump(label); } } } diff --git a/src/UnitTests/StatisticsTests/StatisticsTests.cs b/src/UnitTests/StatisticsTests/StatisticsTests.cs index 92748ea2..5686952b 100644 --- a/src/UnitTests/StatisticsTests/StatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/StatisticsTests.cs @@ -772,6 +772,16 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.AreEqual(0.2d, SortedArrayStatistics.Median(odd), 1e-14); } + [Test] + public void MedianOnLongConstantSequence() + { + var even = Generate.Repeat(100000, 2.0); + Assert.AreEqual(2.0,SortedArrayStatistics.Median(even), 1e-14); + + var odd = Generate.Repeat(100001, 2.0); + Assert.AreEqual(2.0, SortedArrayStatistics.Median(odd), 1e-14); + } + /// /// Validate Median/Variance/StdDev on a longer fixed-random sequence of a, /// large mean but only a very small variance, verifying the numerical stability.