diff --git a/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs b/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs index 40aa8e23..56547af7 100644 --- a/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs +++ b/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs @@ -13,10 +13,9 @@ namespace Benchmark.LinearAlgebra { readonly Dictionary> _data = new Dictionary>(); - readonly ILinearAlgebraProvider _managed; - readonly ILinearAlgebraProvider _managedExperimental; - readonly ILinearAlgebraProvider _mkl; - readonly ILinearAlgebraProvider _experimental; + readonly ILinearAlgebraProvider _mathnetMkl; + readonly ILinearAlgebraProvider _mathnetManaged; + readonly ILinearAlgebraProvider _mathnetExperimental; [Params(8, 64, 128)] public int M { get; set; } @@ -39,26 +38,21 @@ namespace Benchmark.LinearAlgebra } Control.NativeProviderPath = @"..\..\..\..\out\MKL\Windows\"; - _managed = new ManagedLinearAlgebraProvider(); - _managedExperimental = new ManagedLinearAlgebraProvider(Variation.Experimental); - _mkl = new MklLinearAlgebraProvider(); - _experimental = new ExperimentalProvider(); + _mathnetMkl = new MklLinearAlgebraProvider(); + _mathnetManaged = new ManagedLinearAlgebraProvider(); + _mathnetExperimental = new ManagedLinearAlgebraProvider(Variation.Experimental); - _managed.InitializeVerify(); - _managedExperimental.InitializeVerify(); - _mkl.InitializeVerify(); - _experimental.InitializeVerify(); - - //Verify(); + _mathnetMkl.InitializeVerify(); + _mathnetManaged.InitializeVerify(); + _mathnetExperimental.InitializeVerify(); } - private void Verify() + public void Verify() { M = 8; N = 8; - var resultMkl = MathNet().ToRowArrays(); + var resultMkl = MathNetMKL().ToRowArrays(); var resultManaged = MathNetManaged().ToRowArrays(); - var resultManagedExperimental = MathNetManagedExperimental().ToRowArrays(); var resultExperimental = MathNetExperimental().ToRowArrays(); for (int i = 0; i < 8; i++) { @@ -68,10 +62,6 @@ namespace Benchmark.LinearAlgebra { throw new Exception($"Managed [{i}][{j}] {resultManaged[i][j]} != {resultMkl[i][j]}"); } - if (!resultMkl[i][j].AlmostEqual(resultManagedExperimental[i][j], 1e-14)) - { - throw new Exception($"ManagedExperimental [{i}][{j}] {resultManagedExperimental[i][j]} != {resultMkl[i][j]}"); - } if (!resultMkl[i][j].AlmostEqual(resultExperimental[i][j], 1e-14)) { throw new Exception($"Experimental [{i}][{j}] {resultExperimental[i][j]} != {resultMkl[i][j]}"); @@ -80,557 +70,40 @@ namespace Benchmark.LinearAlgebra } } - [Setup] - public void Setup() - { - } - - [Benchmark(OperationsPerInvoke = 1)] - public Matrix MathNet() - { - Control.LinearAlgebraProvider = _mkl; - return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); - } - - [Benchmark(OperationsPerInvoke = 1)] - public Matrix MathNetManaged() - { - Control.LinearAlgebraProvider = _managed; - return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); - } - [Benchmark(OperationsPerInvoke = 1)] - public Matrix MathNetManagedExperimental() - { - Control.LinearAlgebraProvider = _managedExperimental; - return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); - } - - [Benchmark(OperationsPerInvoke = 1)] - public Matrix MathNetExperimental() - { - Control.LinearAlgebraProvider = _experimental; - return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); - } - - public class SafeProvider : ManagedLinearAlgebraProvider + public Matrix MathNetMKL() { - public override void MatrixMultiply( - double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) - { - if (rowsX + columnsY <= Control.MaxDegreeOfParallelism) - { - 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) + Control.LinearAlgebraProvider = _mathnetMkl; + if (M != N) { - 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); + return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); } - 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.MaxDegreeOfParallelism) - { - 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); - } - } + return _data[Key(M, N)]*_data[Key(M, N)]; } - public unsafe class UnsafeProvider : ManagedLinearAlgebraProvider + [Benchmark(OperationsPerInvoke = 1)] + public Matrix MathNetManaged() { - public override void MatrixMultiply( - double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + Control.LinearAlgebraProvider = _mathnetManaged; + if (M != N) { - 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); + return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); } - 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); - } - } + return _data[Key(M, N)] * _data[Key(M, N)]; } - public class ExperimentalProvider : ManagedLinearAlgebraProvider + [Benchmark(OperationsPerInvoke = 1)] + public Matrix MathNetExperimental() { - 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(nameof(a)); - } - - if (b == null) - { - throw new ArgumentNullException(nameof(b)); - } - - if (c == null) - { - throw new ArgumentNullException(nameof(c)); - } - - if (transposeA != Transpose.DontTranspose) - { - var swap = rowsA; - rowsA = columnsA; - columnsA = swap; - } - - if (transposeB != Transpose.DontTranspose) - { - var swap = rowsB; - rowsB = columnsB; - columnsB = swap; - } - - if (columnsA != rowsB) - { - throw new ArgumentOutOfRangeException($"columnsA ({columnsA}) != rowsB ({rowsB})"); - } - - if (rowsA * columnsA != a.Length) - { - throw new ArgumentOutOfRangeException($"rowsA ({rowsA}) * columnsA ({columnsA}) != a.Length ({a.Length})"); - } - - if (rowsB * columnsB != b.Length) - { - throw new ArgumentOutOfRangeException($"rowsB ({rowsB}) * columnsB ({columnsB}) != b.Length ({b.Length})"); - } - - if (rowsA * columnsB != c.Length) - { - throw new ArgumentOutOfRangeException($"rowsA ({rowsA}) * columnsB ({columnsB}) != c.Length ({c.Length})"); - } - - // handle 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++) - { - var column = new double[rowsB]; - GetColumn(transposeB, i, rowsB, columnsB, b, column); - columnDataB[i] = column; - } - - var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; - if (shouldNotParallelize) - { - var row = new double[columnsA]; - for (int i = 0; i < rowsA; i++) - { - GetRow(transposeA, i, rowsA, columnsA, a, row); - 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) => - { - var row = new double[columnsA]; - for (int i = u; i < v; i++) - { - GetRow(transposeA, i, rowsA, columnsA, a, row); - 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; - } - } - }); - } - } - - /// - /// Assumes that and have already been transposed. - /// - static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, double[] matrix, double[] row) + Control.LinearAlgebraProvider = _mathnetExperimental; + if (M != N) { - if (transpose == Transpose.DontTranspose) - { - for (int i = 0; i < numCols; i++) - { - row[i] = matrix[(i * numRows) + rowindx]; - } - } - else - { - Array.Copy(matrix, rowindx * numCols, row, 0, numCols); - } + return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); } - /// - /// Assumes that and have already been transposed. - /// - static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, double[] matrix, double[] column) - { - if (transpose == Transpose.DontTranspose) - { - Array.Copy(matrix, colindx * numRows, column, 0, numRows); - } - else - { - for (int i = 0; i < numRows; i++) - { - column[i] = matrix[(i * numCols) + colindx]; - } - } - } + return _data[Key(M, N)] * _data[Key(M, N)]; } } } diff --git a/src/Benchmark/Program.cs b/src/Benchmark/Program.cs index 73f8abc1..5ab4efb3 100644 --- a/src/Benchmark/Program.cs +++ b/src/Benchmark/Program.cs @@ -15,8 +15,12 @@ namespace Benchmark Console.WriteLine("Linear Algebra: " + Control.LinearAlgebraProvider); Console.WriteLine("FFT: " + Control.FourierTransformProvider); + var subject = new LinearAlgebra.DenseMatrixProduct(); + subject.Verify(); + Console.WriteLine("Verified."); + var config = ManualConfig.Create(DefaultConfig.Instance); - config.Add(Job.RyuJitX64, Job.LegacyJitX64, Job.LegacyJitX86); + config.Add(Job.RyuJitX64, Job.LegacyJitX86); //config.Add(new MemoryDiagnoser()); BenchmarkRunner.Run(config);