diff --git a/build/build.proj b/build/build.proj index 86eddcec..cb636df4 100644 --- a/build/build.proj +++ b/build/build.proj @@ -1,4 +1,4 @@ - + diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj index 987d4026..f0dd7450 100644 --- a/src/FSharp/FSharp.fsproj +++ b/src/FSharp/FSharp.fsproj @@ -33,8 +33,8 @@ MathNet.Numerics.FSharp.XML + - 3.5 diff --git a/src/FSharpExamples/FSharpExamples.fsproj b/src/FSharpExamples/FSharpExamples.fsproj index 0d7941a9..8076204c 100644 --- a/src/FSharpExamples/FSharpExamples.fsproj +++ b/src/FSharpExamples/FSharpExamples.fsproj @@ -31,18 +31,12 @@ - 3.5 - - FSharp - {37e8e802-a354-4114-bfc1-6e1357da605b} - True - Numerics {b7cae5f4-a23f-4438-b5be-41226618b695} diff --git a/src/FSharpUnitTests/FSharpUnitTests.fsproj b/src/FSharpUnitTests/FSharpUnitTests.fsproj index 3a68e243..eab9e601 100644 --- a/src/FSharpUnitTests/FSharpUnitTests.fsproj +++ b/src/FSharpUnitTests/FSharpUnitTests.fsproj @@ -31,7 +31,6 @@ - 3.5 @@ -42,11 +41,6 @@ - - FSharp - {37e8e802-a354-4114-bfc1-6e1357da605b} - True - Numerics {b7cae5f4-a23f-4438-b5be-41226618b695} diff --git a/src/NativeWrappers/Common/blas.c b/src/NativeWrappers/Common/blas.c index b8314040..891ce5c7 100644 --- a/src/NativeWrappers/Common/blas.c +++ b/src/NativeWrappers/Common/blas.c @@ -28,7 +28,7 @@ DLLEXPORT void d_scale(const int n, const double alpha, double x[]){ DLLEXPORT void c_scale(const int n, const Complex8 alpha, Complex8 x[]){ cblas_cscal(n, &alpha, x, 1); } - + DLLEXPORT void z_scale(const int n, const Complex16 alpha, Complex16 x[]){ cblas_zscal(n, &alpha, x, 1); } @@ -53,6 +53,30 @@ DLLEXPORT Complex16 z_dot_product(const int n, const Complex16 x[], const Comple return ret; } +DLLEXPORT void s_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const float alpha, const float x[], const float y[], const float beta, float c[]){ + int lda = transA == CblasNoTrans ? m : k; + int ldb = transB == CblasNoTrans ? k : n; + + cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} +DLLEXPORT void d_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const double alpha, const double x[], const double y[], const double beta, double c[]){ + int lda = transA == CblasNoTrans ? m : k; + int ldb = transB == CblasNoTrans ? k : n; + + cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} +DLLEXPORT void c_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const Complex8 alpha, const Complex8 x[], const Complex8 y[], const Complex8 beta, Complex8 c[]){ + int lda = transA == CblasNoTrans ? m : k; + int ldb = transB == CblasNoTrans ? k : n; + cblas_cgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); +} + +DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const int m, const int n, const int k, const Complex16 alpha, const Complex16 x[], const Complex16 y[], const Complex16 beta, Complex16 c[]){ + int lda = transA == CblasNoTrans ? m : k; + int ldb = transB == CblasNoTrans ? k : n; + + cblas_zgemm(CblasColMajor, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); +} diff --git a/src/NativeWrappers/NativeWrappers.sln b/src/NativeWrappers/NativeWrappers.sln index 24222cd2..676bbbe0 100644 --- a/src/NativeWrappers/NativeWrappers.sln +++ b/src/NativeWrappers/NativeWrappers.sln @@ -56,8 +56,8 @@ Global {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Mixed Platforms.Build.0 = Release|x64 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Win32.ActiveCfg = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Win32.Build.0 = Release|Win32 - {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.ActiveCfg = Release|x64 - {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.Build.0 = Release|x64 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.Build.0 = Release|Win32 {0EFC01B9-1F75-4BFD-ADB6-3FF18B2B9B5E}.Debug|Any CPU.ActiveCfg = Debug|Any CPU {0EFC01B9-1F75-4BFD-ADB6-3FF18B2B9B5E}.Debug|Any CPU.Build.0 = Debug|Any CPU {0EFC01B9-1F75-4BFD-ADB6-3FF18B2B9B5E}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs index dcc8352d..2a2a92a4 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs @@ -24,7 +24,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. - Last generated on: 2/9/2010 12:24:25 PM + Last generated on: 2/19/2010 4:32:24 PM */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { @@ -267,7 +267,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, xRows, xColumns, y, yRows, yColumns, 0.0, result); } /// @@ -834,7 +834,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result); } /// @@ -1401,7 +1401,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result); } /// @@ -1968,7 +1968,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs index 403d7666..02688cf0 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs @@ -28,7 +28,7 @@ /* This file is automatically generated - do not modify it. Change SafeNativeMethods.include instead. - Last generated on: 2/9/2010 12:22:29 PM + Last generated on: 2/19/2010 4:23:26 PM */ using System.Runtime.InteropServices; @@ -85,6 +85,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out]float[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out]double[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex32 alpha, Complex32[] x, Complex32[] y, ref Complex32 beta, [In, Out]Complex32[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex alpha, Complex[] x, Complex[] y, ref Complex beta, [In, Out]Complex[] c); + #endregion BLAS #region LAPACK diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 31a9f8eb..d3a5a9ac 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -330,12 +330,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/ Parallel.For(0, xRows, i => { + int ixIndex = i * xColumns; + int iyIndex = i * yColumns; for (int j = 0; j < yColumns; j++) { + double s = 0; for (int k = 0; k < xColumns; k++) { - result[j + yColumns * i] += xdata[k + xColumns * i] * ydata[j + yColumns * k]; + s += xdata[k + ixIndex] * ydata[j + yColumns * k]; } + result[j + iyIndex] = s; } }); } @@ -357,10 +361,9 @@ 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) { - // Choose nonsensical values for the number of rows and columns in c; fill them in depending + // Choose nonsensical values for the number of rows 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) @@ -386,7 +389,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } cRows = aColumns; - cColumns = bRows; } else if ((int)transposeA > 111) { @@ -401,7 +403,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } cRows = aColumns; - cColumns = bColumns; } else if ((int)transposeB > 111) { @@ -416,7 +417,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } cRows = aRows; - cColumns = bRows; } else { @@ -431,7 +431,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } cRows = aRows; - cColumns = bColumns; } if (alpha == 0.0 && beta == 0.0) diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 7e2771f3..83a7d283 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -24,7 +24,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. - Last generated on: 2/9/2010 12:24:29 PM + Last generated on: 2/19/2010 4:32:27 PM */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { @@ -266,7 +266,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, xRows, xColumns, y, yRows, yColumns, 0.0, result); } /// @@ -833,7 +833,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result); } /// @@ -1400,7 +1400,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result); } /// @@ -1967,7 +1967,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs index b09bf224..ee3e2435 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs @@ -28,7 +28,7 @@ /* This file is automatically generated - do not modify it. Change SafeNativeMethods.include instead. - Last generated on: 2/9/2010 12:22:33 PM + Last generated on: 2/19/2010 4:23:28 PM */ using System.Runtime.InteropServices; @@ -85,6 +85,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out]float[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out]double[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex32 alpha, Complex32[] x, Complex32[] y, ref Complex32 beta, [In, Out]Complex32[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex alpha, Complex[] x, Complex[] y, ref Complex beta, [In, Out]Complex[] c); + #endregion BLAS #region LAPACK diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include index c6bdd1ec..61966dac 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include +++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include @@ -281,7 +281,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, xRows, xColumns, y, yRows, yColumns, 0.0, result); } /// @@ -860,7 +860,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result); } /// @@ -1439,7 +1439,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result); } /// @@ -2018,7 +2018,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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) { - throw new NotImplementedException(); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include b/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include index f662b87c..115ea85b 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include +++ b/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include @@ -85,6 +85,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out]float[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out]double[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex32 alpha, Complex32[] x, Complex32[] y, ref Complex32 beta, [In, Out]Complex32[] c); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, ref Complex alpha, Complex[] x, Complex[] y, ref Complex beta, [In, Out]Complex[] c); + #endregion BLAS #region LAPACK diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index caeac622..ff584ad2 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -29,13 +29,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { using System; - using System.Collections.Generic; - using System.Globalization; - using Algorithms; - using Algorithms.LinearAlgebra; - using NumberTheory; using Properties; - using Threading; /// /// A Matrix class with dense storage. The underlying storage is a one dimensional array in column-major order. @@ -52,7 +46,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public DenseMatrix(int order) : base(order) { - Data = new double[order*order]; + Data = new double[order * order]; } /// @@ -116,13 +110,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double { for (int j = 0; j < array.GetLength(1); j++) { - At(i, j, array[i,j]); + At(i, j, array[i, j]); } } } /// - /// Gets or sets the matrix's data. + /// Gets the matrix's data. /// /// The matrix's data. internal double[] Data @@ -174,7 +168,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override double At(int row, int column) { - return Data[column * RowCount + row]; + return Data[(column * RowCount) + row]; } /// @@ -191,7 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public override void At(int row, int column, double value) { - Data[column * RowCount + row] = value; + Data[(column * RowCount) + row] = value; } /// @@ -299,8 +293,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the multiplication. /// If the other matrix is . /// If the result matrix is . - /// If this.Columns != other.Rows. - /// If the result matrix's dimensions are not the this.Rows x other.Columns. + /// If this.Columns != other.Rows. + /// If the result matrix's dimensions are not the this.Rows x other.Columns. public void Multiply(DenseMatrix other, DenseMatrix result) { if (other == null) @@ -323,17 +317,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentMatrixDimensions); } - if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) - { - Matrix tmp = result.CreateMatrix(result.RowCount, result.ColumnCount); - Multiply(other, tmp); - tmp.CopyTo(result); - } - else - { - Control.LinearAlgebraProvider.MatrixMultiply(this.Data, this.RowCount, this.ColumnCount, - other.Data, other.RowCount, other.ColumnCount, result.Data); - } + Control.LinearAlgebraProvider.MatrixMultiply(this.Data, this.RowCount, this.ColumnCount, other.Data, other.RowCount, other.ColumnCount, result.Data); } /// @@ -343,8 +327,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// choose the representation of either or depending on which /// is denser. /// The matrix to multiply with. - /// If this.Columns != other.Rows. + /// If this.Columns != other.Rows. /// If the other matrix is . + /// The result of multiplication. public Matrix Multiply(DenseMatrix other) { if (other == null) @@ -357,8 +342,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentMatrixDimensions); } - Matrix result = CreateMatrix(RowCount, other.ColumnCount); - Multiply(other, result); + DenseMatrix result = (DenseMatrix)CreateMatrix(RowCount, other.ColumnCount); + this.Multiply(other, result); return result; } @@ -369,7 +354,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The right matrix to multiply. /// The result of multiplication. /// If or is . - /// If the dimensions of or don't conform. + /// If the dimensions of or don't conform. public static DenseMatrix operator *(DenseMatrix leftSide, DenseMatrix rightSide) { if (leftSide == null) diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs index 04f87fe8..db90d14f 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // @@ -212,6 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { s += At(i, j) * rightSide[j]; } + result[i] = s; }); } @@ -280,6 +281,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { s += leftSide[i] * At(i, j); } + result[j] = s; }); } @@ -292,8 +294,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the multiplication. /// If the other matrix is . /// If the result matrix is . - /// If this.Columns != other.Rows. - /// If the result matrix's dimensions are not the this.Rows x other.Columns. + /// If this.Columns != other.Rows. + /// If the result matrix's dimensions are not the this.Rows x other.Columns. public virtual void Multiply(Matrix other, Matrix result) { if (other == null) @@ -324,18 +326,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double } else { - Parallel.For(0, this.RowCount, j => - { - for (int i = 0; i != other.ColumnCount; i++) + Parallel.For( + 0, + this.RowCount, + j => { - double s = 0; - for (int l = 0; l < this.ColumnCount; l++) + for (int i = 0; i != other.ColumnCount; i++) { - s += this.At(j, l) * other.At(l, i); + double s = 0; + for (int l = 0; l < this.ColumnCount; l++) + { + s += this.At(j, l) * other.At(l, i); + } + + result.At(j, i, s); } - result.At(j, i, s); - } - }); + }); } } @@ -346,8 +352,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// choose the representation of either or depending on which /// is denser. /// The matrix to multiply with. - /// If this.Columns != other.Rows. + /// If this.Columns != other.Rows. /// If the other matrix is . + /// The result of the multiplication. public virtual Matrix Multiply(Matrix other) { if (other == null) @@ -543,7 +550,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The right matrix to multiply. /// The result of multiplication. /// If or is . - /// If the dimensions of or don't conform. + /// If the dimensions of or don't conform. public static Matrix operator *(Matrix leftSide, Matrix rightSide) { if (leftSide == null)