diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
index b10d4e0c..4bde927a 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
@@ -176,11 +176,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// Multiples two matrices. result = x * y
///
/// The x matrix.
+ /// The number of rows in the x matrix.
+ /// The number of columns in the x matrix.
/// The y matrix.
+ /// The number of rows in the y matrix.
+ /// The number of columns in the y matrix.
/// Where to store the result of the multiplication.
/// This is a simplified version of the BLAS GEMM routine with alpha
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.
- void MatrixMultiply(T[] x, T[] y, T[] result);
+ void MatrixMultiply(T[] x, int xRows, int xColumns, T[] y, int yRows, int yColumns, T[] result);
///
/// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c
diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
index 6a528b5a..6133df16 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
@@ -104,11 +104,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Computes the dot product between two vectors.
+ /// Computes the dot product of x and y.
///
- /// The first argument of the dot product.
- /// The second argument of the dot product.
- /// The dot product between and .
+ /// The vector x.
+ /// The vector y.
+ /// The dot product of x and y.
+ /// This is equivalent to the DOT BLAS routine.
public double DotProduct(double[] x, double[] y)
{
if (y == null)
@@ -137,11 +138,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Adds two arrays together and writes the result in a third array.
+ /// Does a point wise add of two arrays z = x + y. This can be used
+ /// to add vectors or matrices.
///
- /// The first argument to add.
- /// The second argument to add.
- /// The result to write the addition into.
+ /// The array x.
+ /// The array y.
+ /// The result of the addition.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void AddArrays(double[] x, double[] y, double[] result)
{
if (y == null)
@@ -168,11 +173,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Subtract two arrays and writes the result in a third array.
+ /// Does a point wise subtraction of two arrays z = x - y. This can be used
+ /// to subtract vectors or matrices.
///
- /// The first argument to subtract.
- /// The second argument to subtract.
- /// The result to write the subtraction into.
+ /// The array x.
+ /// The array y.
+ /// The result of the subtraction.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void SubtractArrays(double[] x, double[] y, double[] result)
{
if (y == null)
@@ -199,11 +208,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Pointwise multiplies two arrays and writes the result in a third array.
+ /// Does a point wise multiplication of two arrays z = x * y. This can be used
+ /// to multiple elements of vectors or matrices.
///
- /// The first argument to pointwise multiply.
- /// The second argument to pointwise multiply.
- /// The result to write the pointwise multiplication into.
+ /// The array x.
+ /// The array y.
+ /// The result of the point wise multiplication.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result)
{
if (y == null)
@@ -239,9 +252,92 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new NotImplementedException();
}
- public void MatrixMultiply(double[] x, double[] y, double[] result)
- {
- throw new NotImplementedException();
+ ///
+ /// Multiples two matrices. result = x * y
+ ///
+ /// The x matrix.
+ /// The number of rows in the x matrix.
+ /// The number of columns in the x matrix.
+ /// The y matrix.
+ /// The number of rows in the y matrix.
+ /// The number of columns in the y matrix.
+ /// Where to store the result of the multiplication.
+ /// This is a simplified version of the BLAS GEMM routine with alpha
+ /// set to 1.0 and beta set to 0.0, and x and y are not transposed.
+ public void MatrixMultiply(double[] x, int xRows, int xColumns, double[] y, int yRows, int yColumns, double[] result)
+ {
+ // First check some basic requirement on the parameters of the matrix multiplication.
+ if (x == null)
+ {
+ throw new ArgumentNullException("x");
+ }
+ if (y == null)
+ {
+ throw new ArgumentNullException("y");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (xRows * xColumns != x.Length)
+ {
+ throw new ArgumentException("x.Length != xRows * xColumns");
+ }
+
+ if (yRows * yColumns != y.Length)
+ {
+ throw new ArgumentException("y.Length != yRows * yColumns");
+ }
+
+ if (xColumns != yRows)
+ {
+ throw new ArgumentException("xColumns != yRows");
+ }
+
+ if (xRows * yColumns != result.Length)
+ {
+ throw new ArgumentException("xRows * yColumns != result.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
+ {
+ xdata = x;
+ }
+
+ double[] ydata;
+ if (ReferenceEquals(y, result))
+ {
+ ydata = (double[]) y.Clone();
+ }
+ else
+ {
+ ydata = y;
+ }
+
+ // Start the actual matrix multiplication.
+ // TODO - For small matrices we should get rid of the parallelism because of startup costs.
+ // Perhaps the following implementations would be a good one
+ // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
+ Parallel.For(0, xRows, i =>
+ {
+ for (int j = 0; j < yColumns; j++)
+ {
+ for (int k = 0; k < xColumns; k++)
+ {
+ result[j + yColumns * i] += xdata[k + xColumns * i] * ydata[j + yColumns * k];
+ }
+ }
+ });
}
public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c)
@@ -424,11 +520,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Computes the dot product between two vectors.
+ /// Computes the dot product of x and y.
///
- /// The first argument of the dot product.
- /// The second argument of the dot product.
- /// The dot product between and .
+ /// The vector x.
+ /// The vector y.
+ /// The dot product of x and y.
+ /// This is equivalent to the DOT BLAS routine.
public float DotProduct(float[] x, float[] y)
{
if (y == null)
@@ -457,11 +554,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Adds two arrays together and writes the result in a third array.
+ /// Does a point wise add of two arrays z = x + y. This can be used
+ /// to add vectors or matrices.
///
- /// The first argument to add.
- /// The second argument to add.
- /// The result to write the addition into.
+ /// The array x.
+ /// The array y.
+ /// The result of the addition.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void AddArrays(float[] x, float[] y, float[] result)
{
if (y == null)
@@ -488,11 +589,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Subtract two arrays and writes the result in a third array.
+ /// Does a point wise subtraction of two arrays z = x - y. This can be used
+ /// to subtract vectors or matrices.
///
- /// The first argument to subtract.
- /// The second argument to subtract.
- /// The result to write the subtraction into.
+ /// The array x.
+ /// The array y.
+ /// The result of the subtraction.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void SubtractArrays(float[] x, float[] y, float[] result)
{
if (y == null)
@@ -519,11 +624,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Pointwise multiplies two arrays and writes the result in a third array.
+ /// Does a point wise multiplication of two arrays z = x * y. This can be used
+ /// to multiple elements of vectors or matrices.
///
- /// The first argument to pointwise multiply.
- /// The second argument to pointwise multiply.
- /// The result to write the pointwise multiplication into.
+ /// The array x.
+ /// The array y.
+ /// The result of the point wise multiplication.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result)
{
if (y == null)
@@ -559,9 +668,92 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new NotImplementedException();
}
- public void MatrixMultiply(float[] x, float[] y, float[] result)
- {
- throw new NotImplementedException();
+ ///
+ /// Multiples two matrices. result = x * y
+ ///
+ /// The x matrix.
+ /// The number of rows in the x matrix.
+ /// The number of columns in the x matrix.
+ /// The y matrix.
+ /// The number of rows in the y matrix.
+ /// The number of columns in the y matrix.
+ /// Where to store the result of the multiplication.
+ /// This is a simplified version of the BLAS GEMM routine with alpha
+ /// set to 1.0 and beta set to 0.0, and x and y are not transposed.
+ public void MatrixMultiply(float[] x, int xRows, int xColumns, float[] y, int yRows, int yColumns, float[] result)
+ {
+ // First check some basic requirement on the parameters of the matrix multiplication.
+ if (x == null)
+ {
+ throw new ArgumentNullException("x");
+ }
+ if (y == null)
+ {
+ throw new ArgumentNullException("y");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (xRows * xColumns != x.Length)
+ {
+ throw new ArgumentException("x.Length != xRows * xColumns");
+ }
+
+ if (yRows * yColumns != y.Length)
+ {
+ throw new ArgumentException("y.Length != yRows * yColumns");
+ }
+
+ if (xColumns != yRows)
+ {
+ throw new ArgumentException("xColumns != yRows");
+ }
+
+ if (xRows * yColumns != result.Length)
+ {
+ throw new ArgumentException("xRows * yColumns != result.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
+ {
+ xdata = x;
+ }
+
+ float[] ydata;
+ if (ReferenceEquals(y, result))
+ {
+ ydata = (float[]) y.Clone();
+ }
+ else
+ {
+ ydata = y;
+ }
+
+ // Start the actual matrix multiplication.
+ // TODO - For small matrices we should get rid of the parallelism because of startup costs.
+ // Perhaps the following implementations would be a good one
+ // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
+ Parallel.For(0, xRows, i =>
+ {
+ for (int j = 0; j < yColumns; j++)
+ {
+ for (int k = 0; k < xColumns; k++)
+ {
+ result[j + yColumns * i] += xdata[k + xColumns * i] * ydata[j + yColumns * k];
+ }
+ }
+ });
}
public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, float[] b, float beta, float[] c)
@@ -744,11 +936,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Computes the dot product between two vectors.
+ /// Computes the dot product of x and y.
///
- /// The first argument of the dot product.
- /// The second argument of the dot product.
- /// The dot product between and .
+ /// The vector x.
+ /// The vector y.
+ /// The dot product of x and y.
+ /// This is equivalent to the DOT BLAS routine.
public Complex DotProduct(Complex[] x, Complex[] y)
{
if (y == null)
@@ -777,11 +970,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Adds two arrays together and writes the result in a third array.
+ /// Does a point wise add of two arrays z = x + y. This can be used
+ /// to add vectors or matrices.
///
- /// The first argument to add.
- /// The second argument to add.
- /// The result to write the addition into.
+ /// The array x.
+ /// The array y.
+ /// The result of the addition.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void AddArrays(Complex[] x, Complex[] y, Complex[] result)
{
if (y == null)
@@ -808,11 +1005,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Subtract two arrays and writes the result in a third array.
+ /// Does a point wise subtraction of two arrays z = x - y. This can be used
+ /// to subtract vectors or matrices.
///
- /// The first argument to subtract.
- /// The second argument to subtract.
- /// The result to write the subtraction into.
+ /// The array x.
+ /// The array y.
+ /// The result of the subtraction.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result)
{
if (y == null)
@@ -839,11 +1040,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Pointwise multiplies two arrays and writes the result in a third array.
+ /// Does a point wise multiplication of two arrays z = x * y. This can be used
+ /// to multiple elements of vectors or matrices.
///
- /// The first argument to pointwise multiply.
- /// The second argument to pointwise multiply.
- /// The result to write the pointwise multiplication into.
+ /// The array x.
+ /// The array y.
+ /// The result of the point wise multiplication.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result)
{
if (y == null)
@@ -878,10 +1083,93 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
{
throw new NotImplementedException();
}
+
+ ///
+ /// Multiples two matrices. result = x * y
+ ///
+ /// The x matrix.
+ /// The number of rows in the x matrix.
+ /// The number of columns in the x matrix.
+ /// The y matrix.
+ /// The number of rows in the y matrix.
+ /// The number of columns in the y matrix.
+ /// Where to store the result of the multiplication.
+ /// This is a simplified version of the BLAS GEMM routine with alpha
+ /// set to 1.0 and beta set to 0.0, and x and y are not transposed.
+ public void MatrixMultiply(Complex[] x, int xRows, int xColumns, Complex[] y, int yRows, int yColumns, Complex[] result)
+ {
+ // First check some basic requirement on the parameters of the matrix multiplication.
+ if (x == null)
+ {
+ throw new ArgumentNullException("x");
+ }
+ if (y == null)
+ {
+ throw new ArgumentNullException("y");
+ }
- public void MatrixMultiply(Complex[] x, Complex[] y, Complex[] result)
- {
- throw new NotImplementedException();
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (xRows * xColumns != x.Length)
+ {
+ throw new ArgumentException("x.Length != xRows * xColumns");
+ }
+
+ if (yRows * yColumns != y.Length)
+ {
+ throw new ArgumentException("y.Length != yRows * yColumns");
+ }
+
+ if (xColumns != yRows)
+ {
+ throw new ArgumentException("xColumns != yRows");
+ }
+
+ if (xRows * yColumns != result.Length)
+ {
+ throw new ArgumentException("xRows * yColumns != result.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 = (Complex[]) x.Clone();
+ }
+ else
+ {
+ xdata = x;
+ }
+
+ double[] ydata;
+ if (ReferenceEquals(y, result))
+ {
+ ydata = (Complex[]) y.Clone();
+ }
+ else
+ {
+ ydata = y;
+ }
+
+ // Start the actual matrix multiplication.
+ // TODO - For small matrices we should get rid of the parallelism because of startup costs.
+ // Perhaps the following implementations would be a good one
+ // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
+ Parallel.For(0, xRows, i =>
+ {
+ for (int j = 0; j < yColumns; j++)
+ {
+ for (int k = 0; k < xColumns; k++)
+ {
+ result[j + yColumns * i] += xdata[k + xColumns * i] * ydata[j + yColumns * k];
+ }
+ }
+ });
}
public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, Complex[] b, Complex beta, Complex[] c)
@@ -1064,11 +1352,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Computes the dot product between two vectors.
+ /// Computes the dot product of x and y.
///
- /// The first argument of the dot product.
- /// The second argument of the dot product.
- /// The dot product between and .
+ /// The vector x.
+ /// The vector y.
+ /// The dot product of x and y.
+ /// This is equivalent to the DOT BLAS routine.
public Complex32 DotProduct(Complex32[] x, Complex32[] y)
{
if (y == null)
@@ -1097,11 +1386,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Adds two arrays together and writes the result in a third array.
+ /// Does a point wise add of two arrays z = x + y. This can be used
+ /// to add vectors or matrices.
///
- /// The first argument to add.
- /// The second argument to add.
- /// The result to write the addition into.
+ /// The array x.
+ /// The array y.
+ /// The result of the addition.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
@@ -1128,11 +1421,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Subtract two arrays and writes the result in a third array.
+ /// Does a point wise subtraction of two arrays z = x - y. This can be used
+ /// to subtract vectors or matrices.
///
- /// The first argument to subtract.
- /// The second argument to subtract.
- /// The result to write the subtraction into.
+ /// The array x.
+ /// The array y.
+ /// The result of the subtraction.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
@@ -1159,11 +1456,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
///
- /// Pointwise multiplies two arrays and writes the result in a third array.
+ /// Does a point wise multiplication of two arrays z = x * y. This can be used
+ /// to multiple elements of vectors or matrices.
///
- /// The first argument to pointwise multiply.
- /// The second argument to pointwise multiply.
- /// The result to write the pointwise multiplication into.
+ /// The array x.
+ /// The array y.
+ /// The result of the point wise multiplication.
+ /// There is no equivalent BLAS routine, but many libraries
+ /// provide optimized (parallel and/or vectorized) versions of this
+ /// routine.
public void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
@@ -1199,9 +1500,92 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new NotImplementedException();
}
- public void MatrixMultiply(Complex32[] x, Complex32[] y, Complex32[] result)
- {
- throw new NotImplementedException();
+ ///
+ /// Multiples two matrices. result = x * y
+ ///
+ /// The x matrix.
+ /// The number of rows in the x matrix.
+ /// The number of columns in the x matrix.
+ /// The y matrix.
+ /// The number of rows in the y matrix.
+ /// The number of columns in the y matrix.
+ /// Where to store the result of the multiplication.
+ /// This is a simplified version of the BLAS GEMM routine with alpha
+ /// set to 1.0 and beta set to 0.0, and x and y are not transposed.
+ public void MatrixMultiply(Complex32[] x, int xRows, int xColumns, Complex32[] y, int yRows, int yColumns, Complex32[] result)
+ {
+ // First check some basic requirement on the parameters of the matrix multiplication.
+ if (x == null)
+ {
+ throw new ArgumentNullException("x");
+ }
+ if (y == null)
+ {
+ throw new ArgumentNullException("y");
+ }
+
+ if (result == null)
+ {
+ throw new ArgumentNullException("result");
+ }
+
+ if (xRows * xColumns != x.Length)
+ {
+ throw new ArgumentException("x.Length != xRows * xColumns");
+ }
+
+ if (yRows * yColumns != y.Length)
+ {
+ throw new ArgumentException("y.Length != yRows * yColumns");
+ }
+
+ if (xColumns != yRows)
+ {
+ throw new ArgumentException("xColumns != yRows");
+ }
+
+ if (xRows * yColumns != result.Length)
+ {
+ throw new ArgumentException("xRows * yColumns != result.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))
+ {
+ xdata = (Complex32[]) x.Clone();
+ }
+ else
+ {
+ xdata = x;
+ }
+
+ Complex32[] ydata;
+ if (ReferenceEquals(y, result))
+ {
+ ydata = (Complex32[]) y.Clone();
+ }
+ else
+ {
+ ydata = y;
+ }
+
+ // Start the actual matrix multiplication.
+ // TODO - For small matrices we should get rid of the parallelism because of startup costs.
+ // Perhaps the following implementations would be a good one
+ // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/
+ Parallel.For(0, xRows, i =>
+ {
+ for (int j = 0; j < yColumns; j++)
+ {
+ for (int k = 0; k < xColumns; k++)
+ {
+ result[j + yColumns * i] += xdata[k + xColumns * i] * ydata[j + yColumns * k];
+ }
+ }
+ });
}
public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, Complex32[] b, Complex32 beta, Complex32[] c)