native: added lu

bug:BulirschStoerRationalInterpolation was always returing NaN

build/mathnet-numerics.shfbproj

@@ -21,19 +21,26 @@
     <PresentationStyle>vs2005</PresentationStyle>
     <DocumentationSources>
       <DocumentationSource sourceFile="..\out\lib\Net40\MathNet.Numerics.dll" xmlns="" />
-<DocumentationSource sourceFile="..\out\lib\Net40\MathNet.Numerics.xml" xmlns="" /></DocumentationSources>
+      <DocumentationSource sourceFile="..\out\lib\Net40\MathNet.Numerics.xml" xmlns="" />
+    </DocumentationSources>
     <Preliminary>True</Preliminary>
     <ProjectSummary>Numerical Library for .NET</ProjectSummary>
     <CopyrightText>Copyright %28c%29 2009-2010 Math.NET</CopyrightText>
-    <RootNamespaceContainer>True</RootNamespaceContainer>
+    <RootNamespaceContainer>False</RootNamespaceContainer>
-    <HelpTitle>Math.NET Numerics</HelpTitle>
+    <HelpTitle>Math.NET Numerics Beta 1</HelpTitle>
     <RootNamespaceTitle>API Reference</RootNamespaceTitle>
     <IncludeFavorites>False</IncludeFavorites>
     <PlugInConfigurations>
-      <PlugInConfig id="Hierarchical Table of Contents" enabled="True" xmlns="">
+      <PlugInConfig id="Hierarchical Table of Contents" enabled="True">
         <configuration />
       </PlugInConfig>
-<PlugInConfig id="Additional Content Only" enabled="True" xmlns=""><configuration /></PlugInConfig></PlugInConfigurations>
+<PlugInConfig id="Additional Content Only" enabled="False"><configuration /></PlugInConfig></PlugInConfigurations>
+    <ContentPlacement>AboveNamespaces</ContentPlacement>
+    <ApiFilter>
+      <Filter entryType="Namespace" fullName="MathNet.Numerics.Algorithms.LinearAlgebra.Mkl" isExposed="False" xmlns="" />
+    </ApiFilter>
+    <KeepLogFile>False</KeepLogFile>
+    <HelpFileFormat>HtmlHelp1, MSHelpViewer</HelpFileFormat>
   </PropertyGroup>
   <!-- There are no properties for these groups.  AnyCPU needs to appear in
        order for Visual Studio to perform the build.  The others are optional
@@ -55,10 +62,10 @@
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|Win32' ">
   </PropertyGroup>
   <ItemGroup>
-    <Content Include="License.txt" />
+    <SiteMap Include="SiteMap.sitemap" />
   </ItemGroup>
   <ItemGroup>
-    <SiteMap Include="SiteMap.sitemap" />
+    <Content Include="License.txt" />
   </ItemGroup>
   <!-- Import the SHFB build targets -->
   <Import Project="$(SHFBROOT)\SandcastleHelpFileBuilder.targets" />

src/MSUnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

@@ -73,11 +73,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
         /// </summary>
         private readonly IDictionary<string, DenseMatrix> _matrices = new Dictionary<string, DenseMatrix>
                                                           {
-                                                              { "Singular3x3", new DenseMatrix(new[,] { { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 } }) },
+                                                              { "Singular3x3", new DenseMatrix(new[,] { { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 }, { 1.0, 1.0, 2.0 } }) }, 
-                                                              { "Square3x3",  new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 }, { -4.4, 5.5, 6.6 } }) },
+                                                              { "Square3x3",  new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 }, { -4.4, 5.5, 6.6 } }) }, 
-                                                              { "Square4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { 0.0, 1.1, 2.2, 3.3 }, { 1.0, 2.1, 6.2, 4.3 }, { -4.4, 5.5, 6.6, -7.7 } }) },
+                                                              { "Square4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { 0.0, 1.1, 2.2, 3.3 }, { 1.0, 2.1, 6.2, 4.3 }, { -4.4, 5.5, 6.6, -7.7 } }) }, 
-                                                              { "Singular4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 } }) },
+                                                              { "Singular4x4", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 }, { -1.1, -2.2, -3.3, -4.4 } }) }, 
-                                                              { "Tall3x2", new DenseMatrix(new[,] { { -1.1, -2.2 }, { 0.0, 1.1 }, { -4.4, 5.5 } }) },
+                                                              { "Tall3x2", new DenseMatrix(new[,] { { -1.1, -2.2 }, { 0.0, 1.1 }, { -4.4, 5.5 } }) }, 
                                                               { "Wide2x3", new DenseMatrix(new[,] { { -1.1, -2.2, -3.3 }, { 0.0, 1.1, 2.2 } }) }
                                                           };
 
@@ -406,5 +406,34 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
             Assert.AreEqual(matrix[14], 0);
             Assert.AreEqual(matrix[15], 1);
         }
+
+        /// <summary>
+        /// Can compute the LU factor of a matrix.
+        /// </summary>
+        [TestMethod]
+        public void CanComputeLuFactor()
+        {
+            var matrix = _matrices["Square3x3"];
+            var a = new double[matrix.RowCount * matrix.RowCount];
+            Array.Copy(matrix.Data, a, a.Length);
+
+            var ipiv = new int[matrix.RowCount];
+
+            Provider.LUFactor(a, matrix.RowCount, ipiv);
+
+            AssertHelpers.AlmostEqual(a[0], -4.4, 15);
+            AssertHelpers.AlmostEqual(a[1], 0.25, 15);
+            AssertHelpers.AlmostEqual(a[2], 0, 15);
+            AssertHelpers.AlmostEqual(a[3], 5.5, 15);
+            AssertHelpers.AlmostEqual(a[4], -3.575, 15);
+            AssertHelpers.AlmostEqual(a[5], -0.307692307692308, 15);
+            AssertHelpers.AlmostEqual(a[6], 6.6, 15);
+            AssertHelpers.AlmostEqual(a[7], -4.95, 15);
+            AssertHelpers.AlmostEqual(a[8], 0.676923076923077, 15);
+
+            Assert.AreEqual(ipiv[0], 2);
+            Assert.AreEqual(ipiv[1], 2);
+            Assert.AreEqual(ipiv[2], 2);
+        }
     }
 }

src/NativeWrappers/MKL/lapack.cpp

@@ -4,89 +4,476 @@
 
 extern "C" {
 
-	DLLEXPORT int s_cholesky_factor(int n, float a[]){
+    DLLEXPORT int s_cholesky_factor(int n, float a[]){
-		char uplo = 'L';
+        char uplo = 'L';
-		int info = 0;
+        int info = 0;
-		SPOTRF(&uplo, &n, a, &n, &info);
+        SPOTRF(&uplo, &n, a, &n, &info);
-		for (int i = 0; i < n; ++i)
+        for (int i = 0; i < n; ++i)
-		{
+        {
-			int index = i * n;
+            int index = i * n;
-			for (int j = 0; j < n && i > j; ++j)
+            for (int j = 0; j < n && i > j; ++j)
-			{
+            {
-				a[index + j] = 0;
+                a[index + j] = 0;
-			}
+            }
-		}
+        }
-		return info;
+        return info;
-	}
+    }
-
+
-	DLLEXPORT int d_cholesky_factor(int n, double* a){
+    DLLEXPORT int d_cholesky_factor(int n, double* a){
-		char uplo = 'L';
+        char uplo = 'L';
-		int info = 0;
+        int info = 0;
-		DPOTRF(&uplo, &n, a, &n, &info);
+        DPOTRF(&uplo, &n, a, &n, &info);
-		for (int i = 0; i < n; ++i)
+        for (int i = 0; i < n; ++i)
-		{
+        {
-			int index = i * n;
+            int index = i * n;
-			for (int j = 0; j < n && i > j; ++j)
+            for (int j = 0; j < n && i > j; ++j)
-			{
+            {
-				a[index + j] = 0;
+                a[index + j] = 0;
-			}
+            }
-		}
+        }
-		return info;
+        return info;
-	}
+    }
-
+
-	DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){
+    DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){
-		char uplo = 'L';
+        char uplo = 'L';
-		int info = 0;
+        int info = 0;
-		Complex8 zero;
+        Complex8 zero;
-		zero.real = 0.0;
+        zero.real = 0.0;
-		zero.real = 0.0;
+        zero.real = 0.0;
-		CPOTRF(&uplo, &n, a, &n, &info);
+        CPOTRF(&uplo, &n, a, &n, &info);
-		for (int i = 0; i < n; ++i)
+        for (int i = 0; i < n; ++i)
-		{
+        {
-			int index = i * n;
+            int index = i * n;
-			for (int j = 0; j < n && i > j; ++j)
+            for (int j = 0; j < n && i > j; ++j)
-			{
+            {
-				a[index + j] = zero;
+                a[index + j] = zero;
-			}
+            }
-		}
+        }
-		return info;
+        return info;
-	}
+    }
-
+
-	DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){
+    DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){
-		char uplo = 'L';
+        char uplo = 'L';
-		int info = 0;
+        int info = 0;
-		Complex16 zero;
+        Complex16 zero;
-		zero.real = 0.0;
+        zero.real = 0.0;
-		zero.real = 0.0;
+        zero.real = 0.0;
-		ZPOTRF(&uplo, &n, a, &n, &info);
+        ZPOTRF(&uplo, &n, a, &n, &info);
-		for (int i = 0; i < n; ++i)
+        for (int i = 0; i < n; ++i)
-		{
+        {
-			int index = i * n;
+            int index = i * n;
-			for (int j = 0; j < n && i > j; ++j)
+            for (int j = 0; j < n && i > j; ++j)
-			{
+            {
-				a[index + j] = zero;
+                a[index + j] = zero;
-			}
+            }
-		}
+        }
-		return info;
+        return info;
-	}
+    }
-
+
-	DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[])
+    DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[])
-	{
+    {
-		return SLANGE(&norm, &m, &n, a, &m, work);
+        return SLANGE(&norm, &m, &n, a, &m, work);
-	}
+    }
-
+
-	DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[])
+    DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[])
-	{
+    {
-		return DLANGE(&norm, &m, &n, a, &m, work);
+        return DLANGE(&norm, &m, &n, a, &m, work);
-	}
+    }
-	
+
-	DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[])
+    DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[])
-	{
+    {
-		return CLANGE(&norm, &m, &n, a, &m, work);
+        return CLANGE(&norm, &m, &n, a, &m, work);
-	}
+    }
-	
+
-	DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[])
+    DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[])
-	{
+    {
-		return ZLANGE(&norm, &m, &n, a, &m, work);
+        return ZLANGE(&norm, &m, &n, a, &m, work);
-	}
+    }
+
+    DLLEXPORT void s_lu_factor(int m, float a[], int ipiv[])
+    {
+        int info;
+        SGETRF(&m,&m,a,&m,ipiv,&info);
+        for(int i = 0; i < m; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void d_lu_factor(int m, double a[], int ipiv[])
+    {
+        int info;
+        DGETRF(&m,&m,a,&m,ipiv,&info);
+        for(int i = 0; i < m; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void c_lu_factor(int m, MKL_Complex8 a[], int ipiv[])
+    {
+        int info;
+        CGETRF(&m,&m,a,&m,ipiv,&info);
+        for(int i = 0; i < m; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void z_lu_factor(int m, MKL_Complex16 a[], int ipiv[])
+    {
+        int info;
+        ZGETRF(&m,&m,a,&m,ipiv,&info);
+        for(int i = 0; i < m; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void s_lu_inverse(int n, float a[], int ipiv[], float work[], int lwork)
+    {
+        int i;
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+        int info;
+        SGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void d_lu_inverse(int n, double a[], int ipiv[], double work[], int lwork)
+    {
+        int i;
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        int info;
+        DGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork)
+    {
+        int i;
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        int info;
+        CGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork)
+    {
+        int i;
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        int info;
+        ZGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void s_lu_solve(int n, int nrhs, float a[], int ipiv[], float b[])
+    {
+        int info;
+        int i;    
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        char trans = 'N';
+        SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void  d_lu_solve(int n, int nrhs, double a[], int ipiv[], double b[])
+    {
+        int info;
+        int i;    
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        char trans = 'N';
+        DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[])
+    {
+        int info;
+        int i;    
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        char trans = 'N';
+        CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[])
+    {
+        int info;
+        int i;    
+        for(i = 0; i < n; ++i ){
+            ipiv[i] += 1;
+        }
+
+        char trans = 'N';
+        ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info);
+        for(i = 0; i < n; ++i ){
+            ipiv[i] -= 1;
+        }
+    }
+
+    DLLEXPORT void s_cholesky_solve(int n, int nrhs, float a[], float b[])
+    {
+        char uplo = 'L';
+        int info = 0;
+        SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
+    }
+
+    DLLEXPORT void d_cholesky_solve(int n, int nrhs, double a[], double b[])
+    {
+        char uplo = 'L';
+        int info = 0;
+        DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
+    }
+
+    DLLEXPORT void c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[])
+    {
+        char uplo = 'L';
+        int info = 0;
+        CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
+    }
+
+    DLLEXPORT void z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[])
+    {
+        char uplo = 'L';
+        int info = 0;
+        ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info);
+    }
+
+    DLLEXPORT void s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len)
+    {
+        int info = 0;
+        SGEQRF(&m, &n, r, &m, tau, work, &len, &info);
+
+        for (int i = 0; i < m; ++i)
+        {
+            for (int j = 0; j < m && j < n; ++j)
+            {
+                if (i > j)
+                {
+                    q[j * m + i] = r[j * m + i];
+                }
+            }
+        }
+
+        //compute the q elements explicitly
+        if (m <= n)
+        {
+            SORGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
+        }
+        else
+        {
+            SORGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
+        }
+    }
+
+    DLLEXPORT void d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len)
+    {
+        int info = 0;
+        DGEQRF(&m, &n, r, &m, tau, work, &len, &info);
+
+        for (int i = 0; i < m; ++i)
+        {
+            for (int j = 0; j < m && j < n; ++j)
+            {
+                if (i > j)
+                {
+                    q[j * m + i] = r[j * m + i];
+                }
+            }
+        }
+
+        //compute the q elements explicitly
+        if (m <= n)
+        {
+            DORGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
+        }
+        else
+        {
+            DORGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
+        }
+    }
+
+    DLLEXPORT void c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], int len)
+    {
+        int info = 0;
+        CGEQRF(&m, &n, r, &m, tau, work, &len, &info);
+
+        for (int i = 0; i < m; ++i)
+        {
+            for (int j = 0; j < m && j < n; ++j)
+            {
+                if (i > j)
+                {
+                    q[j * m + i] = r[j * m + i];
+                }
+            }
+        }
+
+        //compute the q elements explicitly
+        if (m <= n)
+        {
+            CUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
+        }
+        else
+        {
+            CUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
+        }
+    }
+
+    DLLEXPORT void z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], int len)
+    {
+        int info = 0;
+        ZGEQRF(&m, &n, r, &m, tau, work, &len, &info);
+
+        for (int i = 0; i < m; ++i)
+        {
+            for (int j = 0; j < m && j < n; ++j)
+            {
+                if (i > j)
+                {
+                    q[j * m + i] = r[j * m + i];
+                }
+            }
+        }
+
+        //compute the q elements explicitly
+        if (m <= n)
+        {
+            ZUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info);
+        }
+        else
+        {
+            ZUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info);
+        }
+    }
+
+    DLLEXPORT void s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len)
+    {
+        char side ='L';
+        char tran = 'T';
+        int info = 0;
+        SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
+        cblas_strsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m);
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < bn; ++j)
+            {
+                x[j * n + i] = b[j * m + i];
+            }
+        }
+    }
+
+    DLLEXPORT void d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len)
+    {
+        char side ='L';
+        char tran = 'T';
+        int info = 0;
+        DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
+        cblas_dtrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m);
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < bn; ++j)
+            {
+                x[j * n + i] = b[j * m + i];
+            }
+        }
+    }
+
+    DLLEXPORT void c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], int len)
+    {
+        char side ='L';
+        char tran = 'T';
+        int info = 0;
+        CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
+        MKL_Complex8 one;
+        one.real = 1.0;
+        cblas_ctrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m);
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < bn; ++j)
+            {
+                x[j * n + i] = b[j * m + i];
+            }
+        }
+    }
+
+    DLLEXPORT void z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], int len)
+    {
+        char side ='L';
+        char tran = 'T';
+        int info = 0;
+        ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info);
+        MKL_Complex16 one;
+        one.real = 1.0;
+        cblas_ztrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m);
+        for (int i = 0; i < n; ++i)
+        {
+            for (int j = 0; j < bn; ++j)
+            {
+                x[j * n + i] = b[j * m + i];
+            }
+        }
+    }
+
+    DLLEXPORT void s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len)
+    {
+        int info = 0;
+        char job = compute_vectors ? 'A' : 'N';
+        SGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
+    }
+
+    DLLEXPORT void d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len)
+    {
+        int info = 0;
+        char job = compute_vectors ? 'A' : 'N';
+        DGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info);
+    }
+
+    DLLEXPORT void c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], int len, float rwork[])
+    {
+        int info = 0;
+        char job = compute_vectors ? 'A' : 'N';
+        CGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info);
+    }
+
+    DLLEXPORT void z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[], MKL_Complex16 work[], int len, double rwork[])
+    {
+        int info = 0;
+        char job = compute_vectors ? 'A' : 'N';
+        ZGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info);
+    }
 }

src/NativeWrappers/Windows/MKLWrapper32Tests/MKLWrapper32Tests.csproj

@@ -55,8 +55,9 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
   </PropertyGroup>
   <ItemGroup>
-    <Reference Include="MathNet.Numerics">
+    <Reference Include="MathNet.Numerics, Version=2010.11.21.836, Culture=neutral, PublicKeyToken=cd8b63ad3d691a37, processorArchitecture=MSIL">
-      <HintPath>..\..\..\..\out\debug\Net40\MathNet.Numerics.dll</HintPath>
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\..\out\lib\Net40\MathNet.Numerics.dll</HintPath>
     </Reference>
     <Reference Include="Microsoft.VisualStudio.QualityTools.UnitTestFramework, Version=10.0.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a, processorArchitecture=MSIL" />
     <Reference Include="System" />

src/Numerics/Algorithms/LinearAlgebra/native.generic.include

@@ -154,7 +154,27 @@
         /// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
         public override void LUFactor(<#=dataType#>[] data, int order, int[] ipiv)
         {
-            throw new NotImplementedException();
+            if (data == null)
+            {
+                throw new ArgumentNullException("data");
+            }
+
+            if (ipiv == null)
+            {
+                throw new ArgumentNullException("ipiv");
+            }
+
+            if (data.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "data");
+            }
+
+            if (ipiv.Length != order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
+            }
+            
+            SafeNativeMethods.<#=prefix#>_lu_factor(order, data, ipiv);
         }
 
         /// <summary>
@@ -284,6 +304,11 @@
                 throw new ArgumentException(Resources.ArgumentMustBePositive, "order");
             }
 
+            if (a.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
+            }
+
             SafeNativeMethods.<#=prefix#>_cholesky_factor(order, a);
         }
 

src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include

@@ -101,6 +101,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
         
         #region LAPACK
 
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work);
+
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
         internal static extern void s_cholesky_factor(int n, [In, Out] float[] a);
 
@@ -112,17 +124,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
         internal static extern void z_cholesky_factor(int n, [In, Out] Complex[] a);
-
+        
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-        internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work);
+        internal static extern void s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-        internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work);
+        internal static extern void d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-        internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work);
+        internal static extern void c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-        internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work);
+        internal static extern void z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv);
 
         #endregion LAPACK

src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs

@@ -168,7 +168,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms
                     double ho = (_points[i] - t) * d[i] / hp;
 
                     double den = ho - c[i + 1];
-                    if (den.AlmostEqual(0.0))
+                    if (den == 0.0)
                     {
                         return double.NaN; // zero-div, singularity
                     }

src/Numerics/Permutation.cs

@@ -110,7 +110,7 @@ namespace MathNet.Numerics
         /// From wikipedia: the permutation 12043 has the inversions (0,2), (1,2) and (3,4). This would be
         /// encoded using the array [22244].
         /// </example>
-        /// <param name="inv">The set of inversions to contruct the permutation from.</param>
+        /// <param name="inv">The set of inversions to construct the permutation from.</param>
         /// <returns>A permutation generated from a sequence of inversions.</returns>
         public static Permutation FromInversions(int[] inv)
         {
@@ -176,7 +176,7 @@ namespace MathNet.Numerics
         /// </summary>
         /// <param name="indices">An array which represents where each integer is permuted too: indices[i] represents that integer i
         /// is permuted to location indices[i].</param>
-        /// <returns>True if <paramref name="indices"/> represents a proper permutation, false otherwise.</returns>
+        /// <returns>True if <paramref name="indices"/> represents a proper permutation, <c>false</c> otherwise.</returns>
         static private bool CheckForProperPermutation(int[] indices)
         {
             var idxCheck = new bool[indices.Length];