From 08f1483e877d5f307c679855ab2a282715c0e674 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Sat, 25 Jun 2011 19:34:35 +0800 Subject: [PATCH] tweaked cholesky solve fixed issue with LU solve for complex types using parallel element copy instead of Array.Copy --- .../ManagedLinearAlgebraProvider.Complex.cs | 152 +++++++++-------- .../ManagedLinearAlgebraProvider.Complex32.cs | 152 +++++++++-------- .../ManagedLinearAlgebraProvider.Double.cs | 157 ++++++++++-------- .../ManagedLinearAlgebraProvider.Single.cs | 154 +++++++++-------- src/Numerics/Control.cs | 17 +- 5 files changed, 353 insertions(+), 279 deletions(-) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 6733931e..5b84df9d 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -67,35 +67,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (alpha.IsZero()) { - Array.Copy(y, result, y.Length); + y.Copy(result); } else if (alpha.IsOne()) { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + x[index]; } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); - } } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + (alpha * x[index]); } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); - } } } @@ -119,21 +119,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } else if (alpha.IsOne()) { - Array.Copy(x, result, x.Length); + x.Copy(result); } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = alpha * x[index]; } } - else - { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); - } } } @@ -202,17 +202,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] + y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); - } } /// @@ -247,17 +247,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] - y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); - } } /// @@ -292,17 +292,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] * y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); - } } /// @@ -337,17 +337,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] / y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); - } } /// @@ -1065,7 +1065,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } LUSolveFactored(order, a, order, ipiv, inverse); - CommonParallel.For(0, a.Length, index => a[index] = inverse[index]); + inverse.Copy(a); } /// @@ -1134,7 +1134,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var ipiv = new int[order]; var clone = new Complex[a.Length]; - Array.Copy(a, 0, clone, 0, a.Length); + a.Copy(clone); LUFactor(clone, order, ipiv); LUSolveFactored(columnsOfB, clone, order, ipiv, b); } @@ -1357,7 +1357,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new Complex[a.Length]; - Array.Copy(a, 0, clone, 0, a.Length); + a.Copy(clone); CholeskyFactor(clone, orderA); CholeskySolveFactored(clone, orderA, b, columnsB); } @@ -1392,39 +1392,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CommonParallel.For( - 0, - columnsB, - c => + if (Control.ParallelizeOperation(columnsB * 10)) + { + CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); + } + else + { + for (var index = 0; index < columnsB; index++) { - var cindex = c * orderA; + DoCholeskySolve(a, orderA, b, index); + } + } + } - // Solve L*Y = B; - Complex sum; - for (var i = 0; i < orderA; i++) - { - sum = b[cindex + i]; - for (var k = i - 1; k >= 0; k--) - { - sum -= a[(k * orderA) + i] * b[cindex + k]; - } + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. Has to be different than . + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The column to solve for. + private static void DoCholeskySolve(Complex[] a, int orderA, Complex[] b, int index) + { + var cindex = index * orderA; - b[cindex + i] = sum / a[(i * orderA) + i]; - } + // Solve L*Y = B; + Complex sum; + for (var i = 0; i < orderA; i++) + { + sum = b[cindex + i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= a[(k * orderA) + i] * b[cindex + k]; + } - // Solve L'*X = Y; - for (var i = orderA - 1; i >= 0; i--) - { - sum = b[cindex + i]; - var iindex = i * orderA; - for (var k = i + 1; k < orderA; k++) - { - sum -= a[iindex + k].Conjugate() * b[cindex + k]; - } + b[cindex + i] = sum / a[(i * orderA) + i]; + } - b[cindex + i] = sum / a[iindex + i]; - } - }); + // Solve L'*X = Y; + for (var i = orderA - 1; i >= 0; i--) + { + sum = b[cindex + i]; + var iindex = i * orderA; + for (var k = i + 1; k < orderA; k++) + { + sum -= a[iindex + k].Conjugate() * b[cindex + k]; + } + + b[cindex + i] = sum / a[iindex + i]; + } } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 23f28698..61c955e4 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (alpha.IsZero()) { - Array.Copy(y, result, y.Length); + y.Copy(result); } else if (alpha.IsOne()) { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + x[index]; } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); - } } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + (alpha * x[index]); } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); - } } } @@ -114,21 +114,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } else if (alpha.IsOne()) { - Array.Copy(x, result, x.Length); + x.Copy(result); } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = alpha * x[index]; } } - else - { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); - } } } @@ -198,17 +198,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] + y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); - } } /// @@ -243,17 +243,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] - y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); - } } /// @@ -288,17 +288,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] * y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); - } } /// @@ -333,17 +333,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] / y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); - } } /// @@ -1061,7 +1061,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } LUSolveFactored(order, a, order, ipiv, inverse); - CommonParallel.For(0, a.Length, index => a[index] = inverse[index]); + inverse.Copy(a); } /// @@ -1130,7 +1130,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var ipiv = new int[order]; var clone = new Complex32[a.Length]; - Array.Copy(a, 0, clone, 0, a.Length); + a.Copy(clone); LUFactor(clone, order, ipiv); LUSolveFactored(columnsOfB, clone, order, ipiv, b); } @@ -1353,7 +1353,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new Complex32[a.Length]; - Array.Copy(a, 0, clone, 0, a.Length); + a.Copy(clone); CholeskyFactor(clone, orderA); CholeskySolveFactored(clone, orderA, b, columnsB); } @@ -1388,39 +1388,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CommonParallel.For( - 0, - columnsB, - c => + if (Control.ParallelizeOperation(columnsB * 10)) + { + CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); + } + else + { + for (var index = 0; index < columnsB; index++) { - var cindex = c * orderA; + DoCholeskySolve(a, orderA, b, index); + } + } + } - // Solve L*Y = B; - Complex32 sum; - for (var i = 0; i < orderA; i++) - { - sum = b[cindex + i]; - for (var k = i - 1; k >= 0; k--) - { - sum -= a[(k * orderA) + i] * b[cindex + k]; - } + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. Has to be different than . + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The column to solve for. + private static void DoCholeskySolve(Complex32[] a, int orderA, Complex32[] b, int index) + { + var cindex = index * orderA; - b[cindex + i] = sum / a[(i * orderA) + i]; - } + // Solve L*Y = B; + Complex32 sum; + for (var i = 0; i < orderA; i++) + { + sum = b[cindex + i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= a[(k * orderA) + i] * b[cindex + k]; + } - // Solve L'*X = Y; - for (var i = orderA - 1; i >= 0; i--) - { - sum = b[cindex + i]; - var iindex = i * orderA; - for (var k = i + 1; k < orderA; k++) - { - sum -= a[iindex + k].Conjugate() * b[cindex + k]; - } + b[cindex + i] = sum / a[(i * orderA) + i]; + } - b[cindex + i] = sum / a[iindex + i]; - } - }); + // Solve L'*X = Y; + for (var i = orderA - 1; i >= 0; i--) + { + sum = b[cindex + i]; + var iindex = i * orderA; + for (var k = i + 1; k < orderA; k++) + { + sum -= a[iindex + k].Conjugate() * b[cindex + k]; + } + + b[cindex + i] = sum / a[iindex + i]; + } } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index f4b2be15..aabfc893 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (alpha == 0.0) { - Buffer.BlockCopy(y, 0, result, 0, y.Length * Constants.SizeOfDouble); + y.Copy(result); } else if (alpha == 1.0) { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + x[index]; } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); - } } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + (alpha * x[index]); } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); - } } } @@ -113,21 +113,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } else if (alpha == 1.0) { - Buffer.BlockCopy(x, 0, result, 0, x.Length * Constants.SizeOfDouble); + x.Copy(result); } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = alpha * x[index]; } } - else - { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); - } } } @@ -197,16 +197,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) { - for (var index = 0; index < x.Length; index++) - { - result[index] = x[index] + y[index]; - } + CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); } else { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + for (var index = 0; index < x.Length; index++) + { + result[index] = x[index] + y[index]; + } } } @@ -242,17 +242,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] - y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); - } } /// @@ -287,17 +287,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] * y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); - } } /// @@ -332,17 +332,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] / y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); - } } /// @@ -950,7 +950,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } LUSolveFactored(order, a, order, ipiv, inverse); - Buffer.BlockCopy(inverse, 0, a, 0, a.Length * Constants.SizeOfDouble); + inverse.Copy(a); } /// @@ -1019,7 +1019,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var ipiv = new int[order]; var clone = new double[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + a.Copy(clone); LUFactor(clone, order, ipiv); LUSolveFactored(columnsOfB, clone, order, ipiv, b); } @@ -1242,7 +1242,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new double[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + a.Copy(clone); CholeskyFactor(clone, orderA); CholeskySolveFactored(clone, orderA, b, columnsB); } @@ -1277,40 +1277,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - if( ) - CommonParallel.For( - 0, - columnsB, - c => + if (Control.ParallelizeOperation(columnsB * 10)) + { + CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); + } + else + { + for (var index = 0; index < columnsB; index++) { - var cindex = c * orderA; + DoCholeskySolve(a, orderA, b, index); + } + } + } - // Solve L*Y = B; - double sum; - for (var i = 0; i < orderA; i++) - { - sum = b[cindex + i]; - for (var k = i - 1; k >= 0; k--) - { - sum -= a[(k * orderA) + i] * b[cindex + k]; - } + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. Has to be different than . + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The column to solve for. + private static void DoCholeskySolve(double[] a, int orderA, double[] b, int index) + { + var cindex = index * orderA; - b[cindex + i] = sum / a[(i * orderA) + i]; - } + // Solve L*Y = B; + double sum; + for (var i = 0; i < orderA; i++) + { + sum = b[cindex + i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= a[(k * orderA) + i] * b[cindex + k]; + } - // Solve L'*X = Y; - for (var i = orderA - 1; i >= 0; i--) - { - sum = b[cindex + i]; - var iindex = i * orderA; - for (var k = i + 1; k < orderA; k++) - { - sum -= a[iindex + k] * b[cindex + k]; - } + b[cindex + i] = sum / a[(i * orderA) + i]; + } - b[cindex + i] = sum / a[iindex + i]; - } - }); + // Solve L'*X = Y; + for (var i = orderA - 1; i >= 0; i--) + { + sum = b[cindex + i]; + var iindex = i * orderA; + for (var k = i + 1; k < orderA; k++) + { + sum -= a[iindex + k] * b[cindex + k]; + } + + b[cindex + i] = sum / a[iindex + i]; + } } /// @@ -1642,7 +1657,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new double[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + a.Copy(clone); var q = new double[rows * rows]; QRFactor(clone, rows, columns, q, work); QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 165f0547..1a6b3162 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (alpha == 0.0) { - Buffer.BlockCopy(y, 0, result, 0, y.Length * Constants.SizeOfFloat); + y.Copy(result); } else if (alpha == 1.0) { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + x[index]; } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); - } } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = y[index] + (alpha * x[index]); } } - else - { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); - } } } @@ -113,21 +113,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } else if (alpha == 1.0) { - Buffer.BlockCopy(x, 0, result, 0, x.Length * Constants.SizeOfFloat); + x.Copy(result); } else { - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = alpha * x[index]; } } - else - { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); - } } } @@ -197,17 +197,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] + y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); - } } /// @@ -242,17 +242,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] - y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); - } } /// @@ -287,17 +287,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] * y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); - } } /// @@ -332,17 +332,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (x.Length < Control.ParallelizeElements || Control.DisableParallelization) + if (Control.ParallelizeOperation(x.Length)) + { + CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + } + else { for (var index = 0; index < x.Length; index++) { result[index] = x[index] / y[index]; } } - else - { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); - } } /// @@ -951,7 +951,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } LUSolveFactored(order, a, order, ipiv, inverse); - Buffer.BlockCopy(inverse, 0, a, 0, a.Length * Constants.SizeOfFloat); + inverse.Copy(a); } /// @@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var ipiv = new int[order]; var clone = new float[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + a.Copy(clone); LUFactor(clone, order, ipiv); LUSolveFactored(columnsOfB, clone, order, ipiv, b); } @@ -1243,7 +1243,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new float[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + a.Copy(clone); CholeskyFactor(clone, orderA); CholeskySolveFactored(clone, orderA, b, columnsB); } @@ -1278,39 +1278,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - CommonParallel.For( - 0, - columnsB, - c => + if (Control.ParallelizeOperation(columnsB * 10)) + { + CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); + } + else + { + for (var index = 0; index < columnsB; index++) { - var cindex = c * orderA; + DoCholeskySolve(a, orderA, b, index); + } + } + } - // Solve L*Y = B; - float sum; - for (var i = 0; i < orderA; i++) - { - sum = b[cindex + i]; - for (var k = i - 1; k >= 0; k--) - { - sum -= a[(k * orderA) + i] * b[cindex + k]; - } + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. Has to be different than . + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The column to solve for. + private static void DoCholeskySolve(float[] a, int orderA, float[] b, int index) + { + var cindex = index * orderA; - b[cindex + i] = sum / a[(i * orderA) + i]; - } + // Solve L*Y = B; + float sum; + for (var i = 0; i < orderA; i++) + { + sum = b[cindex + i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= a[(k * orderA) + i] * b[cindex + k]; + } - // Solve L'*X = Y; - for (var i = orderA - 1; i >= 0; i--) - { - sum = b[cindex + i]; - var iindex = i * orderA; - for (var k = i + 1; k < orderA; k++) - { - sum -= a[iindex + k] * b[cindex + k]; - } + b[cindex + i] = sum / a[(i * orderA) + i]; + } - b[cindex + i] = sum / a[iindex + i]; - } - }); + // Solve L'*X = Y; + for (var i = orderA - 1; i >= 0; i--) + { + sum = b[cindex + i]; + var iindex = i * orderA; + for (var k = i + 1; k < orderA; k++) + { + sum -= a[iindex + k] * b[cindex + k]; + } + + b[cindex + i] = sum / a[iindex + i]; + } } /// @@ -1642,7 +1658,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var clone = new float[a.Length]; - Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + a.Copy(clone); var q = new float[rows * rows]; QRFactor(clone, rows, columns, q, work); QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index b56c6c74..fbcca5a8 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -59,7 +59,7 @@ namespace MathNet.Numerics /// /// The default cutoff point for order size for the matrix multiply in linear algebra provider. /// - private static int _parallelizeElements = 256; + private static int _parallelizeElements = 300; /// /// Initializes static members of the Control class. @@ -99,7 +99,8 @@ namespace MathNet.Numerics /// The linear algebra provider. public static ILinearAlgebraProvider LinearAlgebraProvider { - get; set; + get; + set; } /// @@ -169,7 +170,7 @@ namespace MathNet.Numerics /// /// Gets or sets the number of elements a vector or matrix must contain before we multiply threads. /// - /// Number of elements. Default is 256. + /// Number of elements. Default is 300. public static int ParallelizeElements { get @@ -185,5 +186,15 @@ namespace MathNet.Numerics } } } + + /// + /// Given the number elements, should the operation be parallelized. + /// + /// The number elements to check. + /// true if the operation should be parallelized; false otherwise. + public static bool ParallelizeOperation(int elements) + { + return elements < ParallelizeElements || DisableParallelization || NumberOfParallelWorkerThreads < 2; + } } }