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

LA: managed providers drop variation argument

spatial
Christoph Ruegg 9 years ago
parent
f9e4312c6d
commit
7e24e98f39
10 changed files with 226 additions and 2040 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 482
      src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs
  2. 482
      src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs
  3. 371
      src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs
  4. 371
      src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs
  5. 12
      src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs
  6. 134
      src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs
  7. 134
      src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs
  8. 134
      src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs
  9. 134
      src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs
  10. 12
      src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs

482
src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex.cs
View File

@ -468,480 +468,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
throw new ArgumentNullException("x");
throw new ArgumentNullException("a");
}
if (y == null)
{
throw new ArgumentNullException("y");
throw new ArgumentNullException("b");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (rowsX*columnsX != x.Length)
{
throw new ArgumentException("x.Length != xRows * xColumns");
}
if (rowsY*columnsY != y.Length)
{
throw new ArgumentException("y.Length != yRows * yColumns");
throw new ArgumentNullException("c");
}
if (columnsX != rowsY)
{
throw new ArgumentException("xColumns != yRows");
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY));
}
if (rowsX*columnsY != result.Length)
if (rowsX * columnsX != x.Length)
{
throw new ArgumentException("xRows * yColumns != result.Length");
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length));
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
Complex[] xdata;
if (ReferenceEquals(x, result))
if (rowsY * columnsY != y.Length)
{
xdata = (Complex[]) x.Clone();
}
else
{
xdata = x;
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length));
}
Complex[] ydata;
if (ReferenceEquals(y, result))
{
ydata = (Complex[]) y.Clone();
}
else
if (rowsX * columnsY != result.Length)
{
ydata = y;
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length));
}
// handle degenerate cases
Array.Clear(result, 0, result.Length);
CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
// First check some basic requirement on the parameters of the matrix multiplication.
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if ((int) transposeA > 111 && (int) transposeB > 111)
{
if (rowsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = rowsB;
k = rowsA;
}
else if ((int) transposeA > 111)
{
if (rowsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = columnsB;
k = rowsA;
}
else if ((int) transposeB > 111)
{
if (columnsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = rowsB;
k = columnsA;
}
else
{
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = columnsB;
k = columnsA;
}
if (alpha.IsZero() && beta.IsZero())
{
Array.Clear(c, 0, c.Length);
return;
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
Complex[] adata;
if (ReferenceEquals(a, c))
{
adata = (Complex[]) a.Clone();
}
else
{
adata = a;
}
Complex[] bdata;
if (ReferenceEquals(b, c))
{
bdata = (Complex[]) b.Clone();
}
else
{
bdata = b;
}
if (beta.IsZero())
{
Array.Clear(c, 0, c.Length);
}
else if (!beta.IsOne())
{
ScaleArray(beta, c, c);
}
if (alpha.IsZero())
// Extract column arrays
var columnDataB = new Complex[columnsY][];
for (int i = 0; i < columnDataB.Length; i++)
{
return;
var column = new Complex[rowsY];
GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column);
columnDataB[i] = column;
}
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true);
}
/// <summary>
/// Cache-Oblivious Matrix Multiplication
/// </summary>
/// <param name="transposeA">if set to <c>true</c> transpose matrix A.</param>
/// <param name="transposeB">if set to <c>true</c> transpose matrix B.</param>
/// <param name="alpha">The value to scale the matrix A with.</param>
/// <param name="matrixA">The matrix A.</param>
/// <param name="shiftArow">Row-shift of the left matrix</param>
/// <param name="shiftAcol">Column-shift of the left matrix</param>
/// <param name="matrixB">The matrix B.</param>
/// <param name="shiftBrow">Row-shift of the right matrix</param>
/// <param name="shiftBcol">Column-shift of the right matrix</param>
/// <param name="result">The matrix C.</param>
/// <param name="shiftCrow">Row-shift of the result matrix</param>
/// <param name="shiftCcol">Column-shift of the result matrix</param>
/// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="k">The number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="constM">The constant number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="first">Indicates if this is the first recursion.</param>
static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] matrixA, int shiftArow, int shiftAcol, Complex[] matrixB, int shiftBrow, int shiftBcol, Complex[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
{
if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1)
var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
if ((int) transposeA > 111 && (int) transposeB > 111)
var row = new Complex[columnsX];
for (int i = 0; i < rowsX; i++)
{
if ((int) transposeA > 112 && (int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeA > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
}
else if ((int) transposeA > 111)
{
if ((int) transposeA > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
}
else if ((int) transposeB > 111)
{
if ((int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
var col = columnDataB[j];
Complex sum = Complex.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
sum += row[ii] * col[ii];
}
}
else
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
result[j * rowsX + i] += Complex.One * sum;
}
}
else
}
else
{
CommonParallel.For(0, rowsX, 1, (u, v) =>
{
for (var m1 = 0; m1 < m; m1++)
var row = new Complex[columnsX];
for (int i = u; i < v; i++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex.Zero;
for (var k1 = 0; k1 < k; ++k1)
var column = columnDataB[j];
Complex sum = Complex.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
sum += row[ii] * column[ii];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
result[j * rowsX + i] += Complex.One * sum;
}
}
}
}
else
{
// divide and conquer
int m2 = m/2, n2 = n/2, k2 = k/2;
if (first)
{
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false));
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false));
}
else
{
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false);
}
});
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA,
Complex[] b,
int rowsB, int columnsB, Complex beta, Complex[] c)
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c)
{
if (a == null)
{

482
src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Complex32.cs
View File

@ -470,480 +470,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
throw new ArgumentNullException("x");
throw new ArgumentNullException("a");
}
if (y == null)
{
throw new ArgumentNullException("y");
throw new ArgumentNullException("b");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (rowsX*columnsX != x.Length)
{
throw new ArgumentException("x.Length != xRows * xColumns");
}
if (rowsY*columnsY != y.Length)
{
throw new ArgumentException("y.Length != yRows * yColumns");
throw new ArgumentNullException("c");
}
if (columnsX != rowsY)
{
throw new ArgumentException("xColumns != yRows");
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY));
}
if (rowsX*columnsY != result.Length)
if (rowsX * columnsX != x.Length)
{
throw new ArgumentException("xRows * yColumns != result.Length");
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length));
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
Complex32[] xdata;
if (ReferenceEquals(x, result))
if (rowsY * columnsY != y.Length)
{
xdata = (Complex32[]) x.Clone();
}
else
{
xdata = x;
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length));
}
Complex32[] ydata;
if (ReferenceEquals(y, result))
{
ydata = (Complex32[]) y.Clone();
}
else
if (rowsX * columnsY != result.Length)
{
ydata = y;
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length));
}
// handle degenerate cases
Array.Clear(result, 0, result.Length);
CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
// First check some basic requirement on the parameters of the matrix multiplication.
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if ((int) transposeA > 111 && (int) transposeB > 111)
{
if (rowsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = rowsB;
k = rowsA;
}
else if ((int) transposeA > 111)
{
if (rowsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = columnsB;
k = rowsA;
}
else if ((int) transposeB > 111)
{
if (columnsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = rowsB;
k = columnsA;
}
else
{
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = columnsB;
k = columnsA;
}
if (alpha.IsZero() && beta.IsZero())
{
Array.Clear(c, 0, c.Length);
return;
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
Complex32[] adata;
if (ReferenceEquals(a, c))
{
adata = (Complex32[]) a.Clone();
}
else
{
adata = a;
}
Complex32[] bdata;
if (ReferenceEquals(b, c))
{
bdata = (Complex32[]) b.Clone();
}
else
{
bdata = b;
}
if (beta.IsZero())
{
Array.Clear(c, 0, c.Length);
}
else if (!beta.IsOne())
{
ScaleArray(beta, c, c);
}
if (alpha.IsZero())
// Extract column arrays
var columnDataB = new Complex32[columnsY][];
for (int i = 0; i < columnDataB.Length; i++)
{
return;
var column = new Complex32[rowsY];
GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column);
columnDataB[i] = column;
}
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true);
}
/// <summary>
/// Cache-Oblivious Matrix Multiplication
/// </summary>
/// <param name="transposeA">if set to <c>true</c> transpose matrix A.</param>
/// <param name="transposeB">if set to <c>true</c> transpose matrix B.</param>
/// <param name="alpha">The value to scale the matrix A with.</param>
/// <param name="matrixA">The matrix A.</param>
/// <param name="shiftArow">Row-shift of the left matrix</param>
/// <param name="shiftAcol">Column-shift of the left matrix</param>
/// <param name="matrixB">The matrix B.</param>
/// <param name="shiftBrow">Row-shift of the right matrix</param>
/// <param name="shiftBcol">Column-shift of the right matrix</param>
/// <param name="result">The matrix C.</param>
/// <param name="shiftCrow">Row-shift of the result matrix</param>
/// <param name="shiftCcol">Column-shift of the result matrix</param>
/// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="k">The number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="constM">The constant number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="first">Indicates if this is the first recursion.</param>
static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] matrixA, int shiftArow, int shiftAcol, Complex32[] matrixB, int shiftBrow, int shiftBcol, Complex32[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
{
if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1)
var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
if ((int) transposeA > 111 && (int) transposeB > 111)
var row = new Complex32[columnsX];
for (int i = 0; i < rowsX; i++)
{
if ((int) transposeA > 112 && (int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeA > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
}
else if ((int) transposeA > 111)
{
if ((int) transposeA > 112)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
}
else if ((int) transposeB > 111)
{
if ((int) transposeB > 112)
{
for (var m1 = 0; m1 < m; m1++)
var col = columnDataB[j];
Complex32 sum = Complex32.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
sum += row[ii] * col[ii];
}
}
else
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
result[j * rowsX + i] += Complex32.One * sum;
}
}
else
}
else
{
CommonParallel.For(0, rowsX, 1, (u, v) =>
{
for (var m1 = 0; m1 < m; m1++)
var row = new Complex32[columnsX];
for (int i = u; i < v; i++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matBcolPos = n1 + shiftBcol;
var sum = Complex32.Zero;
for (var k1 = 0; k1 < k; ++k1)
var column = columnDataB[j];
Complex32 sum = Complex32.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
sum += row[ii] * column[ii];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
result[j * rowsX + i] += Complex32.One * sum;
}
}
}
}
else
{
// divide and conquer
int m2 = m/2, n2 = n/2, k2 = k/2;
if (first)
{
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false));
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false));
}
else
{
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false);
}
});
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA,
Complex32[] b,
int rowsB, int columnsB, Complex32 beta, Complex32[] c)
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c)
{
if (a == null)
{

371
src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Double.cs
View File

@ -465,371 +465,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
throw new ArgumentNullException("x");
throw new ArgumentNullException("a");
}
if (y == null)
{
throw new ArgumentNullException("y");
throw new ArgumentNullException("b");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (rowsX * columnsX != x.Length)
{
throw new ArgumentException("x.Length != xRows * xColumns");
}
if (rowsY * columnsY != y.Length)
{
throw new ArgumentException("y.Length != yRows * yColumns");
throw new ArgumentNullException("c");
}
if (columnsX != rowsY)
{
throw new ArgumentException("xColumns != yRows");
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY));
}
if (rowsX * columnsY != result.Length)
if (rowsX * columnsX != x.Length)
{
throw new ArgumentException("xRows * yColumns != result.Length");
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length));
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
double[] xdata;
if (ReferenceEquals(x, result))
{
xdata = (double[])x.Clone();
}
else
if (rowsY * columnsY != y.Length)
{
xdata = x;
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length));
}
double[] ydata;
if (ReferenceEquals(y, result))
{
ydata = (double[])y.Clone();
}
else
if (rowsX * columnsY != result.Length)
{
ydata = y;
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length));
}
// handle degenerate cases
Array.Clear(result, 0, result.Length);
CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
// First check some basic requirement on the parameters of the matrix multiplication.
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if ((int)transposeA > 111 && (int)transposeB > 111)
{
if (rowsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA * rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = rowsB;
k = rowsA;
}
else if ((int)transposeA > 111)
{
if (rowsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = columnsB;
k = rowsA;
}
else if ((int)transposeB > 111)
{
if (columnsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA * rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = rowsB;
k = columnsA;
}
else
{
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = columnsB;
k = columnsA;
}
if (alpha == 0.0 && beta == 0.0)
{
Array.Clear(c, 0, c.Length);
return;
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
double[] adata;
if (ReferenceEquals(a, c))
{
adata = (double[])a.Clone();
}
else
{
adata = a;
}
double[] bdata;
if (ReferenceEquals(b, c))
{
bdata = (double[])b.Clone();
}
else
{
bdata = b;
}
if (beta == 0.0)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != 1.0)
{
ScaleArray(beta, c, c);
}
if (alpha == 0.0)
// Extract column arrays
var columnDataB = new double[columnsY][];
for (int i = 0; i < columnDataB.Length; i++)
{
return;
var column = new double[rowsY];
GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column);
columnDataB[i] = column;
}
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true);
}
/// <summary>
/// Cache-Oblivious Matrix Multiplication
/// </summary>
/// <param name="transposeA">if set to <c>true</c> transpose matrix A.</param>
/// <param name="transposeB">if set to <c>true</c> transpose matrix B.</param>
/// <param name="alpha">The value to scale the matrix A with.</param>
/// <param name="matrixA">The matrix A.</param>
/// <param name="shiftArow">Row-shift of the left matrix</param>
/// <param name="shiftAcol">Column-shift of the left matrix</param>
/// <param name="matrixB">The matrix B.</param>
/// <param name="shiftBrow">Row-shift of the right matrix</param>
/// <param name="shiftBcol">Column-shift of the right matrix</param>
/// <param name="result">The matrix C.</param>
/// <param name="shiftCrow">Row-shift of the result matrix</param>
/// <param name="shiftCcol">Column-shift of the result matrix</param>
/// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="k">The number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="constM">The constant number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="first">Indicates if this is the first recursion.</param>
static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, double alpha, double[] matrixA, int shiftArow, int shiftAcol, double[] matrixB, int shiftBrow, int shiftBcol, double[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
{
if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1)
var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
if ((int) transposeA > 111 && (int) transposeB > 111)
var row = new double[columnsX];
for (int i = 0; i < rowsX; i++)
{
for (var m1 = 0; m1 < m; m1++)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
double sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeA > 111)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
var col = columnDataB[j];
double sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
var matBcolPos = n1 + shiftBcol;
double sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
sum += row[ii] * col[ii];
}
}
}
else if ((int) transposeB > 111)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
double sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
result[j * rowsX + i] += 1.0 * sum;
}
}
else
}
else
{
CommonParallel.For(0, rowsX, 1, (u, v) =>
{
for (var m1 = 0; m1 < m; m1++)
var row = new double[columnsX];
for (int i = u; i < v; i++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matBcolPos = n1 + shiftBcol;
var column = columnDataB[j];
double sum = 0;
for (var k1 = 0; k1 < k; ++k1)
for (int ii = 0; ii < row.Length; ii++)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
sum += row[ii] * column[ii];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
result[j * rowsX + i] += 1.0 * sum;
}
}
}
}
else
{
// divide and conquer
int m2 = m/2, n2 = n/2, k2 = k/2;
if (first)
{
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false));
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false));
}
else
{
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false);
}
});
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA,
double[] b,
int rowsB, int columnsB, double beta, double[] c)
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c)
{
if (a == null)
{

371
src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.Single.cs
View File

@ -465,371 +465,112 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
throw new ArgumentNullException("x");
throw new ArgumentNullException("a");
}
if (y == null)
{
throw new ArgumentNullException("y");
throw new ArgumentNullException("b");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (rowsX*columnsX != x.Length)
{
throw new ArgumentException("x.Length != xRows * xColumns");
}
if (rowsY*columnsY != y.Length)
{
throw new ArgumentException("y.Length != yRows * yColumns");
throw new ArgumentNullException("c");
}
if (columnsX != rowsY)
{
throw new ArgumentException("xColumns != yRows");
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsX, rowsY));
}
if (rowsX*columnsY != result.Length)
if (rowsX * columnsX != x.Length)
{
throw new ArgumentException("xRows * yColumns != result.Length");
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsX, columnsX, x.Length));
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
float[] xdata;
if (ReferenceEquals(x, result))
{
xdata = (float[]) x.Clone();
}
else
if (rowsY * columnsY != y.Length)
{
xdata = x;
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsY, columnsY, y.Length));
}
float[] ydata;
if (ReferenceEquals(y, result))
{
ydata = (float[]) y.Clone();
}
else
if (rowsX * columnsY != result.Length)
{
ydata = y;
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsX, columnsY, result.Length));
}
// handle degenerate cases
Array.Clear(result, 0, result.Length);
CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
// First check some basic requirement on the parameters of the matrix multiplication.
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if ((int) transposeA > 111 && (int) transposeB > 111)
{
if (rowsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = rowsB;
k = rowsA;
}
else if ((int) transposeA > 111)
{
if (rowsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (columnsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = columnsA;
n = columnsB;
k = rowsA;
}
else if ((int) transposeB > 111)
{
if (columnsA != columnsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*rowsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = rowsB;
k = columnsA;
}
else
{
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException();
}
if (rowsA*columnsB != c.Length)
{
throw new ArgumentOutOfRangeException();
}
m = rowsA;
n = columnsB;
k = columnsA;
}
if (alpha == 0.0 && beta == 0.0)
{
Array.Clear(c, 0, c.Length);
return;
}
// Check whether we will be overwriting any of our inputs and make copies if necessary.
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
// as result, we can do it on a row wise basis. We should investigate this.
float[] adata;
if (ReferenceEquals(a, c))
{
adata = (float[]) a.Clone();
}
else
{
adata = a;
}
float[] bdata;
if (ReferenceEquals(b, c))
{
bdata = (float[]) b.Clone();
}
else
{
bdata = b;
}
if (beta == 0.0f)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != 1.0f)
{
ScaleArray(beta, c, c);
}
if (alpha == 0.0f)
// Extract column arrays
var columnDataB = new float[columnsY][];
for (int i = 0; i < columnDataB.Length; i++)
{
return;
var column = new float[rowsY];
GetColumn(Transpose.DontTranspose, i, rowsY, columnsY, y, column);
columnDataB[i] = column;
}
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true);
}
/// <summary>
/// Cache-Oblivious Matrix Multiplication
/// </summary>
/// <param name="transposeA">if set to <c>true</c> transpose matrix A.</param>
/// <param name="transposeB">if set to <c>true</c> transpose matrix B.</param>
/// <param name="alpha">The value to scale the matrix A with.</param>
/// <param name="matrixA">The matrix A.</param>
/// <param name="shiftArow">Row-shift of the left matrix</param>
/// <param name="shiftAcol">Column-shift of the left matrix</param>
/// <param name="matrixB">The matrix B.</param>
/// <param name="shiftBrow">Row-shift of the right matrix</param>
/// <param name="shiftBcol">Column-shift of the right matrix</param>
/// <param name="result">The matrix C.</param>
/// <param name="shiftCrow">Row-shift of the result matrix</param>
/// <param name="shiftCcol">Column-shift of the result matrix</param>
/// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="k">The number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="constM">The constant number of rows of matrix op(A) and of the matrix C.</param>
/// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
/// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
/// <param name="first">Indicates if this is the first recursion.</param>
static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, float alpha, float[] matrixA, int shiftArow, int shiftAcol, float[] matrixB, int shiftBrow, int shiftBcol, float[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
{
if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1)
var shouldNotParallelize = rowsX + columnsY + columnsX < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
if ((int) transposeA > 111 && (int) transposeB > 111)
var row = new float[columnsX];
for (int i = 0; i < rowsX; i++)
{
for (var m1 = 0; m1 < m; m1++)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
float sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
}
}
else if ((int) transposeA > 111)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
var col = columnDataB[j];
float sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
var matBcolPos = n1 + shiftBcol;
float sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
sum += row[ii] * col[ii];
}
}
}
else if ((int) transposeB > 111)
{
for (var m1 = 0; m1 < m; m1++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
{
var matBcolPos = n1 + shiftBcol;
float sum = 0;
for (var k1 = 0; k1 < k; ++k1)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
}
result[j * rowsX + i] += 1.0f * sum;
}
}
else
}
else
{
CommonParallel.For(0, rowsX, 1, (u, v) =>
{
for (var m1 = 0; m1 < m; m1++)
var row = new float[columnsX];
for (int i = u; i < v; i++)
{
var matArowPos = m1 + shiftArow;
var matCrowPos = m1 + shiftCrow;
for (var n1 = 0; n1 < n; ++n1)
GetRow(Transpose.DontTranspose, i, rowsX, columnsX, x, row);
for (int j = 0; j < columnsY; j++)
{
var matBcolPos = n1 + shiftBcol;
var column = columnDataB[j];
float sum = 0;
for (var k1 = 0; k1 < k; ++k1)
for (int ii = 0; ii < row.Length; ii++)
{
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
sum += row[ii] * column[ii];
}
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
result[j * rowsX + i] += 1.0f * sum;
}
}
}
}
else
{
// divide and conquer
int m2 = m/2, n2 = n/2, k2 = k/2;
if (first)
{
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false));
CommonParallel.Invoke(
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false),
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false));
}
else
{
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false);
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false);
}
});
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA,
float[] b,
int rowsB, int columnsB, float beta, float[] c)
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c)
{
if (a == null)
{

12
src/Numerics/Providers/LinearAlgebra/Managed/ManagedLinearAlgebraProvider.cs
View File

@ -42,18 +42,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Managed
/// </summary>
internal partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider
{
private readonly Variation _variation;
internal ManagedLinearAlgebraProvider()
{
_variation = Variation.Experimental;
}
internal ManagedLinearAlgebraProvider(Variation variation)
{
_variation = variation;
}
/// <summary>
/// Try to find out whether the provider is available, at least in principle.
/// Verification may still fail if available, but it will certainly fail if unavailable.

134
src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex.cs
View File

@ -468,12 +468,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
@ -554,12 +548,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@ -938,128 +926,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA,
Complex[] b,
int rowsB, int columnsB, Complex beta, Complex[] c)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (c == null)
{
throw new ArgumentNullException("c");
}
if (transposeA != Transpose.DontTranspose)
{
var swap = rowsA;
rowsA = columnsA;
columnsA = swap;
}
if (transposeB != Transpose.DontTranspose)
{
var swap = rowsB;
rowsB = columnsB;
columnsB = swap;
}
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsA, rowsB));
}
if (rowsA * columnsA != a.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsA, columnsA, a.Length));
}
if (rowsB * columnsB != b.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsB, columnsB, b.Length));
}
if (rowsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsA, columnsB, c.Length));
}
// handle degenerate cases
if (beta == Complex.Zero)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != Complex.One)
{
ScaleArray(beta, c, c);
}
if (alpha == Complex.Zero)
{
return;
}
// Extract column arrays
var columnDataB = new Complex[columnsB][];
for (int i = 0; i < columnDataB.Length; i++)
{
var column = new Complex[rowsB];
GetColumn(transposeB, i, rowsB, columnsB, b, column);
columnDataB[i] = column;
}
var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
var row = new Complex[columnsA];
for (int i = 0; i < rowsA; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var col = columnDataB[j];
Complex sum = Complex.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * col[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
}
else
{
CommonParallel.For(0, rowsA, 1, (u, v) =>
{
var row = new Complex[columnsA];
for (int i = u; i < v; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var column = columnDataB[j];
Complex sum = Complex.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * column[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
});
}
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>

134
src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Complex32.cs
View File

@ -470,12 +470,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
@ -556,12 +550,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@ -940,128 +928,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA,
Complex32[] b,
int rowsB, int columnsB, Complex32 beta, Complex32[] c)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (c == null)
{
throw new ArgumentNullException("c");
}
if (transposeA != Transpose.DontTranspose)
{
var swap = rowsA;
rowsA = columnsA;
columnsA = swap;
}
if (transposeB != Transpose.DontTranspose)
{
var swap = rowsB;
rowsB = columnsB;
columnsB = swap;
}
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsA, rowsB));
}
if (rowsA * columnsA != a.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsA, columnsA, a.Length));
}
if (rowsB * columnsB != b.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsB, columnsB, b.Length));
}
if (rowsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsA, columnsB, c.Length));
}
// handle degenerate cases
if (beta == Complex32.Zero)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != Complex32.One)
{
ScaleArray(beta, c, c);
}
if (alpha == Complex32.Zero)
{
return;
}
// Extract column arrays
var columnDataB = new Complex32[columnsB][];
for (int i = 0; i < columnDataB.Length; i++)
{
var column = new Complex32[rowsB];
GetColumn(transposeB, i, rowsB, columnsB, b, column);
columnDataB[i] = column;
}
var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
var row = new Complex32[columnsA];
for (int i = 0; i < rowsA; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var col = columnDataB[j];
Complex32 sum = Complex32.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * col[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
}
else
{
CommonParallel.For(0, rowsA, 1, (u, v) =>
{
var row = new Complex32[columnsA];
for (int i = u; i < v; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var column = columnDataB[j];
Complex32 sum = Complex32.Zero;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * column[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
});
}
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>

134
src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Double.cs
View File

@ -465,12 +465,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
@ -551,12 +545,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@ -826,128 +814,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA,
double[] b,
int rowsB, int columnsB, double beta, double[] c)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (c == null)
{
throw new ArgumentNullException("c");
}
if (transposeA != Transpose.DontTranspose)
{
var swap = rowsA;
rowsA = columnsA;
columnsA = swap;
}
if (transposeB != Transpose.DontTranspose)
{
var swap = rowsB;
rowsB = columnsB;
columnsB = swap;
}
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsA, rowsB));
}
if (rowsA * columnsA != a.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsA, columnsA, a.Length));
}
if (rowsB * columnsB != b.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsB, columnsB, b.Length));
}
if (rowsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsA, columnsB, c.Length));
}
// handle degenerate cases
if (beta == 0.0)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != 1.0)
{
ScaleArray(beta, c, c);
}
if (alpha == 0.0)
{
return;
}
// Extract column arrays
var columnDataB = new double[columnsB][];
for (int i = 0; i < columnDataB.Length; i++)
{
var column = new double[rowsB];
GetColumn(transposeB, i, rowsB, columnsB, b, column);
columnDataB[i] = column;
}
var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
var row = new double[columnsA];
for (int i = 0; i < rowsA; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var col = columnDataB[j];
double sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * col[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
}
else
{
CommonParallel.For(0, rowsA, 1, (u, v) =>
{
var row = new double[columnsA];
for (int i = u; i < v; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var column = columnDataB[j];
double sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * column[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
});
}
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>

134
src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.Single.cs
View File

@ -465,12 +465,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
public virtual void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result);
return;
}
// First check some basic requirement on the parameters of the matrix multiplication.
if (x == null)
{
@ -551,12 +545,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// <param name="c">The c matrix.</param>
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c)
{
if (_variation == Variation.Experimental)
{
MatrixMultiplyWithUpdateExperimental(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
return;
}
int m; // The number of rows of matrix op(A) and of the matrix C.
int n; // The number of columns of matrix op(B) and of the matrix C.
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
@ -826,128 +814,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
}
}
public void MatrixMultiplyWithUpdateExperimental(
Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA,
float[] b,
int rowsB, int columnsB, float beta, float[] c)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (c == null)
{
throw new ArgumentNullException("c");
}
if (transposeA != Transpose.DontTranspose)
{
var swap = rowsA;
rowsA = columnsA;
columnsA = swap;
}
if (transposeB != Transpose.DontTranspose)
{
var swap = rowsB;
rowsB = columnsB;
columnsB = swap;
}
if (columnsA != rowsB)
{
throw new ArgumentOutOfRangeException(string.Format("columnsA ({0}) != rowsB ({1})", columnsA, rowsB));
}
if (rowsA * columnsA != a.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsA ({1}) != a.Length ({2})", rowsA, columnsA, a.Length));
}
if (rowsB * columnsB != b.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsB ({0}) * columnsB ({1}) != b.Length ({2})", rowsB, columnsB, b.Length));
}
if (rowsA * columnsB != c.Length)
{
throw new ArgumentOutOfRangeException(string.Format("rowsA ({0}) * columnsB ({1}) != c.Length ({2})", rowsA, columnsB, c.Length));
}
// handle degenerate cases
if (beta == 0.0)
{
Array.Clear(c, 0, c.Length);
}
else if (beta != 1.0)
{
ScaleArray(beta, c, c);
}
if (alpha == 0.0)
{
return;
}
// Extract column arrays
var columnDataB = new float[columnsB][];
for (int i = 0; i < columnDataB.Length; i++)
{
var column = new float[rowsB];
GetColumn(transposeB, i, rowsB, columnsB, b, column);
columnDataB[i] = column;
}
var shouldNotParallelize = rowsA + columnsB + columnsA < Control.ParallelizeOrder || Control.MaxDegreeOfParallelism < 2;
if (shouldNotParallelize)
{
var row = new float[columnsA];
for (int i = 0; i < rowsA; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var col = columnDataB[j];
float sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * col[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
}
else
{
CommonParallel.For(0, rowsA, 1, (u, v) =>
{
var row = new float[columnsA];
for (int i = u; i < v; i++)
{
GetRow(transposeA, i, rowsA, columnsA, a, row);
for (int j = 0; j < columnsB; j++)
{
var column = columnDataB[j];
float sum = 0;
for (int ii = 0; ii < row.Length; ii++)
{
sum += row[ii] * column[ii];
}
c[j * rowsA + i] += alpha * sum;
}
}
});
}
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>

12
src/Numerics/Providers/LinearAlgebra/ManagedReference/ManagedReferenceLinearAlgebraProvider.cs
View File

@ -42,18 +42,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.ManagedReference
/// </summary>
internal partial class ManagedReferenceLinearAlgebraProvider : ILinearAlgebraProvider
{
private readonly Variation _variation;
internal ManagedReferenceLinearAlgebraProvider()
{
_variation = Variation.Experimental;
}
internal ManagedReferenceLinearAlgebraProvider(Variation variation)
{
_variation = variation;
}
/// <summary>
/// Try to find out whether the provider is available, at least in principle.
/// Verification may still fail if available, but it will certainly fail if unavailable.

Write Preview
Loading…
Cancel
Save