Perf: matrix product bench with experimental providers

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())

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())

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)

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)

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<double> _a;
+        readonly Matrix<double> _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<double>.Build.Random(size, size);
+            _a = Matrix<double>.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<double> ManagedProvider()
+        {
+            Control.LinearAlgebraProvider = _managed;
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                z = _a*z;
+            }
+            return z;
+        }
+
+        [BenchSharkTask("MklProvider")]
+        public Matrix<double> MklProvider()
+        {
+            Control.LinearAlgebraProvider = _mkl;
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                z = _a*z;
+            }
+            return z;
+        }
+
+        [BenchSharkTask("SafeProvider")]
+        public Matrix<double> SafeProvider()
+        {
+            Control.LinearAlgebraProvider = _safeProvider;
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                z = _a*z;
+            }
+            return z;
+        }
+
+        [BenchSharkTask("UnsafeProvider")]
+        public Matrix<double> UnsafeProvider()
+        {
+            Control.LinearAlgebraProvider = _unsafeProvider;
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                z = _a*z;
+            }
+            return z;
+        }
+
+        [BenchSharkTask("ExperimentalProvider")]
+        public Matrix<double> 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;
+        }
+
+        /// <summary>
+        /// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
+        /// </summary>
+        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;
+        }
+
+        /// <summary>
+        /// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
+        /// </summary>
+        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;
+        }
+    }
+}

src/Performance/LinearAlgebra/DenseVectorAdd.cs

@@ -11,103 +11,139 @@ namespace Performance.LinearAlgebra
 {
     public class DenseVectorAdd
     {
-        readonly Vector<double> a;
-        readonly Vector<double> b;
+        readonly int _rounds;
+        readonly Vector<double> _a;
+        readonly Vector<double> _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<double>.Build.Random(size);
-            a = Vector<double>.Build.Random(size);
+            _rounds = rounds;
 
-            managed.InitializeVerify();
-            Control.LinearAlgebraProvider = managed;
+            _b = Vector<double>.Build.Random(size);
+            _a = Vector<double>.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<double> AddOperator()
         {
-            return a + b;
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                z = _a + z;
+            }
+            return z;
         }
 
         [BenchSharkTask("Map2")]
         public Vector<double> 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<double> Loop()
         {
-            var aa = ((DenseVectorStorage<double>)a.Storage).Data;
-            var ab = ((DenseVectorStorage<double>)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<double>)_a.Storage).Data;
+                var az = ((DenseVectorStorage<double>)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<double>.Build.Dense(ar);
             }
-            return Vector<double>.Build.Dense(ar);
+            return z;
         }
 
         [BenchSharkTask("ParallelLoop4096")]
         public Vector<double> ParallelLoop4096()
         {
-            var aa = ((DenseVectorStorage<double>)a.Storage).Data;
-            var ab = ((DenseVectorStorage<double>)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<double>)_a.Storage).Data;
+                var az = ((DenseVectorStorage<double>)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<double>.Build.Dense(ar);
+                    for (int k = u; k < v; k++)
+                    {
+                        ar[k] = aa[k] + az[k];
+                    }
+                });
+                z = Vector<double>.Build.Dense(ar);
+            }
+            return z;
         }
 
         [BenchSharkTask("ParallelLoop32768")]
         public Vector<double> ParallelLoop32768()
         {
-            var aa = ((DenseVectorStorage<double>)a.Storage).Data;
-            var ab = ((DenseVectorStorage<double>)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<double>)_a.Storage).Data;
+                var az = ((DenseVectorStorage<double>)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<double>.Build.Dense(ar);
+                    for (int k = u; k < v; k++)
+                    {
+                        ar[k] = aa[k] + az[k];
+                    }
+                });
+                z = Vector<double>.Build.Dense(ar);
+            }
+            return z;
         }
 
         [BenchSharkTask("ManagedProvider")]
         public Vector<double> ManagedProvider()
         {
-            var aa = ((DenseVectorStorage<double>)a.Storage).Data;
-            var ab = ((DenseVectorStorage<double>)b.Storage).Data;
-            var ar = new Double[aa.Length];
-            managed.AddArrays(aa, ab, ar);
-            return Vector<double>.Build.Dense(ar);
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                var aa = ((DenseVectorStorage<double>)_a.Storage).Data;
+                var az = ((DenseVectorStorage<double>)z.Storage).Data;
+                var ar = new Double[aa.Length];
+                _managed.AddArrays(aa, az, ar);
+                z = Vector<double>.Build.Dense(ar);
+            }
+            return z;
         }
 
 #if NATIVEMKL
         [BenchSharkTask("MklProvider")]
         public Vector<double> MklProvider()
         {
-            var aa = ((DenseVectorStorage<double>)a.Storage).Data;
-            var ab = ((DenseVectorStorage<double>)b.Storage).Data;
-            var ar = new Double[aa.Length];
-            mkl.AddArrays(aa, ab, ar);
-            return Vector<double>.Build.Dense(ar);
+            var z = _b;
+            for (int i = 0; i < _rounds; i++)
+            {
+                var aa = ((DenseVectorStorage<double>)_a.Storage).Data;
+                var az = ((DenseVectorStorage<double>)z.Storage).Data;
+                var ar = new Double[aa.Length];
+                _mkl.AddArrays(aa, az, ar);
+                z = Vector<double>.Build.Dense(ar);
+            }
+            return z;
         }
 #endif
     }

src/Performance/Performance.csproj

@@ -20,6 +20,8 @@
     <DefineConstants>DEBUG;TRACE</DefineConstants>
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
+    <AllowUnsafeBlocks>true</AllowUnsafeBlocks>
+    <PlatformTarget>x64</PlatformTarget>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
@@ -29,6 +31,7 @@
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
     <PlatformTarget>x64</PlatformTarget>
+    <AllowUnsafeBlocks>true</AllowUnsafeBlocks>
   </PropertyGroup>
   <PropertyGroup>
     <StartupObject />
@@ -60,6 +63,7 @@
   </ItemGroup>
   <ItemGroup>
     <Compile Include="LinearAlgebra\DenseVectorAdd.cs" />
+    <Compile Include="LinearAlgebra\DenseMatrixProduct.cs" />
     <Compile Include="Program.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
   </ItemGroup>

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<T>(uint iterations, string suffix = null) where T:new()
-        {
-            var bench = new BenchShark();
-            var result = bench.EvaluateDecoratedTasks<T>(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);
         }
     }
 }

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);
+        }
+
         /// <summary>
         /// Validate Median/Variance/StdDev on a longer fixed-random sequence of a,
         /// large mean but only a very small variance, verifying the numerical stability.