diff --git a/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs b/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs index b6d88a30..40aa8e23 100644 --- a/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs +++ b/src/Benchmark/LinearAlgebra/DenseMatrixProduct.cs @@ -14,6 +14,7 @@ namespace Benchmark.LinearAlgebra readonly Dictionary> _data = new Dictionary>(); 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 MathNetManagedExperimental() + { + Control.LinearAlgebraProvider = _managedExperimental; + return _data[Key(M, N)].TransposeAndMultiply(_data[Key(M, N)]); + } + [Benchmark(OperationsPerInvoke = 1)] public Matrix 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 /// /// Assumes that and have already been transposed. /// - 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; } } } diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 30d887d7..8ebc6ba8 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/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. 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 /// 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). @@ -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; + } + } + }); + } + } + /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index a0351631..9ddc97de 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/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. 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 /// 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). @@ -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; + } + } + }); + } + } + /// /// Computes the LUP factorization of A. P*A = L*U. /// diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 05a0da68..b8456706 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/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. 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 /// 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). @@ -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; + } + } + }); + } + } + /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -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++) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 5435eeda..a877f72b 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/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. 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 /// 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). @@ -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; + } + } + }); + } + } + /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -1778,7 +1912,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new float[rowsA]; - + const int maxiter = 1000; var e = new float[columnsA]; diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 6fbc98a8..7d4a942b 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -27,13 +27,33 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; + namespace MathNet.Numerics.Providers.LinearAlgebra { + internal enum Variation + { + Original, + Experimental + } + /// /// The managed linear algebra provider. /// public partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider { + private Variation _variation; + + public ManagedLinearAlgebraProvider() + { + _variation = Variation.Original; + } + + 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. @@ -54,5 +74,42 @@ namespace MathNet.Numerics.Providers.LinearAlgebra { return "Managed"; } + + + /// + /// Assumes that and have already been transposed. + /// + protected static void GetRow(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); + } + } + + /// + /// Assumes that and have already been transposed. + /// + protected static void GetColumn(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]; + } + } + } } }