From d3dc54b4c128a44fc2fa3b16d0298fd7ba84256d Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 4 Apr 2013 19:52:09 +0200 Subject: [PATCH] Reworking common parallelization routines #92 --- .../ManagedLinearAlgebraProvider.Complex.cs | 182 +++++++++----- .../ManagedLinearAlgebraProvider.Complex32.cs | 178 ++++++++----- .../ManagedLinearAlgebraProvider.Double.cs | 180 ++++++++----- .../ManagedLinearAlgebraProvider.Single.cs | 178 ++++++++----- .../DiscreteFourierTransform.Naive.cs | 34 +-- .../DiscreteFourierTransform.RadixN.cs | 26 +- .../DiscreteHartleyTransform.Naive.cs | 28 +-- .../LinearAlgebra/Complex/DenseVector.cs | 45 ++-- .../Complex/Factorization/DenseCholesky.cs | 4 +- .../Complex/Factorization/DenseLU.cs | 4 +- .../LinearAlgebra/Complex/SparseMatrix.cs | 2 +- .../LinearAlgebra/Complex/SparseVector.cs | 10 +- .../LinearAlgebra/Complex32/DenseVector.cs | 45 ++-- .../Complex32/Factorization/DenseCholesky.cs | 4 +- .../Complex32/Factorization/DenseLU.cs | 4 +- .../LinearAlgebra/Complex32/SparseMatrix.cs | 2 +- .../LinearAlgebra/Complex32/SparseVector.cs | 10 +- .../LinearAlgebra/Double/DenseMatrix.cs | 26 +- .../LinearAlgebra/Double/DenseVector.cs | 34 +-- .../LinearAlgebra/Double/SparseMatrix.cs | 2 +- .../LinearAlgebra/Single/DenseMatrix.cs | 24 +- .../LinearAlgebra/Single/DenseVector.cs | 34 +-- .../LinearAlgebra/Single/SparseMatrix.cs | 2 +- src/Numerics/Random/Palf.cs | 78 +++--- src/Numerics/Threading/CommonParallel.cs | 236 ++++++++++-------- 25 files changed, 842 insertions(+), 530 deletions(-) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 9afc6360..26e4c909 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -76,7 +76,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + x[i]; + } + }); } else { @@ -90,7 +96,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + (alpha*x[i]); + } + }); } else { @@ -128,7 +140,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + CommonParallel.For(0, x.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = alpha*x[i]; + } + }); } else { @@ -207,7 +225,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] + y[i]; + } + }); } else { @@ -252,7 +276,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] - y[i]; + } + }); } else { @@ -297,7 +327,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]*y[i]; + } + }); } else { @@ -342,7 +378,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]/y[i]; + } + }); } else { @@ -1395,17 +1437,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - if (Control.ParallelizeOperation(columnsB * 10)) - { - CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); - } - else - { - for (var index = 0; index < columnsB; index++) + CommonParallel.For(0, columnsB, (u, v) => { - DoCholeskySolve(a, orderA, b, index); - } - } + for (int i = u; i < v; i++) + { + DoCholeskySolve(a, orderA, b, i); + } + }); } /// @@ -1553,7 +1591,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - CommonParallel.For(0, rowsR, i => q[(i * rowsR) + i] = Complex.One); + CommonParallel.For(0, rowsR, (a, b) => + { + for (int i = a; i < b; i++) + { + q[(i*rowsR) + i] = Complex.One; + } + }); var minmn = Math.Min(rowsR, columnsR); for (var i = 0; i < minmn; i++) @@ -1762,14 +1806,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var tmp = column * rowCount; var index = tmp + row; - CommonParallel.For( - row, - rowCount, - i => + CommonParallel.For(row, rowCount, (u, v) => { - var iIndex = tmp + i; - work[iIndex - row] = a[iIndex]; - a[iIndex] = Complex.Zero; + for (int i = u; i < v; i++) + { + var iIndex = tmp + i; + work[iIndex - row] = a[iIndex]; + a[iIndex] = Complex.Zero; + } }); var norm = Complex.Zero; @@ -1793,11 +1837,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } a[index] = -norm; - CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] /= norm; + } + }); work[tmp] += 1.0; var s = (1.0 / work[tmp]).SquareRoot(); - CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] = work[tmp + i].Conjugate()*s; + } + }); } #endregion @@ -1965,7 +2021,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } int rowsQ, columnsQ, rowsR, columnsR; - if( method == QRMethod.Full) + if (method == QRMethod.Full) { rowsQ = columnsQ = rowsR = rowsA; columnsR = columnsA; @@ -1976,24 +2032,24 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra columnsQ = rowsR = columnsR = columnsA; } - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR*columnsR), "r"); } - if (q.Length != rowsQ * columnsQ) + if (q.Length != rowsQ*columnsQ) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ*columnsQ), "q"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA*columnsB), "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA*columnsB), "x"); } var sol = new Complex[b.Length]; @@ -2005,56 +2061,62 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var column = new Complex[rowsA]; for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - CommonParallel.For(0, rowsA, k => column[k] = sol[jm + k]); - CommonParallel.For( - 0, - columnsA, - i => + var jm = j*rowsA; + CommonParallel.For(0, rowsA, (u, v) => { - var im = i * rowsA; - - var sum = Complex.Zero; - for (var k = 0; k < rowsA; k++) + for (int k = u; k < v; k++) { - sum += q[im + k].Conjugate() * column[k]; + column[k] = sol[jm + k]; } + }); + CommonParallel.For(0, columnsA, (u, v) => + { + for (int i = u; i < v; i++) + { + var im = i*rowsA; - sol[jm + i] = sum; + var sum = Complex.Zero; + for (var k = 0; k < rowsA; k++) + { + sum += q[im + k].Conjugate()*column[k]; + } + + sol[jm + i] = sum; + } }); } // Solve R*X = Y; for (var k = columnsA - 1; k >= 0; k--) { - var km = k * rowsR; + var km = k*rowsR; for (var j = 0; j < columnsB; j++) { - sol[(j * rowsA) + k] /= r[km + k]; + sol[(j*rowsA) + k] /= r[km + k]; } for (var i = 0; i < k; i++) { for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - sol[jm + i] -= sol[jm + k] * r[km + i]; + var jm = j*rowsA; + sol[jm + i] -= sol[jm + k]*r[km + i]; } } } // Fill result matrix - CommonParallel.For( - 0, - columnsR, - row => + CommonParallel.For(0, columnsR, (u, v) => { - for (var col = 0; col < columnsB; col++) + for (int row = u; row < v; row++) { - x[(col * columnsA) + row] = sol[row + (col * rowsA)]; + for (var col = 0; col < columnsB; col++) + { + x[(col*columnsA) + row] = sol[row + (col*rowsA)]; + } } }); - } + } /// /// Computes the singular value decomposition of A. @@ -2787,7 +2849,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // Copy stemp to s with size adjustment. We are using ported copy of linpack's svd code and it uses // a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed. // We should port lapack's svd routine to remove this problem. - CommonParallel.For(0, Math.Min(rowsA, columnsA), index => s[index] = stemp[index]); + Array.Copy(stemp, s, Math.Min(rowsA, columnsA)); // On return the first element of the work array stores the min size of the work array could have been // work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns)); diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 5a7a0216..c420ca1a 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -71,7 +71,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + x[i]; + } + }); } else { @@ -85,7 +91,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + (alpha*x[i]); + } + }); } else { @@ -124,7 +136,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + CommonParallel.For(0, x.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = alpha*x[i]; + } + }); } else { @@ -204,7 +222,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] + y[i]; + } + }); } else { @@ -249,7 +273,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] - y[i]; + } + }); } else { @@ -294,7 +324,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]*y[i]; + } + }); } else { @@ -339,7 +375,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]/y[i]; + } + }); } else { @@ -1392,17 +1434,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - if (Control.ParallelizeOperation(columnsB * 10)) - { - CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); - } - else - { - for (var index = 0; index < columnsB; index++) + CommonParallel.For(0, columnsB, (u, v) => { - DoCholeskySolve(a, orderA, b, index); - } - } + for (int i = u; i < v; i++) + { + DoCholeskySolve(a, orderA, b, i); + } + }); } /// @@ -1550,7 +1588,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - CommonParallel.For(0, rowsR, i => q[(i * rowsR) + i] = Complex32.One); + CommonParallel.For(0, rowsR, (a, b) => + { + for (int i = a; i < b; i++) + { + q[(i*rowsR) + i] = Complex32.One; + } + }); var minmn = Math.Min(rowsR, columnsR); for (var i = 0; i < minmn; i++) @@ -1759,14 +1803,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var tmp = column * rowCount; var index = tmp + row; - CommonParallel.For( - row, - rowCount, - i => + CommonParallel.For(row, rowCount, (u, v) => { - var iIndex = tmp + i; - work[iIndex - row] = a[iIndex]; - a[iIndex] = Complex32.Zero; + for (int i = u; i < v; i++) + { + var iIndex = tmp + i; + work[iIndex - row] = a[iIndex]; + a[iIndex] = Complex32.Zero; + } }); var norm = Complex32.Zero; @@ -1790,11 +1834,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } a[index] = -norm; - CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] /= norm; + } + }); work[tmp] += 1.0f; var s = (1.0f / work[tmp]).SquareRoot(); - CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] = work[tmp + i].Conjugate()*s; + } + }); } #endregion @@ -1973,24 +2029,24 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra columnsQ = rowsR = columnsR = columnsA; } - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR*columnsR), "r"); } - if (q.Length != rowsQ * columnsQ) + if (q.Length != rowsQ*columnsQ) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ*columnsQ), "q"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA*columnsB), "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA*columnsB), "x"); } var sol = new Complex32[b.Length]; @@ -2002,53 +2058,59 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var column = new Complex32[rowsA]; for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - CommonParallel.For(0, rowsA, k => column[k] = sol[jm + k]); - CommonParallel.For( - 0, - columnsA, - i => + var jm = j*rowsA; + CommonParallel.For(0, rowsA, (u, v) => { - var im = i * rowsA; - - var sum = Complex32.Zero; - for (var k = 0; k < rowsA; k++) + for (int k = u; k < v; k++) { - sum += q[im + k].Conjugate() * column[k]; + column[k] = sol[jm + k]; } + }); + CommonParallel.For(0, columnsA, (u, v) => + { + for (int i = u; i < v; i++) + { + var im = i*rowsA; - sol[jm + i] = sum; + var sum = Complex32.Zero; + for (var k = 0; k < rowsA; k++) + { + sum += q[im + k].Conjugate()*column[k]; + } + + sol[jm + i] = sum; + } }); } // Solve R*X = Y; for (var k = columnsA - 1; k >= 0; k--) { - var km = k * rowsR; + var km = k*rowsR; for (var j = 0; j < columnsB; j++) { - sol[(j * rowsA) + k] /= r[km + k]; + sol[(j*rowsA) + k] /= r[km + k]; } for (var i = 0; i < k; i++) { for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - sol[jm + i] -= sol[jm + k] * r[km + i]; + var jm = j*rowsA; + sol[jm + i] -= sol[jm + k]*r[km + i]; } } } // Fill result matrix - CommonParallel.For( - 0, - columnsR, - row => + CommonParallel.For(0, columnsR, (u, v) => { - for (var col = 0; col < columnsB; col++) + for (int row = u; row < v; row++) { - x[(col * columnsA) + row] = sol[row + (col * rowsA)]; + for (var col = 0; col < columnsB; col++) + { + x[(col*columnsA) + row] = sol[row + (col*rowsA)]; + } } }); } @@ -2784,7 +2846,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // Copy stemp to s with size adjustment. We are using ported copy of linpack's svd code and it uses // a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed. // We should port lapack's svd routine to remove this problem. - CommonParallel.For(0, Math.Min(rowsA, columnsA), index => s[index] = stemp[index]); + Array.Copy(stemp, s, Math.Min(rowsA, columnsA)); // On return the first element of the work array stores the min size of the work array could have been // work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns)); diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 615010f9..a475ce47 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -70,7 +70,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + x[i]; + } + }); } else { @@ -84,7 +90,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + (alpha*x[i]); + } + }); } else { @@ -122,7 +134,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + CommonParallel.For(0, x.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = alpha*x[i]; + } + }); } else { @@ -202,7 +220,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] + y[i]; + } + }); } else { @@ -247,7 +271,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] - y[i]; + } + }); } else { @@ -292,7 +322,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]*y[i]; + } + }); } else { @@ -337,7 +373,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]/y[i]; + } + }); } else { @@ -1280,17 +1322,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - if (Control.ParallelizeOperation(columnsB * 10)) - { - CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); - } - else - { - for (var index = 0; index < columnsB; index++) + CommonParallel.For(0, columnsB, (u, v) => { - DoCholeskySolve(a, orderA, b, index); - } - } + for (int i = u; i < v; i++) + { + DoCholeskySolve(a, orderA, b, i); + } + }); } /// @@ -1439,7 +1477,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - CommonParallel.For(0, rowsR, i => q[(i * rowsR) + i] = 1.0); + CommonParallel.For(0, rowsR, (a, b) => + { + for (int i = a; i < b; i++) + { + q[(i*rowsR) + i] = 1.0; + } + }); var minmn = Math.Min(rowsR, columnsR); for (var i = 0; i < minmn; i++) @@ -1647,14 +1691,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var tmp = column * rowCount; var index = tmp + row; - CommonParallel.For( - row, - rowCount, - i => + CommonParallel.For(row, rowCount, (u, v) => { - var iIndex = tmp + i; - work[iIndex - row] = a[iIndex]; - a[iIndex] = 0.0; + for (int i = u; i < v; i++) + { + var iIndex = tmp + i; + work[iIndex - row] = a[iIndex]; + a[iIndex] = 0.0; + } }); var norm = 0.0; @@ -1679,11 +1723,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } a[index] = -1.0 / scale; - CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] *= scale; + } + }); work[tmp] += 1.0; var s = Math.Sqrt(1.0 / work[tmp]); - CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] *= s; + } + }); } #endregion @@ -1850,7 +1906,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } int rowsQ, columnsQ, rowsR, columnsR; - if( method == QRMethod.Full) + if (method == QRMethod.Full) { rowsQ = columnsQ = rowsR = rowsA; columnsR = columnsA; @@ -1861,82 +1917,88 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra columnsQ = rowsR = columnsR = columnsA; } - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR*columnsR), "r"); } - if (q.Length != rowsQ * columnsQ) + if (q.Length != rowsQ*columnsQ) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ*columnsQ), "q"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA*columnsB), "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA*columnsB), "x"); } var sol = new double[b.Length]; // Copy B matrix to "sol", so B data will not be changed - Buffer.BlockCopy(b, 0, sol, 0, b.Length * Constants.SizeOfDouble); + Buffer.BlockCopy(b, 0, sol, 0, b.Length*Constants.SizeOfDouble); // Compute Y = transpose(Q)*B var column = new double[rowsA]; for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - CommonParallel.For(0, rowsA, k => column[k] = sol[jm + k]); - CommonParallel.For( - 0, - columnsA, - i => + var jm = j*rowsA; + CommonParallel.For(0, rowsA, (u, v) => { - var im = i * rowsA; - - var sum = 0.0; - for (var k = 0; k < rowsA; k++) + for (int k = u; k < v; k++) { - sum += q[im + k] * column[k]; + column[k] = sol[jm + k]; } + }); + CommonParallel.For(0, columnsA, (u, v) => + { + for (int i = u; i < v; i++) + { + var im = i*rowsA; + + var sum = 0.0; + for (var k = 0; k < rowsA; k++) + { + sum += q[im + k]*column[k]; + } - sol[jm + i] = sum; + sol[jm + i] = sum; + } }); } // Solve R*X = Y; for (var k = columnsA - 1; k >= 0; k--) { - var km = k * rowsR; + var km = k*rowsR; for (var j = 0; j < columnsB; j++) { - sol[(j * rowsA) + k] /= r[km + k]; + sol[(j*rowsA) + k] /= r[km + k]; } for (var i = 0; i < k; i++) { for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - sol[jm + i] -= sol[jm + k] * r[km + i]; + var jm = j*rowsA; + sol[jm + i] -= sol[jm + k]*r[km + i]; } } } // Fill result matrix - CommonParallel.For( - 0, - columnsR, - row => + CommonParallel.For(0, columnsR, (u, v) => { - for (var col = 0; col < columnsB; col++) + for (int row = u; row < v; row++) { - x[(col * columnsA) + row] = sol[row + (col * rowsA)]; + for (var col = 0; col < columnsB; col++) + { + x[(col*columnsA) + row] = sol[row + (col*rowsA)]; + } } }); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index ad3cb9b9..356f5561 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -70,7 +70,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + x[i]; + } + }); } else { @@ -84,7 +90,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index])); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = y[i] + (alpha*x[i]); + } + }); } else { @@ -122,7 +134,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, x.Length, index => { result[index] = alpha * x[index]; }); + CommonParallel.For(0, x.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = alpha*x[i]; + } + }); } else { @@ -202,7 +220,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] + y[i]; + } + }); } else { @@ -247,7 +271,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i] - y[i]; + } + }); } else { @@ -292,7 +322,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]*y[i]; + } + }); } else { @@ -337,7 +373,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (Control.ParallelizeOperation(x.Length)) { - CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; }); + CommonParallel.For(0, y.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + result[i] = x[i]/y[i]; + } + }); } else { @@ -1281,17 +1323,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentReferenceDifferent); } - if (Control.ParallelizeOperation(columnsB * 10)) - { - CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c)); - } - else - { - for (var index = 0; index < columnsB; index++) + CommonParallel.For(0, columnsB, (u, v) => { - DoCholeskySolve(a, orderA, b, index); - } - } + for (int i = u; i < v; i++) + { + DoCholeskySolve(a, orderA, b, i); + } + }); } /// @@ -1439,7 +1477,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - CommonParallel.For(0, rowsR, i => q[(i * rowsR) + i] = 1.0f); + CommonParallel.For(0, rowsR, (a, b) => + { + for (int i = a; i < b; i++) + { + q[(i*rowsR) + i] = 1.0f; + } + }); var minmn = Math.Min(rowsR, columnsR); for (var i = 0; i < minmn; i++) @@ -1648,14 +1692,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var tmp = column * rowCount; var index = tmp + row; - CommonParallel.For( - row, - rowCount, - i => + CommonParallel.For(row, rowCount, (u, v) => { - var iIndex = tmp + i; - work[iIndex - row] = a[iIndex]; - a[iIndex] = 0.0f; + for (int i = u; i < v; i++) + { + var iIndex = tmp + i; + work[iIndex - row] = a[iIndex]; + a[iIndex] = 0.0f; + } }); var norm = 0.0; @@ -1680,11 +1724,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } a[index] = -1.0f / scale; - CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] *= scale; + } + }); work[tmp] += 1.0f; var s = (float)Math.Sqrt(1.0 / work[tmp]); - CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s); + CommonParallel.For(0, rowCount - row, (u, v) => + { + for (int i = u; i < v; i++) + { + work[tmp + i] *= s; + } + }); } #endregion @@ -1863,82 +1919,88 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra columnsQ = rowsR = columnsR = columnsA; } - if (r.Length != rowsR * columnsR) + if (r.Length != rowsR*columnsR) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR*columnsR), "r"); } - if (q.Length != rowsQ * columnsQ) + if (q.Length != rowsQ*columnsQ) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ*columnsQ), "q"); } - if (b.Length != rowsA * columnsB) + if (b.Length != rowsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA*columnsB), "b"); } - if (x.Length != columnsA * columnsB) + if (x.Length != columnsA*columnsB) { - throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA*columnsB), "x"); } var sol = new float[b.Length]; // Copy B matrix to "sol", so B data will not be changed - Buffer.BlockCopy(b, 0, sol, 0, b.Length * Constants.SizeOfFloat); + Buffer.BlockCopy(b, 0, sol, 0, b.Length*Constants.SizeOfFloat); // Compute Y = transpose(Q)*B var column = new float[rowsA]; for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - CommonParallel.For(0, rowsA, k => column[k] = sol[jm + k]); - CommonParallel.For( - 0, - columnsA, - i => + var jm = j*rowsA; + CommonParallel.For(0, rowsA, (u, v) => { - var im = i * rowsA; - - var sum = 0.0f; - for (var k = 0; k < rowsA; k++) + for (int k = u; k < v; k++) { - sum += q[im + k] * column[k]; + column[k] = sol[jm + k]; } + }); + CommonParallel.For(0, columnsA, (u, v) => + { + for (int i = u; i < v; i++) + { + var im = i*rowsA; + + var sum = 0.0f; + for (var k = 0; k < rowsA; k++) + { + sum += q[im + k]*column[k]; + } - sol[jm + i] = sum; + sol[jm + i] = sum; + } }); } // Solve R*X = Y; for (var k = columnsA - 1; k >= 0; k--) { - var km = k * rowsR; + var km = k*rowsR; for (var j = 0; j < columnsB; j++) { - sol[(j * rowsA) + k] /= r[km + k]; + sol[(j*rowsA) + k] /= r[km + k]; } for (var i = 0; i < k; i++) { for (var j = 0; j < columnsB; j++) { - var jm = j * rowsA; - sol[jm + i] -= sol[jm + k] * r[km + i]; + var jm = j*rowsA; + sol[jm + i] -= sol[jm + k]*r[km + i]; } } } // Fill result matrix - CommonParallel.For( - 0, - columnsR, - row => + CommonParallel.For(0, columnsR, (u, v) => { - for (var col = 0; col < columnsB; col++) + for (int row = u; row < v; row++) { - x[(col * columnsA) + row] = sol[row + (col * rowsA)]; + for (var col = 0; col < columnsB; col++) + { + x[(col*columnsA) + row] = sol[row + (col*rowsA)]; + } } }); } diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs index e84bacc9..3db327da 100644 --- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs +++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.Naive.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -47,24 +47,24 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms /// Corresponding frequency-space vector. internal static Complex[] Naive(Complex[] samples, int exponentSign) { - var w0 = exponentSign * Constants.Pi2 / samples.Length; + var w0 = exponentSign*Constants.Pi2/samples.Length; var spectrum = new Complex[samples.Length]; - CommonParallel.For( - 0, - samples.Length, - index => - { - var wk = w0 * index; - var sum = Complex.Zero; - for (var n = 0; n < samples.Length; n++) - { - var w = n * wk; - sum += samples[n] * new Complex(Math.Cos(w), Math.Sin(w)); - } + CommonParallel.For(0, samples.Length, (u, v) => + { + for (int i = u; i < v; i++) + { + var wk = w0*i; + var sum = Complex.Zero; + for (var n = 0; n < samples.Length; n++) + { + var w = n*wk; + sum += samples[n]*new Complex(Math.Cos(w), Math.Sin(w)); + } - spectrum[index] = sum; - }); + spectrum[i] = sum; + } + }); return spectrum; } @@ -95,4 +95,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms return timeSpace; } } -} \ No newline at end of file +} diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs index 4a6b396b..044c9588 100644 --- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs +++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteFourierTransform.RadixN.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -46,7 +46,7 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms /// /// Sample type /// Sample vector - private static void Radix2Reorder(T[] samples) + static void Radix2Reorder(T[] samples) { var j = 0; for (var i = 0; i < samples.Length - 1; i++) @@ -64,8 +64,7 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms { m >>= 1; j ^= m; - } - while ((j & m) == 0); + } while ((j & m) == 0); } } @@ -76,17 +75,17 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms /// Fourier series exponent sign. /// Level Group Size. /// Index inside of the level. - private static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) + static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) { // Twiddle Factor - var exponent = (exponentSign * k) * Constants.Pi / levelSize; + var exponent = (exponentSign*k)*Constants.Pi/levelSize; var w = new Complex(Math.Cos(exponent), Math.Sin(exponent)); var step = levelSize << 1; for (var i = k; i < samples.Length; i += step) { var ai = samples[i]; - var t = w * samples[i + levelSize]; + var t = w*samples[i + levelSize]; samples[i] = ai + t; samples[i + levelSize] = ai - t; } @@ -133,10 +132,13 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms { var size = levelSize; - CommonParallel.For( - 0, - size, - index => Radix2Step(samples, exponentSign, size, index)); + CommonParallel.For(0, size, (u, v) => + { + for (int i = u; i < v; i++) + { + Radix2Step(samples, exponentSign, size, i); + } + }); } } @@ -164,4 +166,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms InverseScaleByOptions(options, samples); } } -} \ No newline at end of file +} diff --git a/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs b/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs index 19dd487c..6c5905d9 100644 --- a/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs +++ b/src/Numerics/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -45,23 +45,23 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms /// Corresponding frequency-space vector. internal static double[] Naive(double[] samples) { - var w0 = Constants.Pi2 / samples.Length; + var w0 = Constants.Pi2/samples.Length; var spectrum = new double[samples.Length]; - CommonParallel.For( - 0, - samples.Length, - index => + CommonParallel.For(0, samples.Length, (u, v) => { - var wk = w0 * index; - var sum = 0.0; - for (var n = 0; n < samples.Length; n++) + for (int i = u; i < v; i++) { - var w = n * wk; - sum += samples[n] * Constants.Sqrt2 * Math.Cos(w - Constants.PiOver4); - } + var wk = w0*i; + var sum = 0.0; + for (var n = 0; n < samples.Length; n++) + { + var w = n*wk; + sum += samples[n]*Constants.Sqrt2*Math.Cos(w - Constants.PiOver4); + } - spectrum[index] = sum; + spectrum[i] = sum; + } }); return spectrum; @@ -93,4 +93,4 @@ namespace MathNet.Numerics.IntegralTransforms.Algorithms return timeSpace; } } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs index c358d07d..2501da9c 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs @@ -233,10 +233,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] + scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] + scalar; + } + }); } } @@ -292,10 +295,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] - scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] - scalar; + } + }); } } @@ -625,14 +631,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } var matrix = new DenseMatrix(u.Count, v.Count); - CommonParallel.For( - 0, - u.Count, - i => + CommonParallel.For(0, u.Count, (a, b) => { - for (var j = 0; j < v.Count; j++) + for (int i = a; i < b; i++) { - matrix.At(i, j, u._values[i] * v._values[j]); + for (var j = 0; j < v.Count; j++) + { + matrix.At(i, j, u._values[i]*v._values[j]); + } } }); return matrix; @@ -860,10 +866,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } - CommonParallel.For( - 0, - _length, - index => resultDense._values[index] = _values[index].Conjugate()); + CommonParallel.For(0, _length, (a, b) => + { + for (int i = a; i < b; i++) + { + resultDense._values[i] = _values[i].Conjugate(); + } + }); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs index fc32af9f..9d09c5c7 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs @@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; @@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs index 76307419..9bde144e 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs @@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; @@ -170,7 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index ea851311..602f44f9 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -805,7 +805,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._storage.Values[index] *= scalar); + Control.LinearAlgebraProvider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values); } } diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs index 6accc36e..0ecb330b 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs @@ -202,8 +202,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex if (_storage.ValueCount != 0) { - CommonParallel.For(0, _storage.ValueCount, index => targetSparse._storage.Values[index] = _storage.Values[index].Conjugate()); - Buffer.BlockCopy(_storage.Indices, 0, targetSparse._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt); + CommonParallel.For(0, _storage.ValueCount, (a, b) => + { + for (int i = a; i < b; i++) + { + targetSparse._storage.Values[i] = _storage.Values[i].Conjugate(); + } + }); + Buffer.BlockCopy(_storage.Indices, 0, targetSparse._storage.Indices, 0, _storage.ValueCount*Constants.SizeOfInt); } } diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs index 61a9d978..2ea199f4 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs @@ -232,10 +232,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] + scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] + scalar; + } + }); } } @@ -291,10 +294,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] - scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] - scalar; + } + }); } } @@ -624,14 +630,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } var matrix = new DenseMatrix(u.Count, v.Count); - CommonParallel.For( - 0, - u.Count, - i => + CommonParallel.For(0, u.Count, (a, b) => { - for (var j = 0; j < v.Count; j++) + for (int i = a; i < b; i++) { - matrix.At(i, j, u._values[i] * v._values[j]); + for (var j = 0; j < v.Count; j++) + { + matrix.At(i, j, u._values[i]*v._values[j]); + } } }); return matrix; @@ -859,10 +865,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } - CommonParallel.For( - 0, - _length, - index => resultDense._values[index] = _values[index].Conjugate()); + CommonParallel.For(0, _length, (a, b) => + { + for (int i = a; i < b; i++) + { + resultDense._values[i] = _values[i].Conjugate(); + } + }); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs index e7403170..6842128f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs @@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; @@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. var dfactor = (DenseMatrix)CholeskyFactor; diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs index 5f805c85..efa1b6b8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs @@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; @@ -170,7 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization } // Copy the contents of input to result. - CommonParallel.For(0, dinput.Values.Length, index => dresult.Values[index] = dinput.Values[index]); + Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index 472c3ea4..e4b68103 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -804,7 +804,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._storage.Values[index] *= scalar); + Control.LinearAlgebraProvider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values); } } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs index cd84d739..a3f731f5 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs @@ -202,8 +202,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 if (_storage.ValueCount != 0) { - CommonParallel.For(0, _storage.ValueCount, index => targetSparse._storage.Values[index] = _storage.Values[index].Conjugate()); - Buffer.BlockCopy(_storage.Indices, 0, targetSparse._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt); + CommonParallel.For(0, _storage.ValueCount, (a, b) => + { + for (int i = a; i < b; i++) + { + targetSparse._storage.Values[i] = _storage.Values[i].Conjugate(); + } + }); + Buffer.BlockCopy(_storage.Indices, 0, targetSparse._storage.Indices, 0, _storage.ValueCount*Constants.SizeOfInt); } } diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index cce7ea7c..17da0bd7 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -559,23 +559,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double protected override void DoModulus(double divisor, Matrix result) { var denseResult = result as DenseMatrix; - if (denseResult == null) { - base.DoModulus(divisor, result); + base.DoModulus(divisor, result); + return; } - else - { - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - CommonParallel.For( - 0, - _values.Length, - index => denseResult._values[index] %= divisor); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); } + + CommonParallel.For(0, _values.Length, (a, b) => + { + var v = denseResult._values; + for (int i = a; i < b; i++) + { + v[i] %= divisor; + } + }); } /// diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs index 64aafe68..fb719ec7 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -233,10 +233,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] + scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] + scalar; + } + }); } } @@ -302,10 +305,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] - scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] - scalar; + } + }); } } @@ -716,14 +722,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var matrix = new DenseMatrix(u.Count, v.Count); - CommonParallel.For( - 0, - u.Count, - i => + CommonParallel.For(0, u.Count, (a, b) => { - for (var j = 0; j < v.Count; j++) + for (int i = a; i < b; i++) { - matrix.At(i, j, u._values[i] * v._values[j]); + for (var j = 0; j < v.Count; j++) + { + matrix.At(i, j, u._values[i]*v._values[j]); + } } }); return matrix; diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 3487214c..938a6c12 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -803,7 +803,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CopyTo(sparseResult); } - CommonParallel.For(0, _storage.ValueCount, index => sparseResult._storage.Values[index] *= scalar); + Control.LinearAlgebraProvider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values); } } diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index 5ddd280d..a4c84e08 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -559,23 +559,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single protected override void DoModulus(float divisor, Matrix result) { var denseResult = result as DenseMatrix; - if (denseResult == null) { base.DoModulus(divisor, result); + return; } - else - { - if (!ReferenceEquals(this, result)) - { - CopyTo(result); - } - CommonParallel.For( - 0, - _values.Length, - index => denseResult._values[index] %= divisor); + if (!ReferenceEquals(this, result)) + { + CopyTo(result); } + + CommonParallel.For(0, _values.Length, (a, b) => + { + var v = denseResult._values; + for (int i = a; i < b; i++) + { + v[i] %= divisor; + } + }); } /// diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs index 745dcf5d..55c0dcad 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs @@ -232,10 +232,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] + scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] + scalar; + } + }); } } @@ -291,10 +294,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single } else { - CommonParallel.For( - 0, - _values.Length, - index => dense._values[index] = _values[index] - scalar); + CommonParallel.For(0, _values.Length, (a, b) => + { + for (int i = a; i < b; i++) + { + dense._values[i] = _values[i] - scalar; + } + }); } } @@ -705,14 +711,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single } var matrix = new DenseMatrix(u.Count, v.Count); - CommonParallel.For( - 0, - u.Count, - i => + CommonParallel.For(0, u.Count, (a, b) => { - for (var j = 0; j < v.Count; j++) + for (int i = a; i < b; i++) { - matrix.At(i, j, u._values[i] * v._values[j]); + for (var j = 0; j < v.Count; j++) + { + matrix.At(i, j, u._values[i]*v._values[j]); + } } }); return matrix; diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index 292c2dd2..363d7c5b 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -802,7 +802,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single CopyTo(sparseResult); } - CommonParallel.For(0, NonZerosCount, index => sparseResult._storage.Values[index] *= scalar); + Control.LinearAlgebraProvider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values); } } diff --git a/src/Numerics/Random/Palf.cs b/src/Numerics/Random/Palf.cs index 4ebc7282..115c0066 100644 --- a/src/Numerics/Random/Palf.cs +++ b/src/Numerics/Random/Palf.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -48,17 +48,17 @@ namespace MathNet.Numerics.Random /// /// Default value for the ShortLag /// - private const int DefaultShortLag = 418; + const int DefaultShortLag = 418; /// /// Default value for the LongLag /// - private const int DefaultLongLag = 1279; + const int DefaultLongLag = 1279; /// /// The multiplier to compute a double-precision floating point number [0, 1) /// - private const double IntToDoubleMultiplier = 1.0 / (int.MaxValue + 1.0); + const double IntToDoubleMultiplier = 1.0/(int.MaxValue + 1.0); /// /// Initializes a new instance of the class using @@ -67,7 +67,8 @@ namespace MathNet.Numerics.Random /// If the seed value is zero, it is set to one. Uses the /// value of to /// set whether the instance is thread safe. - public Palf() : this((int)DateTime.Now.Ticks) + public Palf() + : this((int) DateTime.Now.Ticks) { } @@ -76,7 +77,8 @@ namespace MathNet.Numerics.Random /// the current time as the seed. /// /// if set to true , the class is thread safe. - public Palf(bool threadSafe) : this((int)DateTime.Now.Ticks, threadSafe, DefaultShortLag, DefaultLongLag) + public Palf(bool threadSafe) + : this((int) DateTime.Now.Ticks, threadSafe, DefaultShortLag, DefaultLongLag) { } @@ -87,7 +89,8 @@ namespace MathNet.Numerics.Random /// If the seed value is zero, it is set to one. Uses the /// value of to /// set whether the instance is thread safe. - public Palf(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators, DefaultShortLag, DefaultLongLag) + public Palf(int seed) + : this(seed, Control.ThreadSafeRandomNumberGenerators, DefaultShortLag, DefaultLongLag) { } @@ -98,7 +101,8 @@ namespace MathNet.Numerics.Random /// if set to true, the class is thread safe. /// The ShortLag value /// TheLongLag value - public Palf(int seed, bool threadSafe, int shortLag, int longLag) : base(threadSafe) + public Palf(int seed, bool threadSafe, int shortLag, int longLag) + : base(threadSafe) { if (shortLag < 1) { @@ -116,54 +120,46 @@ namespace MathNet.Numerics.Random } ShortLag = shortLag; - + // Align LongLag to number of worker threads. - if (longLag % Control.NumberOfParallelWorkerThreads == 0) + if (longLag%Control.NumberOfParallelWorkerThreads == 0) { LongLag = longLag; } else { - LongLag = ((longLag / Control.NumberOfParallelWorkerThreads) + 1) * Control.NumberOfParallelWorkerThreads; + LongLag = ((longLag/Control.NumberOfParallelWorkerThreads) + 1)*Control.NumberOfParallelWorkerThreads; } - + _x = new uint[LongLag]; var gen = new MersenneTwister(seed, threadSafe); for (var j = 0; j < LongLag; ++j) { - _x[j] = (uint)(gen.NextDouble() * uint.MaxValue); + _x[j] = (uint) (gen.NextDouble()*uint.MaxValue); } _i = LongLag; } - + /// /// Gets the short lag of the Lagged Fibonacci pseudo-random number generator. /// - public int ShortLag - { - get; - private set; - } - + public int ShortLag { get; private set; } + /// /// Gets the long lag of the Lagged Fibonacci pseudo-random number generator. /// - public int LongLag - { - get; - private set; - } + public int LongLag { get; private set; } /// /// Stores an array of random numbers /// - private readonly uint[] _x; + readonly uint[] _x; /// /// Stores an index for the random number array element that will be accessed next. /// - private int _i; + int _i; /// /// Fills the array with new unsigned random numbers. @@ -172,22 +168,22 @@ namespace MathNet.Numerics.Random /// Generated random numbers are 32-bit unsigned integers greater than or equal to 0 /// and less than or equal to . /// - private void Fill() + void Fill() { - CommonParallel.For( - 0, - Control.NumberOfParallelWorkerThreads, - index => + CommonParallel.For(0, Control.NumberOfParallelWorkerThreads, (u, v) => { - // Two loops to avoid costly modulo operations - for (var j = index; j < ShortLag; j = j + Control.NumberOfParallelWorkerThreads) - { - _x[j] += _x[j + (LongLag - ShortLag)]; - } - - for (var j = ShortLag + index; j < LongLag; j = j + Control.NumberOfParallelWorkerThreads) + for (int index = u; index < v; index++) { - _x[j] += _x[j - ShortLag - index]; + // Two loops to avoid costly modulo operations + for (var j = index; j < ShortLag; j = j + Control.NumberOfParallelWorkerThreads) + { + _x[j] += _x[j + (LongLag - ShortLag)]; + } + + for (var j = ShortLag + index; j < LongLag; j = j + Control.NumberOfParallelWorkerThreads) + { + _x[j] += _x[j - ShortLag - index]; + } } }); _i = 0; @@ -207,7 +203,7 @@ namespace MathNet.Numerics.Random } var x = _x[_i++]; - return (int)(x >> 1) * IntToDoubleMultiplier; + return (int) (x >> 1)*IntToDoubleMultiplier; } } } diff --git a/src/Numerics/Threading/CommonParallel.cs b/src/Numerics/Threading/CommonParallel.cs index 70050612..9e4399ac 100644 --- a/src/Numerics/Threading/CommonParallel.cs +++ b/src/Numerics/Threading/CommonParallel.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2012 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -38,7 +38,7 @@ namespace MathNet.Numerics.Threading using System.Collections.Generic; #else using System.Linq; - using MathNet.Numerics.Properties; + using Properties; #endif /// @@ -54,80 +54,91 @@ namespace MathNet.Numerics.Threading /// The body to be invoked for each iteration. /// The argument is null. /// At least one invocation of the body threw an exception. + [Obsolete("Scheduled for removal in v3.0.")] public static void For(int fromInclusive, int toExclusive, Action body) { - if (body == null) - { - throw new ArgumentNullException("body"); - } + if (body == null) throw new ArgumentNullException("body"); + if (toExclusive <= fromInclusive) throw new ArgumentOutOfRangeException("toExclusive"); - // Special case: no action - if (fromInclusive >= toExclusive) - { - return; - } + int rangeSize = (toExclusive - fromInclusive)/(Control.NumberOfParallelWorkerThreads*2); + rangeSize = Math.Max(rangeSize, 1); - // Special case: single action, inline - if (fromInclusive == (toExclusive - 1)) + For(fromInclusive, + toExclusive, + rangeSize, + (start, stop) => + { + for (var i = start; i < stop; i++) + { + body(i); + } + }); + } + + /// + /// Executes a for loop in which iterations may run in parallel. + /// + /// The start index, inclusive. + /// The end index, exclusive. + /// The body to be invoked for each iteration range. + public static void For(int fromInclusive, int toExclusive, Action body) + { + For(fromInclusive, toExclusive, Math.Max(1, (toExclusive - fromInclusive)/Control.NumberOfParallelWorkerThreads), body); + } + + /// + /// Executes a for loop in which iterations may run in parallel. + /// + /// The start index, inclusive. + /// The end index, exclusive. + /// The body to be invoked for each iteration range. + public static void For(int fromInclusive, int toExclusive, int rangeSize, Action body) + { + if (body == null) throw new ArgumentNullException("body"); + if (fromInclusive < 0) throw new ArgumentOutOfRangeException("fromInclusive"); + if (fromInclusive > toExclusive) throw new ArgumentOutOfRangeException("toExclusive"); + if (rangeSize < 1) throw new ArgumentOutOfRangeException("rangeSize"); + + var length = toExclusive - fromInclusive; + + // Special case: nothing to do + if (length <= 0) { - body(fromInclusive); return; } - // Special case: straight execution without parallelism - if (Control.DisableParallelization || Control.NumberOfParallelWorkerThreads < 2) + var maxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads; + + // Special case: not worth to parallelize, inline + if (Control.DisableParallelization || maxDegreeOfParallelism < 2 || (rangeSize*2) > length) { - for (var index = fromInclusive; index < toExclusive; index++) - { - body(index); - } + body(fromInclusive, toExclusive); return; } - // Common case #if PORTABLE - var tasks = new Task[Control.NumberOfParallelWorkerThreads]; - var size = (toExclusive - fromInclusive) / tasks.Length; + var tasks = new Task[Math.Min(maxDegreeOfParallelism, length/rangeSize)]; + rangeSize = (toExclusive - fromInclusive)/tasks.Length; // partition the jobs into separate sets for each but the last worked thread for (var i = 0; i < tasks.Length - 1; i++) { - var start = fromInclusive + (i * size); - var stop = fromInclusive + ((i + 1) * size); + var start = fromInclusive + (i*rangeSize); + var stop = fromInclusive + ((i + 1)*rangeSize); - tasks[i] = Task.Factory.StartNew(() => - { - for (int j = start; j < stop; j++) - { - body(j); - } - }); + tasks[i] = Task.Factory.StartNew(() => body(start, stop)); } // add another set for last worker thread - tasks[tasks.Length - 1] = Task.Factory.StartNew(() => - { - for (int j = fromInclusive + ((tasks.Length - 1) * size); j < toExclusive; j++) - { - body(j); - } - }); + tasks[tasks.Length - 1] = + Task.Factory.StartNew(() => body(fromInclusive + ((tasks.Length - 1)*rangeSize), toExclusive)); Task.WaitAll(tasks); #else Parallel.ForEach( - Partitioner.Create(fromInclusive, toExclusive), - new ParallelOptions - { - MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads - }, - (range, loopState) => - { - for (var i = range.Item1; i < range.Item2; i++) - { - body(i); - } - }); + Partitioner.Create(fromInclusive, toExclusive, rangeSize), + new ParallelOptions {MaxDegreeOfParallelism = maxDegreeOfParallelism}, + (range, loopState) => body(range.Item1, range.Item2)); #endif } @@ -203,16 +214,16 @@ namespace MathNet.Numerics.Threading Parallel.ForEach( Partitioner.Create(0, array.Length), new ParallelOptions - { - MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads - }, + { + MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads + }, (range, loopState) => - { - for (var i = range.Item1; i < range.Item2; i++) { - body(i, array[i]); - } - }); + for (var i = range.Item1; i < range.Item2; i++) + { + body(i, array[i]); + } + }); #endif } @@ -237,6 +248,16 @@ namespace MathNet.Numerics.Threading return; } + // Special case: straight execution without parallelism + if (Control.DisableParallelization || Control.NumberOfParallelWorkerThreads < 2) + { + for (int i = 0; i < actions.Length; i++) + { + actions[i](); + } + return; + } + // Common case #if PORTABLE var tasks = new Task[actions.Length]; @@ -252,12 +273,11 @@ namespace MathNet.Numerics.Threading } Task.WaitAll(tasks); #else - var maxThreads = Control.DisableParallelization ? 1 : Control.NumberOfParallelWorkerThreads; Parallel.Invoke( new ParallelOptions - { - MaxDegreeOfParallelism = maxThreads - }, + { + MaxDegreeOfParallelism = Control.NumberOfParallelWorkerThreads + }, actions); #endif } @@ -290,7 +310,7 @@ namespace MathNet.Numerics.Threading // Special case: single action, inline if (fromInclusive == (toExclusive - 1)) { - return reduce(new [] { select(fromInclusive) }); + return reduce(new[] {select(fromInclusive)}); } // Special case: straight execution without parallelism @@ -346,25 +366,25 @@ namespace MathNet.Numerics.Threading var maxThreads = Control.DisableParallelization ? 1 : Control.NumberOfParallelWorkerThreads; Parallel.ForEach( Partitioner.Create(fromInclusive, toExclusive), - new ParallelOptions { MaxDegreeOfParallelism = maxThreads }, + new ParallelOptions {MaxDegreeOfParallelism = maxThreads}, () => new List(), (range, loop, localData) => - { - var mapped = new T[range.Item2 - range.Item1]; - for (int k = 0; k < mapped.Length; k++) { - mapped[k] = select(k + range.Item1); - } - localData.Add(reduce(mapped)); - return localData; - }, + var mapped = new T[range.Item2 - range.Item1]; + for (int k = 0; k < mapped.Length; k++) + { + mapped[k] = select(k + range.Item1); + } + localData.Add(reduce(mapped)); + return localData; + }, localResult => - { - lock (syncLock) { - intermediateResults.Add(reduce(localResult.ToArray())); - } - }); + lock (syncLock) + { + intermediateResults.Add(reduce(localResult.ToArray())); + } + }); return reduce(intermediateResults.ToArray()); #endif } @@ -396,7 +416,7 @@ namespace MathNet.Numerics.Threading // Special case: single action, inline if (array.Length == 1) { - return reduce(new[] { select(0, array[0]) }); + return reduce(new[] {select(0, array[0])}); } // Special case: straight execution without parallelism @@ -452,25 +472,25 @@ namespace MathNet.Numerics.Threading var maxThreads = Control.DisableParallelization ? 1 : Control.NumberOfParallelWorkerThreads; Parallel.ForEach( Partitioner.Create(0, array.Length), - new ParallelOptions { MaxDegreeOfParallelism = maxThreads }, + new ParallelOptions {MaxDegreeOfParallelism = maxThreads}, () => new List(), (range, loop, localData) => - { - var mapped = new U[range.Item2 - range.Item1]; - for (int k = 0; k < mapped.Length; k++) { - mapped[k] = select(k + range.Item1, array[k + range.Item1]); - } - localData.Add(reduce(mapped)); - return localData; - }, + var mapped = new U[range.Item2 - range.Item1]; + for (int k = 0; k < mapped.Length; k++) + { + mapped[k] = select(k + range.Item1, array[k + range.Item1]); + } + localData.Add(reduce(mapped)); + return localData; + }, localResult => - { - lock (syncLock) { - intermediateResults.Add(reduce(localResult.ToArray())); - } - }); + lock (syncLock) + { + intermediateResults.Add(reduce(localResult.ToArray())); + } + }); return reduce(intermediateResults.ToArray()); #endif } @@ -518,24 +538,24 @@ namespace MathNet.Numerics.Threading public static U Aggregate(T[] array, Func select, Func reducePair, U reduceDefault) { return Aggregate(array, select, results => - { - if (results == null || results.Length == 0) { - return reduceDefault; - } + if (results == null || results.Length == 0) + { + return reduceDefault; + } - if (results.Length == 1) - { - return results[0]; - } + if (results.Length == 1) + { + return results[0]; + } - U result = results[0]; - for (int i = 1; i < results.Length; i++) - { - result = reducePair(result, results[i]); - } - return result; - }); + U result = results[0]; + for (int i = 1; i < results.Length; i++) + { + result = reducePair(result, results[i]); + } + return result; + }); } } } \ No newline at end of file