From 4f694b0cf2ba2b36f3ba4a3ea52952f6853b826c Mon Sep 17 00:00:00 2001 From: Jurgen Van Gael Date: Sat, 28 Nov 2009 17:39:20 +0800 Subject: [PATCH] Ported MatrixMultiplyWithUpdate implementation Ported CholeskyFactor implementation for double and float. --- .../Atlas/AtlasLinearAlgebraProvider.cs | 18 +- .../ILinearAlgebraProviderOfT.cs | 3 +- .../ManagedLinearAlgebraProvider.cs | 1404 ++++++++++++++++- .../Mkl/MklLinearAlgebraProvider.cs | 18 +- .../NativeAlgebraProvider.include | 12 +- src/Numerics/Constants.cs | 5 + src/Numerics/Properties/Resources.Designer.cs | 13 +- src/Numerics/Properties/Resources.resx | 6 + 8 files changed, 1444 insertions(+), 35 deletions(-) diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs index 4ede1bc8..ff6cd573 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs @@ -24,11 +24,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. -<<<<<<< HEAD:src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs - Last generated on: 14/11/2009 20:09:22 -======= - Last generated on: 11/13/2009 9:49:35 AM ->>>>>>> c1c4de3... adding dot product:src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs + Last generated on: 28/11/2009 09:32:52 */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { @@ -350,8 +346,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(double[] a) + public void CholeskyFactor(double[] a, int order) { throw new NotImplementedException(); } @@ -846,8 +843,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(float[] a) + public void CholeskyFactor(float[] a, int order) { throw new NotImplementedException(); } @@ -1342,8 +1340,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex[] a) + public void CholeskyFactor(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1838,8 +1837,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex32[] a) + public void CholeskyFactor(Complex32[] a, int order) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 6327e3de..6adad9a3 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -294,8 +294,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - void CholeskyFactor(T[] a); + void CholeskyFactor(T[] a, int order); /// /// Solves A*X=B for X using Cholesky factorization. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 81bbf501..6bf604ec 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -357,7 +357,330 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int aRows, int aColumns, double[] b, int bRows, int bColumns, double beta, double[] c) { - throw new NotImplementedException(); + // Choose nonsensical values for the number of rows and columns in c; fill them in depending + // on the operations on a and b. + int cRows = -1; + int cColumns = -1; + + // 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 (aRows != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bRows; + } + else if ((int)transposeA > 111) + { + if (aRows != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bColumns; + } + else if ((int)transposeB > 111) + { + if (aColumns != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bRows; + } + else + { + if (aColumns != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bColumns; + } + + if (Precision.AlmostEqual(0.0, alpha) + && Precision.AlmostEqual(0.0, beta)) + { + 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 (Precision.AlmostEqual(1.0, alpha)) + { + if (Precision.AlmostEqual(0.0, beta)) + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + double s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + double s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + } } public void LUFactor(double[] a, int[] ipiv) @@ -405,9 +728,48 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void CholeskyFactor(double[] a) + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + public void CholeskyFactor(double[] a, int order) { - throw new NotImplementedException(); + double[] factor = new double[a.Length]; + + for (int j = 0; j < order; j++) + { + double d = 0.0; + int index; + for (int k = 0; k < j; k++) + { + double s = 0.0; + int i; + for (i = 0; i < k; i++) + { + s += factor[i * order + k] * factor[i * order + j]; + } + int tmp = k * order; + index = tmp + j; + factor[index] = s = (a[index] - s) / factor[tmp + k]; + d += s * s; + } + index = j * order + j; + d = a[index] - d; + if (d <= 0.0) + { + throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); + } + factor[index] = System.Math.Sqrt(d); + for (int k = j + 1; k < order; k++) + { + factor[k * order + j] = 0.0; + } + } + + Buffer.BlockCopy(factor, 0, a, 0, factor.Length * Constants.SizeOfDouble); } public void CholeskySolve(int columnsOfB, double[] a, double[] b) @@ -788,7 +1150,330 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int aRows, int aColumns, float[] b, int bRows, int bColumns, float beta, float[] c) { - throw new NotImplementedException(); + // Choose nonsensical values for the number of rows and columns in c; fill them in depending + // on the operations on a and b. + int cRows = -1; + int cColumns = -1; + + // 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 (aRows != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bRows; + } + else if ((int)transposeA > 111) + { + if (aRows != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bColumns; + } + else if ((int)transposeB > 111) + { + if (aColumns != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bRows; + } + else + { + if (aColumns != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bColumns; + } + + if (Precision.AlmostEqual(0.0, alpha) + && Precision.AlmostEqual(0.0, beta)) + { + 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 (Precision.AlmostEqual(1.0, alpha)) + { + if (Precision.AlmostEqual(0.0, beta)) + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + float s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + float s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + } } public void LUFactor(float[] a, int[] ipiv) @@ -836,9 +1521,48 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void CholeskyFactor(float[] a) + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + public void CholeskyFactor(float[] a, int order) { - throw new NotImplementedException(); + float[] factor = new float[a.Length]; + + for (int j = 0; j < order; j++) + { + float d = 0.0F; + int index; + for (int k = 0; k < j; k++) + { + float s = 0.0F; + int i; + for (i = 0; i < k; i++) + { + s += factor[i * order + k] * factor[i * order + j]; + } + int tmp = k * order; + index = tmp + j; + factor[index] = s = (a[index] - s) / factor[tmp + k]; + d += s * s; + } + index = j * order + j; + d = a[index] - d; + if (d <= 0.0F) + { + throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); + } + factor[index] = (float) System.Math.Sqrt(d); + for (int k = j + 1; k < order; k++) + { + factor[k * order + j] = 0.0F; + } + } + + Buffer.BlockCopy(factor, 0, a, 0, factor.Length * Constants.SizeOfFloat); } public void CholeskySolve(int columnsOfB, float[] a, float[] b) @@ -1219,7 +1943,330 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int aRows, int aColumns, Complex[] b, int bRows, int bColumns, Complex beta, Complex[] c) { - throw new NotImplementedException(); + // Choose nonsensical values for the number of rows and columns in c; fill them in depending + // on the operations on a and b. + int cRows = -1; + int cColumns = -1; + + // 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 (aRows != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bRows; + } + else if ((int)transposeA > 111) + { + if (aRows != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bColumns; + } + else if ((int)transposeB > 111) + { + if (aColumns != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bRows; + } + else + { + if (aColumns != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bColumns; + } + + if (Precision.AlmostEqual(0.0, alpha) + && Precision.AlmostEqual(0.0, beta)) + { + 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 (Precision.AlmostEqual(1.0, alpha)) + { + if (Precision.AlmostEqual(0.0, beta)) + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + } } public void LUFactor(Complex[] a, int[] ipiv) @@ -1267,7 +2314,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void CholeskyFactor(Complex[] a) + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + public void CholeskyFactor(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1650,7 +2704,330 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int aRows, int aColumns, Complex32[] b, int bRows, int bColumns, Complex32 beta, Complex32[] c) { - throw new NotImplementedException(); + // Choose nonsensical values for the number of rows and columns in c; fill them in depending + // on the operations on a and b. + int cRows = -1; + int cColumns = -1; + + // 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 (aRows != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bRows; + } + else if ((int)transposeA > 111) + { + if (aRows != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aColumns * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aColumns; + cColumns = bColumns; + } + else if ((int)transposeB > 111) + { + if (aColumns != bColumns) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bRows != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bRows; + } + else + { + if (aColumns != bRows) + { + throw new ArgumentOutOfRangeException(); + } + + if (aRows * bColumns != c.Length) + { + throw new ArgumentOutOfRangeException(); + } + + cRows = aRows; + cColumns = bColumns; + } + + if (Precision.AlmostEqual((Complex32)0.0, alpha) + && Precision.AlmostEqual((Complex32)0.0, beta)) + { + 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 (Precision.AlmostEqual((Complex32)1.0, alpha)) + { + if (Precision.AlmostEqual((Complex32)0.0, beta)) + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s; + } + }); + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); + } + } + } + else + { + if ((int)transposeA > 111 && (int)transposeB > 111) + { + Parallel.For(0, aColumns, j => + { + int jIndex = j * cRows; + for (int i = 0; i != bRows; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); + } + else if ((int)transposeA > 111) + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aColumns; i++) + { + int iIndex = i * aRows; + Complex32 s = 0; + for (int l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + else if ((int)transposeB > 111) + { + Parallel.For(0, bRows, j => + { + int jIndex = j * cRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); + } + else + { + Parallel.For(0, bColumns, j => + { + int jcIndex = j * cRows; + int jbIndex = j * bRows; + for (int i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (int l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); + } + } } public void LUFactor(Complex32[] a, int[] ipiv) @@ -1698,7 +3075,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void CholeskyFactor(Complex32[] a) + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + public void CholeskyFactor(Complex32[] a, int order) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 3d8ce29a..9150aaad 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -24,11 +24,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. -<<<<<<< HEAD:src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs - Last generated on: 14/11/2009 20:09:20 -======= - Last generated on: 11/6/2009 3:03:50 PM ->>>>>>> c1c4de3... adding dot product:src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs + Last generated on: 28/11/2009 09:32:55 */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { @@ -350,8 +346,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(double[] a) + public void CholeskyFactor(double[] a, int order) { throw new NotImplementedException(); } @@ -846,8 +843,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(float[] a) + public void CholeskyFactor(float[] a, int order) { throw new NotImplementedException(); } @@ -1342,8 +1340,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex[] a) + public void CholeskyFactor(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1838,8 +1837,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex32[] a) + public void CholeskyFactor(Complex32[] a, int order) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include index a916be82..60be35f5 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include +++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include @@ -346,8 +346,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(double[] a) + public void CholeskyFactor(double[] a, int order) { throw new NotImplementedException(); } @@ -842,8 +843,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(float[] a) + public void CholeskyFactor(float[] a, int order) { throw new NotImplementedException(); } @@ -1338,8 +1340,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex[] a) + public void CholeskyFactor(Complex[] a, int order) { throw new NotImplementedException(); } @@ -1834,8 +1837,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the /// the Cholesky factorization. + /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex32[] a) + public void CholeskyFactor(Complex32[] a, int order) { throw new NotImplementedException(); } diff --git a/src/Numerics/Constants.cs b/src/Numerics/Constants.cs index 67718368..fce41ded 100644 --- a/src/Numerics/Constants.cs +++ b/src/Numerics/Constants.cs @@ -162,6 +162,11 @@ namespace MathNet.Numerics /// public const int SizeOfDouble = sizeof(double); + /// + /// The size of a float in bytes. + /// + public const int SizeOfFloat = sizeof(float); + #endregion #region UNIVERSAL CONSTANTS diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 6d990350..cc006619 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // // This code was generated by a tool. -// Runtime Version:2.0.50727.4927 +// Runtime Version:2.0.50727.4200 // // Changes to this file may cause incorrect behavior and will be lost if // the code is regenerated. @@ -61,7 +61,7 @@ namespace MathNet.Numerics.Properties { } /// - /// Looks up a localized string similar to All arrays must have the same length.. + /// Looks up a localized string similar to The array arguments must have the same length.. /// internal static string ArgumentArraysSameLength { get { @@ -168,6 +168,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to Matrix must be positive definite.. + /// + internal static string ArgumentMatrixPositiveDefinite { + get { + return ResourceManager.GetString("ArgumentMatrixPositiveDefinite", resourceCulture); + } + } + /// /// Looks up a localized string similar to Matrix column dimensions must agree.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index f8f60969..88081ee6 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -288,4 +288,10 @@ The sampler's proposal distribution is not upper bounding the target density. + + Matrix must be positive definite. + + + The array arguments must have the same length. + \ No newline at end of file