diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs index 2d41f1fe..525d75ee 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs @@ -468,480 +468,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { - throw new ArgumentNullException("x"); + throw new ArgumentNullException("a"); } if (y == null) { - throw new ArgumentNullException("y"); + throw new ArgumentNullException("b"); } if (result == null) { - throw new ArgumentNullException("result"); - } - - if (rowsX*columnsX != x.Length) - { - throw new ArgumentException("x.Length != xRows * xColumns"); - } - - if (rowsY*columnsY != y.Length) - { - throw new ArgumentException("y.Length != yRows * yColumns"); + throw new ArgumentNullException("c"); } if (columnsX != rowsY) { - throw new ArgumentException("xColumns != yRows"); + throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY)); } - if (rowsX*columnsY != result.Length) + if (rowsX * columnsX != x.Length) { - throw new ArgumentException("xRows * yColumns != result.Length"); + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length)); } - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - Complex[] xdata; - if (ReferenceEquals(x, result)) + if (rowsY * columnsY != y.Length) { - xdata = (Complex[]) x.Clone(); - } - else - { - xdata = x; + throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length)); } - Complex[] ydata; - if (ReferenceEquals(y, result)) - { - ydata = (Complex[]) y.Clone(); - } - else + if (rowsX * columnsY != result.Length) { - ydata = y; + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length)); } + // handle degenerate cases Array.Clear(result, 0, result.Length); - CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true); - } - - /// - /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c - /// - /// How to transpose the matrix. - /// How to transpose the matrix. - /// The value to scale matrix. - /// The a matrix. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The b matrix - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The value to scale the matrix. - /// The c matrix. - 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). - - // First check some basic requirement on the parameters of the matrix multiplication. - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if ((int) transposeA > 111 && (int) transposeB > 111) - { - if (rowsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = rowsB; - k = rowsA; - } - else if ((int) transposeA > 111) - { - if (rowsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = columnsB; - k = rowsA; - } - else if ((int) transposeB > 111) - { - if (columnsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = rowsB; - k = columnsA; - } - else - { - if (columnsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = columnsB; - k = columnsA; - } - - if (alpha.IsZero() && beta.IsZero()) - { - Array.Clear(c, 0, c.Length); - return; - } - - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - Complex[] adata; - if (ReferenceEquals(a, c)) - { - adata = (Complex[]) a.Clone(); - } - else - { - adata = a; - } - - Complex[] bdata; - if (ReferenceEquals(b, c)) - { - bdata = (Complex[]) b.Clone(); - } - else - { - bdata = b; - } - - if (beta.IsZero()) - { - Array.Clear(c, 0, c.Length); - } - else if (!beta.IsOne()) - { - ScaleArray(beta, c, c); - } - - if (alpha.IsZero()) + // Extract column arrays + var columnDataB = new Complex[columnsY][]; + for (int i = 0; i < columnDataB.Length; i++) { - return; + var column = new Complex[rowsY]; + GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column); + columnDataB[i] = column; } - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true); - } - - /// - /// Cache-Oblivious Matrix Multiplication - /// - /// if set to true transpose matrix A. - /// if set to true transpose matrix B. - /// The value to scale the matrix A with. - /// The matrix A. - /// Row-shift of the left matrix - /// Column-shift of the left matrix - /// The matrix B. - /// Row-shift of the right matrix - /// Column-shift of the right matrix - /// The matrix C. - /// Row-shift of the result matrix - /// Column-shift of the result matrix - /// The number of rows of matrix op(A) and of the matrix C. - /// The number of columns of matrix op(B) and of the matrix C. - /// The number of columns of matrix op(A) and the rows of the matrix op(B). - /// The constant number of rows of matrix op(A) and of the matrix C. - /// The constant number of columns of matrix op(B) and of the matrix C. - /// The constant number of columns of matrix op(A) and the rows of the matrix op(B). - /// Indicates if this is the first recursion. - static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] matrixA, int shiftArow, int shiftAcol, Complex[] matrixB, int shiftBrow, int shiftBcol, Complex[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first) - { - if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1) + var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; + if (shouldNotParallelize) { - if ((int) transposeA > 111 && (int) transposeB > 111) + var row = new Complex[columnsX]; + for (int i = 0; i < rowsX; i++) { - if ((int) transposeA > 112 && (int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeA > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - } - else if ((int) transposeA > 111) - { - if ((int) transposeA > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - } - else if ((int) transposeB > 111) - { - if ((int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) + var col = columnDataB[j]; + Complex sum = Complex.Zero; + for (int ii = 0; ii < row.Length; ii++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } + sum += row[ii] * col[ii]; } - } - else - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } + result[j * rowsX + i] += Complex.One * sum; } } - else + } + else + { + CommonParallel.For(0, rowsX, 1, (u, v) => { - for (var m1 = 0; m1 < m; m1++) + var row = new Complex[columnsX]; + for (int i = u; i < v; i++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matBcolPos = n1 + shiftBcol; - var sum = Complex.Zero; - for (var k1 = 0; k1 < k; ++k1) + var column = columnDataB[j]; + Complex sum = Complex.Zero; + for (int ii = 0; ii < row.Length; ii++) { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; + sum += row[ii] * column[ii]; } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + result[j * rowsX + i] += Complex.One * sum; } } - } - } - else - { - // divide and conquer - int m2 = m/2, n2 = n/2, k2 = k/2; - - if (first) - { - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false)); - - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false)); - } - else - { - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false); - } + }); } } - 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) + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + 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 (a == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs index 013d5d7e..6a8d7603 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs @@ -470,480 +470,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { - throw new ArgumentNullException("x"); + throw new ArgumentNullException("a"); } if (y == null) { - throw new ArgumentNullException("y"); + throw new ArgumentNullException("b"); } if (result == null) { - throw new ArgumentNullException("result"); - } - - if (rowsX*columnsX != x.Length) - { - throw new ArgumentException("x.Length != xRows * xColumns"); - } - - if (rowsY*columnsY != y.Length) - { - throw new ArgumentException("y.Length != yRows * yColumns"); + throw new ArgumentNullException("c"); } if (columnsX != rowsY) { - throw new ArgumentException("xColumns != yRows"); + throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY)); } - if (rowsX*columnsY != result.Length) + if (rowsX * columnsX != x.Length) { - throw new ArgumentException("xRows * yColumns != result.Length"); + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length)); } - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - Complex32[] xdata; - if (ReferenceEquals(x, result)) + if (rowsY * columnsY != y.Length) { - xdata = (Complex32[]) x.Clone(); - } - else - { - xdata = x; + throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length)); } - Complex32[] ydata; - if (ReferenceEquals(y, result)) - { - ydata = (Complex32[]) y.Clone(); - } - else + if (rowsX * columnsY != result.Length) { - ydata = y; + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length)); } + // handle degenerate cases Array.Clear(result, 0, result.Length); - CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true); - } - - /// - /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c - /// - /// How to transpose the matrix. - /// How to transpose the matrix. - /// The value to scale matrix. - /// The a matrix. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The b matrix - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The value to scale the matrix. - /// The c matrix. - 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). - - // First check some basic requirement on the parameters of the matrix multiplication. - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if ((int) transposeA > 111 && (int) transposeB > 111) - { - if (rowsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = rowsB; - k = rowsA; - } - else if ((int) transposeA > 111) - { - if (rowsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = columnsB; - k = rowsA; - } - else if ((int) transposeB > 111) - { - if (columnsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = rowsB; - k = columnsA; - } - else - { - if (columnsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = columnsB; - k = columnsA; - } - - if (alpha.IsZero() && beta.IsZero()) - { - Array.Clear(c, 0, c.Length); - return; - } - - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - Complex32[] adata; - if (ReferenceEquals(a, c)) - { - adata = (Complex32[]) a.Clone(); - } - else - { - adata = a; - } - - Complex32[] bdata; - if (ReferenceEquals(b, c)) - { - bdata = (Complex32[]) b.Clone(); - } - else - { - bdata = b; - } - - if (beta.IsZero()) - { - Array.Clear(c, 0, c.Length); - } - else if (!beta.IsOne()) - { - ScaleArray(beta, c, c); - } - - if (alpha.IsZero()) + // Extract column arrays + var columnDataB = new Complex32[columnsY][]; + for (int i = 0; i < columnDataB.Length; i++) { - return; + var column = new Complex32[rowsY]; + GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column); + columnDataB[i] = column; } - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true); - } - - /// - /// Cache-Oblivious Matrix Multiplication - /// - /// if set to true transpose matrix A. - /// if set to true transpose matrix B. - /// The value to scale the matrix A with. - /// The matrix A. - /// Row-shift of the left matrix - /// Column-shift of the left matrix - /// The matrix B. - /// Row-shift of the right matrix - /// Column-shift of the right matrix - /// The matrix C. - /// Row-shift of the result matrix - /// Column-shift of the result matrix - /// The number of rows of matrix op(A) and of the matrix C. - /// The number of columns of matrix op(B) and of the matrix C. - /// The number of columns of matrix op(A) and the rows of the matrix op(B). - /// The constant number of rows of matrix op(A) and of the matrix C. - /// The constant number of columns of matrix op(B) and of the matrix C. - /// The constant number of columns of matrix op(A) and the rows of the matrix op(B). - /// Indicates if this is the first recursion. - static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] matrixA, int shiftArow, int shiftAcol, Complex32[] matrixB, int shiftBrow, int shiftBcol, Complex32[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first) - { - if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1) + var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; + if (shouldNotParallelize) { - if ((int) transposeA > 111 && (int) transposeB > 111) + var row = new Complex32[columnsX]; + for (int i = 0; i < rowsX; i++) { - if ((int) transposeA > 112 && (int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeA > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - } - else if ((int) transposeA > 111) - { - if ((int) transposeA > 112) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - } - else if ((int) transposeB > 111) - { - if ((int) transposeB > 112) - { - for (var m1 = 0; m1 < m; m1++) + var col = columnDataB[j]; + Complex32 sum = Complex32.Zero; + for (int ii = 0; ii < row.Length; ii++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate(); - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } + sum += row[ii] * col[ii]; } - } - else - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } + result[j * rowsX + i] += Complex32.One * sum; } } - else + } + else + { + CommonParallel.For(0, rowsX, 1, (u, v) => { - for (var m1 = 0; m1 < m; m1++) + var row = new Complex32[columnsX]; + for (int i = u; i < v; i++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matBcolPos = n1 + shiftBcol; - var sum = Complex32.Zero; - for (var k1 = 0; k1 < k; ++k1) + var column = columnDataB[j]; + Complex32 sum = Complex32.Zero; + for (int ii = 0; ii < row.Length; ii++) { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; + sum += row[ii] * column[ii]; } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + result[j * rowsX + i] += Complex32.One * sum; } } - } - } - else - { - // divide and conquer - int m2 = m/2, n2 = n/2, k2 = k/2; - - if (first) - { - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false)); - - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false)); - } - else - { - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false); - } + }); } } - 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) + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + 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 (a == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs index 2a79ca49..1f29674b 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs @@ -465,371 +465,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { - throw new ArgumentNullException("x"); + throw new ArgumentNullException("a"); } if (y == null) { - throw new ArgumentNullException("y"); + throw new ArgumentNullException("b"); } if (result == null) { - throw new ArgumentNullException("result"); - } - - if (rowsX * columnsX != x.Length) - { - throw new ArgumentException("x.Length != xRows * xColumns"); - } - - if (rowsY * columnsY != y.Length) - { - throw new ArgumentException("y.Length != yRows * yColumns"); + throw new ArgumentNullException("c"); } if (columnsX != rowsY) { - throw new ArgumentException("xColumns != yRows"); + throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY)); } - if (rowsX * columnsY != result.Length) + if (rowsX * columnsX != x.Length) { - throw new ArgumentException("xRows * yColumns != result.Length"); + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length)); } - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - double[] xdata; - if (ReferenceEquals(x, result)) - { - xdata = (double[])x.Clone(); - } - else + if (rowsY * columnsY != y.Length) { - xdata = x; + throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length)); } - double[] ydata; - if (ReferenceEquals(y, result)) - { - ydata = (double[])y.Clone(); - } - else + if (rowsX * columnsY != result.Length) { - ydata = y; + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length)); } + // handle degenerate cases Array.Clear(result, 0, result.Length); - CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true); - } - - /// - /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c - /// - /// How to transpose the matrix. - /// How to transpose the matrix. - /// The value to scale matrix. - /// The a matrix. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The b matrix - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The value to scale the matrix. - /// The c matrix. - 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). - - // First check some basic requirement on the parameters of the matrix multiplication. - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if ((int)transposeA > 111 && (int)transposeB > 111) - { - if (rowsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA * rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = rowsB; - k = rowsA; - } - else if ((int)transposeA > 111) - { - if (rowsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA * columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = columnsB; - k = rowsA; - } - else if ((int)transposeB > 111) - { - if (columnsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA * rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = rowsB; - k = columnsA; - } - else - { - if (columnsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA * columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = columnsB; - k = columnsA; - } - - if (alpha == 0.0 && beta == 0.0) - { - Array.Clear(c, 0, c.Length); - return; - } - - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - double[] adata; - if (ReferenceEquals(a, c)) - { - adata = (double[])a.Clone(); - } - else - { - adata = a; - } - - double[] bdata; - if (ReferenceEquals(b, c)) - { - bdata = (double[])b.Clone(); - } - else - { - bdata = b; - } - - if (beta == 0.0) - { - Array.Clear(c, 0, c.Length); - } - else if (beta != 1.0) - { - ScaleArray(beta, c, c); - } - - if (alpha == 0.0) + // Extract column arrays + var columnDataB = new double[columnsY][]; + for (int i = 0; i < columnDataB.Length; i++) { - return; + var column = new double[rowsY]; + GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column); + columnDataB[i] = column; } - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true); - } - - /// - /// Cache-Oblivious Matrix Multiplication - /// - /// if set to true transpose matrix A. - /// if set to true transpose matrix B. - /// The value to scale the matrix A with. - /// The matrix A. - /// Row-shift of the left matrix - /// Column-shift of the left matrix - /// The matrix B. - /// Row-shift of the right matrix - /// Column-shift of the right matrix - /// The matrix C. - /// Row-shift of the result matrix - /// Column-shift of the result matrix - /// The number of rows of matrix op(A) and of the matrix C. - /// The number of columns of matrix op(B) and of the matrix C. - /// The number of columns of matrix op(A) and the rows of the matrix op(B). - /// The constant number of rows of matrix op(A) and of the matrix C. - /// The constant number of columns of matrix op(B) and of the matrix C. - /// The constant number of columns of matrix op(A) and the rows of the matrix op(B). - /// Indicates if this is the first recursion. - static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, double alpha, 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, bool first) - { - if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1) + var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; + if (shouldNotParallelize) { - if ((int) transposeA > 111 && (int) transposeB > 111) + var row = new double[columnsX]; + for (int i = 0; i < rowsX; i++) { - for (var m1 = 0; m1 < m; m1++) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - double sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeA > 111) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + var col = columnDataB[j]; + double sum = 0; + for (int ii = 0; ii < row.Length; ii++) { - var matBcolPos = n1 + shiftBcol; - double sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + sum += row[ii] * col[ii]; } - } - } - else if ((int) transposeB > 111) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - double sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } + result[j * rowsX + i] += 1.0 * sum; } } - else + } + else + { + CommonParallel.For(0, rowsX, 1, (u, v) => { - for (var m1 = 0; m1 < m; m1++) + var row = new double[columnsX]; + for (int i = u; i < v; i++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matBcolPos = n1 + shiftBcol; + var column = columnDataB[j]; double sum = 0; - for (var k1 = 0; k1 < k; ++k1) + for (int ii = 0; ii < row.Length; ii++) { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; + sum += row[ii] * column[ii]; } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + result[j * rowsX + i] += 1.0 * sum; } } - } - } - else - { - // divide and conquer - int m2 = m/2, n2 = n/2, k2 = k/2; - - if (first) - { - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false)); - - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false)); - } - else - { - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false); - } + }); } } - 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) + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + 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 (a == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs index 136f848d..1e8fbd3c 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs @@ -465,371 +465,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { - throw new ArgumentNullException("x"); + throw new ArgumentNullException("a"); } if (y == null) { - throw new ArgumentNullException("y"); + throw new ArgumentNullException("b"); } if (result == null) { - throw new ArgumentNullException("result"); - } - - if (rowsX*columnsX != x.Length) - { - throw new ArgumentException("x.Length != xRows * xColumns"); - } - - if (rowsY*columnsY != y.Length) - { - throw new ArgumentException("y.Length != yRows * yColumns"); + throw new ArgumentNullException("c"); } if (columnsX != rowsY) { - throw new ArgumentException("xColumns != yRows"); + throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY)); } - if (rowsX*columnsY != result.Length) + if (rowsX * columnsX != x.Length) { - throw new ArgumentException("xRows * yColumns != result.Length"); + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length)); } - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - float[] xdata; - if (ReferenceEquals(x, result)) - { - xdata = (float[]) x.Clone(); - } - else + if (rowsY * columnsY != y.Length) { - xdata = x; + throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length)); } - float[] ydata; - if (ReferenceEquals(y, result)) - { - ydata = (float[]) y.Clone(); - } - else + if (rowsX * columnsY != result.Length) { - ydata = y; + throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length)); } + // handle degenerate cases Array.Clear(result, 0, result.Length); - CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true); - } - - /// - /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c - /// - /// How to transpose the matrix. - /// How to transpose the matrix. - /// The value to scale matrix. - /// The a matrix. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The b matrix - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The value to scale the matrix. - /// The c matrix. - 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). - - // First check some basic requirement on the parameters of the matrix multiplication. - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if ((int) transposeA > 111 && (int) transposeB > 111) - { - if (rowsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = rowsB; - k = rowsA; - } - else if ((int) transposeA > 111) - { - if (rowsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (columnsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = columnsA; - n = columnsB; - k = rowsA; - } - else if ((int) transposeB > 111) - { - if (columnsA != columnsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*rowsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = rowsB; - k = columnsA; - } - else - { - if (columnsA != rowsB) - { - throw new ArgumentOutOfRangeException(); - } - - if (rowsA*columnsB != c.Length) - { - throw new ArgumentOutOfRangeException(); - } - - m = rowsA; - n = columnsB; - k = columnsA; - } - - if (alpha == 0.0 && beta == 0.0) - { - Array.Clear(c, 0, c.Length); - return; - } - - // Check whether we will be overwriting any of our inputs and make copies if necessary. - // TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory - // as result, we can do it on a row wise basis. We should investigate this. - float[] adata; - if (ReferenceEquals(a, c)) - { - adata = (float[]) a.Clone(); - } - else - { - adata = a; - } - - float[] bdata; - if (ReferenceEquals(b, c)) - { - bdata = (float[]) b.Clone(); - } - else - { - bdata = b; - } - - if (beta == 0.0f) - { - Array.Clear(c, 0, c.Length); - } - else if (beta != 1.0f) - { - ScaleArray(beta, c, c); - } - - if (alpha == 0.0f) + // Extract column arrays + var columnDataB = new float[columnsY][]; + for (int i = 0; i < columnDataB.Length; i++) { - return; + var column = new float[rowsY]; + GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column); + columnDataB[i] = column; } - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true); - } - - /// - /// Cache-Oblivious Matrix Multiplication - /// - /// if set to true transpose matrix A. - /// if set to true transpose matrix B. - /// The value to scale the matrix A with. - /// The matrix A. - /// Row-shift of the left matrix - /// Column-shift of the left matrix - /// The matrix B. - /// Row-shift of the right matrix - /// Column-shift of the right matrix - /// The matrix C. - /// Row-shift of the result matrix - /// Column-shift of the result matrix - /// The number of rows of matrix op(A) and of the matrix C. - /// The number of columns of matrix op(B) and of the matrix C. - /// The number of columns of matrix op(A) and the rows of the matrix op(B). - /// The constant number of rows of matrix op(A) and of the matrix C. - /// The constant number of columns of matrix op(B) and of the matrix C. - /// The constant number of columns of matrix op(A) and the rows of the matrix op(B). - /// Indicates if this is the first recursion. - static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, float alpha, float[] matrixA, int shiftArow, int shiftAcol, float[] matrixB, int shiftBrow, int shiftBcol, float[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first) - { - if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1) + var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2; + if (shouldNotParallelize) { - if ((int) transposeA > 111 && (int) transposeB > 111) + var row = new float[columnsX]; + for (int i = 0; i < rowsX; i++) { - for (var m1 = 0; m1 < m; m1++) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - float sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } - } - } - else if ((int) transposeA > 111) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + var col = columnDataB[j]; + float sum = 0; + for (int ii = 0; ii < row.Length; ii++) { - var matBcolPos = n1 + shiftBcol; - float sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; - } - - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + sum += row[ii] * col[ii]; } - } - } - else if ((int) transposeB > 111) - { - for (var m1 = 0; m1 < m; m1++) - { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) - { - var matBcolPos = n1 + shiftBcol; - float sum = 0; - for (var k1 = 0; k1 < k; ++k1) - { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[((k1 + shiftBrow)*constN) + matBcolPos]; - } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; - } + result[j * rowsX + i] += 1.0f * sum; } } - else + } + else + { + CommonParallel.For(0, rowsX, 1, (u, v) => { - for (var m1 = 0; m1 < m; m1++) + var row = new float[columnsX]; + for (int i = u; i < v; i++) { - var matArowPos = m1 + shiftArow; - var matCrowPos = m1 + shiftCrow; - for (var n1 = 0; n1 < n; ++n1) + GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row); + for (int j = 0; j < columnsY; j++) { - var matBcolPos = n1 + shiftBcol; + var column = columnDataB[j]; float sum = 0; - for (var k1 = 0; k1 < k; ++k1) + for (int ii = 0; ii < row.Length; ii++) { - sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]* - matrixB[(matBcolPos*constK) + k1 + shiftBrow]; + sum += row[ii] * column[ii]; } - result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum; + result[j * rowsX + i] += 1.0f * sum; } } - } - } - else - { - // divide and conquer - int m2 = m/2, n2 = n/2, k2 = k/2; - - if (first) - { - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false)); - - CommonParallel.Invoke( - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false), - () => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false)); - } - else - { - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false); - - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false); - CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false); - } + }); } } - 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) + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + 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 (a == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs index 2b912c88..681f957b 100644 --- a/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs @@ -42,18 +42,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed /// internal partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider { - private readonly Variation _variation; - - internal ManagedLinearAlgebraProvider() - { - _variation = Variation.Experimental; - } - - internal ManagedLinearAlgebraProvider(Variation variation) - { - _variation = variation; - } - /// /// 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. diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs index dac721f2..970d2a8d 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs @@ -468,12 +468,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { @@ -554,12 +548,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// The c matrix. 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). @@ -938,128 +926,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } } - 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; - } - } - }); - } - } - /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs index ed95b2a6..f283f812 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs @@ -470,12 +470,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { @@ -556,12 +550,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// The c matrix. 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). @@ -940,128 +928,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } } - 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; - } - } - }); - } - } - /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs index f5bba15e..c18a489c 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs @@ -465,12 +465,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { @@ -551,12 +545,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// The c matrix. 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). @@ -826,128 +814,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } } - 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; - } - } - }); - } - } - /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs index 7fe3c4fb..4cfffa48 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs @@ -465,12 +465,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// set to 1.0 and beta set to 0.0, and x and y are not transposed. 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) { @@ -551,12 +545,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// The c matrix. 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). @@ -826,128 +814,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference } } - 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; - } - } - }); - } - } - /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs index 0c18cd1c..f0b4952b 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs @@ -42,18 +42,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference /// internal partial class ManagedReferenceLinearAlgebraProvider : ILinearAlgebraProvider { - private readonly Variation _variation; - - internal ManagedReferenceLinearAlgebraProvider() - { - _variation = Variation.Experimental; - } - - internal ManagedReferenceLinearAlgebraProvider(Variation variation) - { - _variation = variation; - } - /// /// 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.