LA: managed provider experimental dense multiplication variation

src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs

@@ -14,6 +14,7 @@ namespace Benchmark.LinearAlgebra
         readonly Dictionary<string, Matrix<double>> _data = new Dictionary<string, Matrix<double>>();
 
         readonly ILinearAlgebraProvider _managed;
+        readonly ILinearAlgebraProvider _managedExperimental;
         readonly ILinearAlgebraProvider _mkl;
         readonly ILinearAlgebraProvider _experimental;
 
@@ -39,10 +40,12 @@ namespace Benchmark.LinearAlgebra
 
             Control.NativeProviderPath = @"..\..\..\..\out\MKL\Windows\";
             _managed = new ManagedLinearAlgebraProvider();
+            _managedExperimental = new ManagedLinearAlgebraProvider(Variation.Experimental);
             _mkl = new MklLinearAlgebraProvider();
             _experimental = new ExperimentalProvider();
 
             _managed.InitializeVerify();
+            _managedExperimental.InitializeVerify();
             _mkl.InitializeVerify();
             _experimental.InitializeVerify();
 
@@ -55,6 +58,7 @@ namespace Benchmark.LinearAlgebra
             N = 8;
             var resultMkl = MathNet().ToRowArrays();
             var resultManaged = MathNetManaged().ToRowArrays();
+            var resultManagedExperimental = MathNetManagedExperimental().ToRowArrays();
             var resultExperimental = MathNetExperimental().ToRowArrays();
             for (int i = 0; i < 8; i++)
             {
@@ -64,6 +68,10 @@ 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]}");
@@ -91,6 +99,13 @@ namespace Benchmark.LinearAlgebra
             return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]);
         }
 
+        [Benchmark(OperationsPerInvoke = 1)]
+        public Matrix<double> MathNetManagedExperimental()
+        {
+            Control.LinearAlgebraProvider = _managedExperimental;
+            return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]);
+        }
+
         [Benchmark(OperationsPerInvoke = 1)]
         public Matrix<double> MathNetExperimental()
         {
@@ -532,7 +547,9 @@ namespace Benchmark.LinearAlgebra
                 var columnDataB = new double[columnsB][];
                 for (int i = 0; i < columnDataB.Length; i++)
                 {
-                    columnDataB[i] = GetColumn(transposeB, i, rowsB, columnsB, b);
+                    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;
@@ -600,22 +617,19 @@ namespace Benchmark.LinearAlgebra
             /// <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)
+            static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, double[] matrix, double[] column)
             {
-                var ret = new double[numRows];
                 if (transpose == Transpose.DontTranspose)
                 {
-                    Array.Copy(matrix, colindx * numRows, ret, 0, numRows);
+                    Array.Copy(matrix, colindx * numRows, column, 0, numRows);
                 }
                 else
                 {
                     for (int i = 0; i < numRows; i++)
                     {
-                        ret[i] = matrix[(i * numCols) + colindx];
+                        column[i] = matrix[(i * numCols) + colindx];
                     }
                 }
-
-                return ret;
             }
         }
     }

src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@@ -436,6 +436,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
         public virtual void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result);
+                return;
+            }
+
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (x == null)
             {
@@ -516,6 +522,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// <param name="c">The c matrix.</param>
         public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
+                return;
+            }
+
             int m; // The number of rows of matrix op(A) and of the matrix C.
             int n; // The number of columns of matrix op(B) and of the matrix C.
             int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@@ -894,6 +906,128 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
         }
 
+        public void MatrixMultiplyWithUpdateExperimental(
+                Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA,
+                Complex[] b,
+                int rowsB, int columnsB, Complex beta, Complex[] 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)
+            {
+                var swap = rowsA;
+                rowsA = columnsA;
+                columnsA = swap;
+            }
+
+            if (transposeB != Transpose.DontTranspose)
+            {
+                var swap = rowsB;
+                rowsB = columnsB;
+                columnsB = swap;
+            }
+
+            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 degenerate cases
+            if (beta == Complex.Zero)
+            {
+                Array.Clear(c, 0, c.Length);
+            }
+            else if (beta != Complex.One)
+            {
+                ScaleArray(beta, c, c);
+            }
+
+            if (alpha == Complex.Zero)
+            {
+                return;
+            }
+
+            // Extract column arrays
+            var columnDataB = new Complex[columnsB][];
+            for (int i = 0; i < columnDataB.Length; i++)
+            {
+                var column = new Complex[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 Complex[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];
+                        Complex sum = Complex.Zero;
+                        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 Complex[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];
+                            Complex sum = Complex.Zero;
+                            for (int ii = 0; ii < row.Length; ii++)
+                            {
+                                sum += row[ii] * column[ii];
+                            }
+
+                            c[j * rowsA + i] += alpha * sum;
+                        }
+                    }
+                });
+            }
+        }
+
         /// <summary>
         /// Computes the LUP factorization of A. P*A = L*U.
         /// </summary>

src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@@ -433,6 +433,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
         public virtual void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result);
+                return;
+            }
+
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (x == null)
             {
@@ -513,6 +519,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// <param name="c">The c matrix.</param>
         public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
+                return;
+            }
+
             int m; // The number of rows of matrix op(A) and of the matrix C.
             int n; // The number of columns of matrix op(B) and of the matrix C.
             int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@@ -891,6 +903,128 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
         }
 
+        public void MatrixMultiplyWithUpdateExperimental(
+                Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA,
+                Complex32[] b,
+                int rowsB, int columnsB, Complex32 beta, Complex32[] 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)
+            {
+                var swap = rowsA;
+                rowsA = columnsA;
+                columnsA = swap;
+            }
+
+            if (transposeB != Transpose.DontTranspose)
+            {
+                var swap = rowsB;
+                rowsB = columnsB;
+                columnsB = swap;
+            }
+
+            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 degenerate cases
+            if (beta == Complex32.Zero)
+            {
+                Array.Clear(c, 0, c.Length);
+            }
+            else if (beta != Complex32.One)
+            {
+                ScaleArray(beta, c, c);
+            }
+
+            if (alpha == Complex32.Zero)
+            {
+                return;
+            }
+
+            // Extract column arrays
+            var columnDataB = new Complex32[columnsB][];
+            for (int i = 0; i < columnDataB.Length; i++)
+            {
+                var column = new Complex32[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 Complex32[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];
+                        Complex32 sum = Complex32.Zero;
+                        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 Complex32[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];
+                            Complex32 sum = Complex32.Zero;
+                            for (int ii = 0; ii < row.Length; ii++)
+                            {
+                                sum += row[ii] * column[ii];
+                            }
+
+                            c[j * rowsA + i] += alpha * sum;
+                        }
+                    }
+                });
+            }
+        }
+
         /// <summary>
         /// Computes the LUP factorization of A. P*A = L*U.
         /// </summary>

src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@@ -428,6 +428,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
         public virtual void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result);
+                return;
+            }
+
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (x == null)
             {
@@ -508,6 +514,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// <param name="c">The c matrix.</param>
         public virtual 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 (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
+                return;
+            }
+
             int m; // The number of rows of matrix op(A) and of the matrix C.
             int n; // The number of columns of matrix op(B) and of the matrix C.
             int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@@ -777,6 +789,128 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
         }
 
+        public void MatrixMultiplyWithUpdateExperimental(
+                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)
+            {
+                var swap = rowsA;
+                rowsA = columnsA;
+                columnsA = swap;
+            }
+
+            if (transposeB != Transpose.DontTranspose)
+            {
+                var swap = rowsB;
+                rowsB = columnsB;
+                columnsB = swap;
+            }
+
+            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 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;
+                        }
+                    }
+                });
+            }
+        }
+
         /// <summary>
         /// Computes the LUP factorization of A. P*A = L*U.
         /// </summary>
@@ -1329,7 +1463,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
                     }
                 });
 
-            var work = columnsR > rowsR ? new double[rowsR * rowsR] : new double[rowsR * columnsR]; 
+            var work = columnsR > rowsR ? new double[rowsR * rowsR] : new double[rowsR * columnsR];
             var minmn = Math.Min(rowsR, columnsR);
             for (var i = 0; i < minmn; i++)
             {
@@ -1383,7 +1517,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
 
             var work = new double[rowsA*columnsA];
-  
+
             var minmn = Math.Min(rowsA, columnsA);
             for (var i = 0; i < minmn; i++)
             {

src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@@ -428,6 +428,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
         public virtual void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result);
+                return;
+            }
+
             // First check some basic requirement on the parameters of the matrix multiplication.
             if (x == null)
             {
@@ -508,6 +514,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         /// <param name="c">The c matrix.</param>
         public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c)
         {
+            if (_variation == Variation.Experimental)
+            {
+                MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
+                return;
+            }
+
             int m; // The number of rows of matrix op(A) and of the matrix C.
             int n; // The number of columns of matrix op(B) and of the matrix C.
             int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@@ -777,6 +789,128 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
         }
 
+        public void MatrixMultiplyWithUpdateExperimental(
+                Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA,
+                float[] b,
+                int rowsB, int columnsB, float beta, float[] 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)
+            {
+                var swap = rowsA;
+                rowsA = columnsA;
+                columnsA = swap;
+            }
+
+            if (transposeB != Transpose.DontTranspose)
+            {
+                var swap = rowsB;
+                rowsB = columnsB;
+                columnsB = swap;
+            }
+
+            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 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 float[columnsB][];
+            for (int i = 0; i < columnDataB.Length; i++)
+            {
+                var column = new float[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 float[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];
+                        float 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 float[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];
+                            float sum = 0;
+                            for (int ii = 0; ii < row.Length; ii++)
+                            {
+                                sum += row[ii] * column[ii];
+                            }
+
+                            c[j * rowsA + i] += alpha * sum;
+                        }
+                    }
+                });
+            }
+        }
+
         /// <summary>
         /// Computes the LUP factorization of A. P*A = L*U.
         /// </summary>
@@ -1778,7 +1912,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
             }
 
             var work = new float[rowsA];
- 
+
             const int maxiter = 1000;
 
             var e = new float[columnsA];

src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs

@@ -27,13 +27,33 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System;
+
 namespace MathNet.Numerics.Providers.LinearAlgebra
 {
+    internal enum Variation
+    {
+        Original,
+        Experimental
+    }
+
     /// <summary>
     /// The managed linear algebra provider.
     /// </summary>
     public partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider
     {
+        private Variation _variation;
+
+        public ManagedLinearAlgebraProvider()
+        {
+            _variation = Variation.Original;
+        }
+
+        internal ManagedLinearAlgebraProvider(Variation variation)
+        {
+            _variation = variation;
+        }
+
         /// <summary>
         /// Try to find out whether the provider is available, at least in principle.
         /// Verification may still fail if available, but it will certainly fail if unavailable.
@@ -54,5 +74,42 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
         {
             return "Managed";
         }
+
+
+        /// <summary>
+        /// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
+        /// </summary>
+        protected static void GetRow<T>(Transpose transpose, int rowindx, int numRows, int numCols, T[] matrix, T[] row)
+        {
+            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);
+            }
+        }
+
+        /// <summary>
+        /// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
+        /// </summary>
+        protected static void GetColumn<T>(Transpose transpose, int colindx, int numRows, int numCols, T[] matrix, T[] 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];
+                }
+            }
+        }
     }
 }