matrix: fixed matrix multiplication bug where the call was doing an multiply and update, not just a multiply

build: corrected the default build task

fixed some stylecop errors

build/build.proj

@@ -1,4 +1,4 @@
-<Project xmlns="http://schemas.microsoft.com/developer/msbuild/2003" DefaultTargets="BuildAndTest" ToolsVersion="3.5">
+<Project xmlns="http://schemas.microsoft.com/developer/msbuild/2003" DefaultTargets="Test" ToolsVersion="3.5">
   <UsingTask TaskName="MSBuild.ExtensionPack.CodeQuality.StyleCop" AssemblyFile=".\MSBuild_Extension_Pack_3.5.5\StyleCop\MSBuild.ExtensionPack.StyleCop.dll" />
   <UsingTask TaskName="MSBuild.ExtensionPack.Framework.Assembly" AssemblyFile=".\MSBuild_Extension_Pack_3.5.5\MSBuild.ExtensionPack.dll" />
   <UsingTask TaskName="MSBuild.ExtensionPack.Communication.Twitter" AssemblyFile=".\MSBuild_Extension_Pack_3.5.5\MMSBuild.ExtensionPack.Communication.dll" />

src/FSharp/FSharp.fsproj

@@ -33,8 +33,8 @@
     <DocumentationFile>MathNet.Numerics.FSharp.XML</DocumentationFile>
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="FSharp.Core, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" />
     <Reference Include="mscorlib" />
-    <Reference Include="FSharp.Core" />
     <Reference Include="System" />
     <Reference Include="System.Core">
       <RequiredTargetFramework>3.5</RequiredTargetFramework>

src/FSharpExamples/FSharpExamples.fsproj

@@ -31,18 +31,12 @@
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="mscorlib" />
-    <Reference Include="FSharp.Core" />
     <Reference Include="System" />
     <Reference Include="System.Core">
       <RequiredTargetFramework>3.5</RequiredTargetFramework>
     </Reference>
   </ItemGroup>
   <ItemGroup>
-    <ProjectReference Include="..\FSharp\FSharp.fsproj">
-      <Name>FSharp</Name>
-      <Project>{37e8e802-a354-4114-bfc1-6e1357da605b}</Project>
-      <Private>True</Private>
-    </ProjectReference>
     <ProjectReference Include="..\Numerics\Numerics.csproj">
       <Name>Numerics</Name>
       <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>

src/FSharpUnitTests/FSharpUnitTests.fsproj

@@ -31,7 +31,6 @@
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="mscorlib" />
-    <Reference Include="FSharp.Core" />
     <Reference Include="System" />
     <Reference Include="System.Core">
       <RequiredTargetFramework>3.5</RequiredTargetFramework>
@@ -42,11 +41,6 @@
     <Compile Include="Program.fs" />
   </ItemGroup>
   <ItemGroup>
-    <ProjectReference Include="..\FSharp\FSharp.fsproj">
-      <Name>FSharp</Name>
-      <Project>{37e8e802-a354-4114-bfc1-6e1357da605b}</Project>
-      <Private>True</Private>
-    </ProjectReference>
     <ProjectReference Include="..\Numerics\Numerics.csproj">
       <Name>Numerics</Name>
       <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>

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);
+}

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.ActiveCfg = Release|Win32
-		{C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.Build.0 = Release|x64
+		{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

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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>

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

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)

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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>

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

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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>
@@ -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.</remarks>
         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);
         }
 
         /// <summary>

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

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;
 
     /// <summary>
     /// 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];
         }
 
         /// <summary>
@@ -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]);
                 }
             }
         }
 
         /// <summary>
-        /// Gets or sets the matrix's data.
+        /// Gets the matrix's data.
         /// </summary>
         /// <value>The matrix's data.</value>
         internal double[] Data
@@ -174,7 +168,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public override double At(int row, int column)
         {
-            return Data[column * RowCount + row];
+            return Data[(column * RowCount) + row];
         }
 
         /// <summary>
@@ -191,7 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </param>
         public override void At(int row, int column, double value)
         {
-            Data[column * RowCount + row] = value;
+            Data[(column * RowCount) + row] = value;
         }
 
         /// <summary>
@@ -299,8 +293,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="result">The result of the multiplication.</param>
         /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
         /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>
+        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
+        /// <exception cref="ArgumentException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
         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))
+            Control.LinearAlgebraProvider.MatrixMultiply(this.Data, this.RowCount, this.ColumnCount, other.Data, other.RowCount, other.ColumnCount, result.Data);
-            {
-                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);
-            }
         }
 
         /// <summary>
@@ -343,8 +327,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// choose the representation of either <paramref name="leftSide"/> or <paramref name="rightSide"/> depending on which
         /// is denser.</remarks>
         /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>        
+        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>        
         /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
+        /// <returns>The result of multiplication.</returns>
         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);
+            DenseMatrix result = (DenseMatrix)CreateMatrix(RowCount, other.ColumnCount);
-            Multiply(other, result);
+            this.Multiply(other, result);
             return result;
         }
 
@@ -369,7 +354,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="rightSide">The right matrix to multiply.</param>
         /// <returns>The result of multiplication.</returns>
         /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
-        /// <exception cref="NotConformableException">If the dimensions of <paramref name="leftSide"/> or <paramref name="rightSide"/> don't conform.</exception>
+        /// <exception cref="ArgumentException">If the dimensions of <paramref name="leftSide"/> or <paramref name="rightSide"/> don't conform.</exception>
         public static DenseMatrix operator *(DenseMatrix leftSide, DenseMatrix rightSide)
         {
             if (leftSide == null)

src/Numerics/LinearAlgebra/Double/Matrix.Arithmetic.cs

@@ -1,4 +1,4 @@
-// <copyright file="Matrix.cs" company="Math.NET">
+// <copyright file="Matrix.Arithmetic.cs" company="Math.NET">
 // 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
         /// <param name="result">The result of the multiplication.</param>
         /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
         /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
-        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>
+        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>
-        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
+        /// <exception cref="ArgumentException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
         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 =>
+                Parallel.For(
-                {
+                    0,
-                    for (int i = 0; i != other.ColumnCount; i++)
+                    this.RowCount,
+                    j =>
                     {
-                        double s = 0;
+                        for (int i = 0; i != other.ColumnCount; i++)
-                        for (int l = 0; l < this.ColumnCount; l++)
                         {
-                            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 <paramref name="leftSide"/> or <paramref name="rightSide"/> depending on which
         /// is denser.</remarks>
         /// <param name="other">The matrix to multiply with.</param>
-        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>        
+        /// <exception cref="ArgumentException">If <strong>this.Columns != other.Rows</strong>.</exception>        
         /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
+        /// <returns>The result of the multiplication.</returns>
         public virtual Matrix Multiply(Matrix other)
         {
             if (other == null)
@@ -543,7 +550,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="rightSide">The right matrix to multiply.</param>
         /// <returns>The result of multiplication.</returns>
         /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
-        /// <exception cref="NotConformableException">If the dimensions of <paramref name="leftSide"/> or <paramref name="rightSide"/> don't conform.</exception>
+        /// <exception cref="ArgumentException">If the dimensions of <paramref name="leftSide"/> or <paramref name="rightSide"/> don't conform.</exception>
         public static Matrix operator *(Matrix leftSide, Matrix rightSide)
         {
             if (leftSide == null)