From d27833fd5e2eb151ffe2ea352133a483d60d4906 Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Fri, 13 Nov 2009 04:50:01 +0800 Subject: [PATCH] Added a matrix multiplication algorithm. --- .../ILinearAlgebraProviderOfT.cs | 6 +- .../ManagedLinearAlgebraProvider.cs | 536 +++++++++++++++--- 2 files changed, 465 insertions(+), 77 deletions(-) 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)