uo9xpb74
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

tweaked cholesky solve 

fixed issue with LU solve for complex types using parallel element copy instead of Array.Copy
la-knuth
Marcus Cuda 15 years ago
parent
57da4a187e
commit
08f1483e87
5 changed files with 353 additions and 279 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 152
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  2. 152
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  3. 157
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  4. 154
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  5. 17
      src/Numerics/Control.cs

152
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
View File

@ -67,35 +67,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
if (alpha.IsZero()) if (alpha.IsZero())
{ {
Array.Copy(y, result, y.Length); y.Copy(result);
} }
else if (alpha.IsOne()) 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + x[index]; result[index] = y[index] + x[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
}
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + (alpha * x[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()) else if (alpha.IsOne())
{ {
Array.Copy(x, result, x.Length); x.Copy(result);
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = alpha * x[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); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] + y[index]; result[index] = x[index] + y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
}
} }
/// <summary> /// <summary>
@ -247,17 +247,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] - y[index]; result[index] = x[index] - y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; });
}
} }
/// <summary> /// <summary>
@ -292,17 +292,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] * y[index]; result[index] = x[index] * y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
}
} }
/// <summary> /// <summary>
@ -337,17 +337,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] / y[index]; result[index] = x[index] / y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
}
} }
/// <summary> /// <summary>
@ -1065,7 +1065,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
LUSolveFactored(order, a, order, ipiv, inverse); LUSolveFactored(order, a, order, ipiv, inverse);
CommonParallel.For(0, a.Length, index => a[index] = inverse[index]); inverse.Copy(a);
} }
/// <summary> /// <summary>
@ -1134,7 +1134,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var ipiv = new int[order]; var ipiv = new int[order];
var clone = new Complex[a.Length]; var clone = new Complex[a.Length];
Array.Copy(a, 0, clone, 0, a.Length); a.Copy(clone);
LUFactor(clone, order, ipiv); LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b); LUSolveFactored(columnsOfB, clone, order, ipiv, b);
} }
@ -1357,7 +1357,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new Complex[a.Length]; var clone = new Complex[a.Length];
Array.Copy(a, 0, clone, 0, a.Length); a.Copy(clone);
CholeskyFactor(clone, orderA); CholeskyFactor(clone, orderA);
CholeskySolveFactored(clone, orderA, b, columnsB); CholeskySolveFactored(clone, orderA, b, columnsB);
} }
@ -1392,39 +1392,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentReferenceDifferent); throw new ArgumentException(Resources.ArgumentReferenceDifferent);
} }
CommonParallel.For( if (Control.ParallelizeOperation(columnsB * 10))
0, {
columnsB, CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c));
c => }
else
{
for (var index = 0; index < columnsB; index++)
{ {
var cindex = c * orderA; DoCholeskySolve(a, orderA, b, index);
}
}
}
// Solve L*Y = B; /// <summary>
Complex sum; /// Solves A*X=B for X using a previously factored A matrix.
for (var i = 0; i < orderA; i++) /// </summary>
{ /// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
sum = b[cindex + i]; /// <param name="orderA">The number of rows and columns in A.</param>
for (var k = i - 1; k >= 0; k--) /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
{ /// <param name="index">The column to solve for.</param>
sum -= a[(k * orderA) + i] * b[cindex + k]; 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; b[cindex + i] = sum / a[(i * orderA) + i];
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]; // 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];
}
} }
/// <summary> /// <summary>

152
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
View File

@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
if (alpha.IsZero()) if (alpha.IsZero())
{ {
Array.Copy(y, result, y.Length); y.Copy(result);
} }
else if (alpha.IsOne()) 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + x[index]; result[index] = y[index] + x[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
}
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + (alpha * x[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()) else if (alpha.IsOne())
{ {
Array.Copy(x, result, x.Length); x.Copy(result);
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = alpha * x[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); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] + y[index]; result[index] = x[index] + y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
}
} }
/// <summary> /// <summary>
@ -243,17 +243,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] - y[index]; result[index] = x[index] - y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; });
}
} }
/// <summary> /// <summary>
@ -288,17 +288,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] * y[index]; result[index] = x[index] * y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
}
} }
/// <summary> /// <summary>
@ -333,17 +333,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] / y[index]; result[index] = x[index] / y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
}
} }
/// <summary> /// <summary>
@ -1061,7 +1061,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
LUSolveFactored(order, a, order, ipiv, inverse); LUSolveFactored(order, a, order, ipiv, inverse);
CommonParallel.For(0, a.Length, index => a[index] = inverse[index]); inverse.Copy(a);
} }
/// <summary> /// <summary>
@ -1130,7 +1130,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var ipiv = new int[order]; var ipiv = new int[order];
var clone = new Complex32[a.Length]; var clone = new Complex32[a.Length];
Array.Copy(a, 0, clone, 0, a.Length); a.Copy(clone);
LUFactor(clone, order, ipiv); LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b); LUSolveFactored(columnsOfB, clone, order, ipiv, b);
} }
@ -1353,7 +1353,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new Complex32[a.Length]; var clone = new Complex32[a.Length];
Array.Copy(a, 0, clone, 0, a.Length); a.Copy(clone);
CholeskyFactor(clone, orderA); CholeskyFactor(clone, orderA);
CholeskySolveFactored(clone, orderA, b, columnsB); CholeskySolveFactored(clone, orderA, b, columnsB);
} }
@ -1388,39 +1388,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentReferenceDifferent); throw new ArgumentException(Resources.ArgumentReferenceDifferent);
} }
CommonParallel.For( if (Control.ParallelizeOperation(columnsB * 10))
0, {
columnsB, CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c));
c => }
else
{
for (var index = 0; index < columnsB; index++)
{ {
var cindex = c * orderA; DoCholeskySolve(a, orderA, b, index);
}
}
}
// Solve L*Y = B; /// <summary>
Complex32 sum; /// Solves A*X=B for X using a previously factored A matrix.
for (var i = 0; i < orderA; i++) /// </summary>
{ /// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
sum = b[cindex + i]; /// <param name="orderA">The number of rows and columns in A.</param>
for (var k = i - 1; k >= 0; k--) /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
{ /// <param name="index">The column to solve for.</param>
sum -= a[(k * orderA) + i] * b[cindex + k]; 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; b[cindex + i] = sum / a[(i * orderA) + i];
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]; // 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];
}
} }
/// <summary> /// <summary>

157
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
View File

@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
if (alpha == 0.0) if (alpha == 0.0)
{ {
Buffer.BlockCopy(y, 0, result, 0, y.Length * Constants.SizeOfDouble); y.Copy(result);
} }
else if (alpha == 1.0) 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + x[index]; result[index] = y[index] + x[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
}
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + (alpha * x[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) else if (alpha == 1.0)
{ {
Buffer.BlockCopy(x, 0, result, 0, x.Length * Constants.SizeOfDouble); x.Copy(result);
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = alpha * x[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); 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++) CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
{
result[index] = x[index] + y[index];
}
} }
else 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); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] - y[index]; result[index] = x[index] - y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; });
}
} }
/// <summary> /// <summary>
@ -287,17 +287,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] * y[index]; result[index] = x[index] * y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
}
} }
/// <summary> /// <summary>
@ -332,17 +332,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] / y[index]; result[index] = x[index] / y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
}
} }
/// <summary> /// <summary>
@ -950,7 +950,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
LUSolveFactored(order, a, order, ipiv, inverse); LUSolveFactored(order, a, order, ipiv, inverse);
Buffer.BlockCopy(inverse, 0, a, 0, a.Length * Constants.SizeOfDouble); inverse.Copy(a);
} }
/// <summary> /// <summary>
@ -1019,7 +1019,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var ipiv = new int[order]; var ipiv = new int[order];
var clone = new double[a.Length]; var clone = new double[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); a.Copy(clone);
LUFactor(clone, order, ipiv); LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b); LUSolveFactored(columnsOfB, clone, order, ipiv, b);
} }
@ -1242,7 +1242,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new double[a.Length]; var clone = new double[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); a.Copy(clone);
CholeskyFactor(clone, orderA); CholeskyFactor(clone, orderA);
CholeskySolveFactored(clone, orderA, b, columnsB); CholeskySolveFactored(clone, orderA, b, columnsB);
} }
@ -1277,40 +1277,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentReferenceDifferent); throw new ArgumentException(Resources.ArgumentReferenceDifferent);
} }
if( ) if (Control.ParallelizeOperation(columnsB * 10))
CommonParallel.For( {
0, CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c));
columnsB, }
c => else
{
for (var index = 0; index < columnsB; index++)
{ {
var cindex = c * orderA; DoCholeskySolve(a, orderA, b, index);
}
}
}
// Solve L*Y = B; /// <summary>
double sum; /// Solves A*X=B for X using a previously factored A matrix.
for (var i = 0; i < orderA; i++) /// </summary>
{ /// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
sum = b[cindex + i]; /// <param name="orderA">The number of rows and columns in A.</param>
for (var k = i - 1; k >= 0; k--) /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
{ /// <param name="index">The column to solve for.</param>
sum -= a[(k * orderA) + i] * b[cindex + k]; 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; b[cindex + i] = sum / a[(i * orderA) + i];
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]; // 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];
}
} }
/// <summary> /// <summary>
@ -1642,7 +1657,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new double[a.Length]; 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]; var q = new double[rows * rows];
QRFactor(clone, rows, columns, q, work); QRFactor(clone, rows, columns, q, work);
QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x);

154
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
View File

@ -61,35 +61,35 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
if (alpha == 0.0) if (alpha == 0.0)
{ {
Buffer.BlockCopy(y, 0, result, 0, y.Length * Constants.SizeOfFloat); y.Copy(result);
} }
else if (alpha == 1.0) 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + x[index]; result[index] = y[index] + x[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
}
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = y[index] + (alpha * x[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) else if (alpha == 1.0)
{ {
Buffer.BlockCopy(x, 0, result, 0, x.Length * Constants.SizeOfFloat); x.Copy(result);
} }
else 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = alpha * x[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); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] + y[index]; result[index] = x[index] + y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
}
} }
/// <summary> /// <summary>
@ -242,17 +242,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] - y[index]; result[index] = x[index] - y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; });
}
} }
/// <summary> /// <summary>
@ -287,17 +287,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] * y[index]; result[index] = x[index] * y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
}
} }
/// <summary> /// <summary>
@ -332,17 +332,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentVectorsSameLength); 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++) for (var index = 0; index < x.Length; index++)
{ {
result[index] = x[index] / y[index]; result[index] = x[index] / y[index];
} }
} }
else
{
CommonParallel.For(0, y.Length, index => { result[index] = x[index] / y[index]; });
}
} }
/// <summary> /// <summary>
@ -951,7 +951,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
LUSolveFactored(order, a, order, ipiv, inverse); LUSolveFactored(order, a, order, ipiv, inverse);
Buffer.BlockCopy(inverse, 0, a, 0, a.Length * Constants.SizeOfFloat); inverse.Copy(a);
} }
/// <summary> /// <summary>
@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var ipiv = new int[order]; var ipiv = new int[order];
var clone = new float[a.Length]; var clone = new float[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); a.Copy(clone);
LUFactor(clone, order, ipiv); LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b); LUSolveFactored(columnsOfB, clone, order, ipiv, b);
} }
@ -1243,7 +1243,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new float[a.Length]; var clone = new float[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); a.Copy(clone);
CholeskyFactor(clone, orderA); CholeskyFactor(clone, orderA);
CholeskySolveFactored(clone, orderA, b, columnsB); CholeskySolveFactored(clone, orderA, b, columnsB);
} }
@ -1278,39 +1278,55 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentException(Resources.ArgumentReferenceDifferent); throw new ArgumentException(Resources.ArgumentReferenceDifferent);
} }
CommonParallel.For( if (Control.ParallelizeOperation(columnsB * 10))
0, {
columnsB, CommonParallel.For(0, columnsB, c => DoCholeskySolve(a, orderA, b, c));
c => }
else
{
for (var index = 0; index < columnsB; index++)
{ {
var cindex = c * orderA; DoCholeskySolve(a, orderA, b, index);
}
}
}
// Solve L*Y = B; /// <summary>
float sum; /// Solves A*X=B for X using a previously factored A matrix.
for (var i = 0; i < orderA; i++) /// </summary>
{ /// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
sum = b[cindex + i]; /// <param name="orderA">The number of rows and columns in A.</param>
for (var k = i - 1; k >= 0; k--) /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
{ /// <param name="index">The column to solve for.</param>
sum -= a[(k * orderA) + i] * b[cindex + k]; 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; b[cindex + i] = sum / a[(i * orderA) + i];
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]; // 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];
}
} }
/// <summary> /// <summary>
@ -1642,7 +1658,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
} }
var clone = new float[a.Length]; 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]; var q = new float[rows * rows];
QRFactor(clone, rows, columns, q, work); QRFactor(clone, rows, columns, q, work);
QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x);

17
src/Numerics/Control.cs
View File

@ -59,7 +59,7 @@ namespace MathNet.Numerics
/// <summary> /// <summary>
/// The default cutoff point for order size for the matrix multiply in linear algebra provider. /// The default cutoff point for order size for the matrix multiply in linear algebra provider.
/// </summary> /// </summary>
private static int _parallelizeElements = 256; private static int _parallelizeElements = 300;
/// <summary> /// <summary>
/// Initializes static members of the Control class. /// Initializes static members of the Control class.
@ -99,7 +99,8 @@ namespace MathNet.Numerics
/// <value>The linear algebra provider.</value> /// <value>The linear algebra provider.</value>
public static ILinearAlgebraProvider LinearAlgebraProvider public static ILinearAlgebraProvider LinearAlgebraProvider
{ {
get; set; get;
set;
} }
/// <summary> /// <summary>
@ -169,7 +170,7 @@ namespace MathNet.Numerics
/// <summary> /// <summary>
/// Gets or sets the number of elements a vector or matrix must contain before we multiply threads. /// Gets or sets the number of elements a vector or matrix must contain before we multiply threads.
/// </summary> /// </summary>
/// <value>Number of elements. Default is 256.</value> /// <value>Number of elements. Default is 300.</value>
public static int ParallelizeElements public static int ParallelizeElements
{ {
get get
@ -185,5 +186,15 @@ namespace MathNet.Numerics
} }
} }
} }
/// <summary>
/// Given the number elements, should the operation be parallelized.
/// </summary>
/// <param name="elements">The number elements to check.</param>
/// <returns><c>true</c> if the operation should be parallelized; <c>false</c> otherwise.</returns>
public static bool ParallelizeOperation(int elements)
{
return elements < ParallelizeElements || DisableParallelization || NumberOfParallelWorkerThreads < 2;
}
} }
} }

Write Preview
Loading…
Cancel
Save