native: added lu inverse

src/MathNet.Numerics.5.1.ReSharper

@@ -79,7 +79,10 @@ Matlab
 &lt
 &gt
 Wikipedia
-kronecker</UserWords>
+kronecker
+blocksize
+ipiv
+dll</UserWords>
             </CustomDictionary>
           </Dictionaries>
         </CustomDictionaries>

src/NativeWrappers/MKL/lapack.cpp

@@ -126,7 +126,43 @@ extern "C" {
         }
     }
 
-    DLLEXPORT void s_lu_inverse(int n, float a[], int ipiv[], float work[], int lwork)
+    DLLEXPORT void s_lu_inverse(int n, float a[], float work[], int lwork)
+    {
+        int* ipiv = new int[n];
+        int info;
+        SGETRF(&n,&n,a,&n,ipiv,&info);
+        SGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+		delete[] ipiv;
+    }
+
+    DLLEXPORT void d_lu_inverse(int n, double a[], double work[], int lwork)
+    {
+        int* ipiv = new int[n];
+        int info;
+        DGETRF(&n,&n,a,&n,ipiv,&info);
+        DGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+		delete[] ipiv;
+    }
+
+    DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], MKL_Complex8 work[], int lwork)
+    {
+        int* ipiv = new int[n];
+        int info;
+        CGETRF(&n,&n,a,&n,ipiv,&info);
+        CGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+		delete[] ipiv;
+    }
+
+    DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], MKL_Complex16 work[], int lwork)
+    {
+        int* ipiv = new int[n];
+        int info;
+        ZGETRF(&n,&n,a,&n,ipiv,&info);
+        ZGETRI(&n,a,&n,ipiv,work,&lwork,&info);
+		delete[] ipiv;
+    }
+
+	DLLEXPORT void s_lu_inverse_factored(int n, float a[], int ipiv[], float work[], int lwork)
     {
         int i;
         for(i = 0; i < n; ++i ){
@@ -140,7 +176,7 @@ extern "C" {
         }
     }
 
-    DLLEXPORT void d_lu_inverse(int n, double a[], int ipiv[], double work[], int lwork)
+    DLLEXPORT void d_lu_inverse_factored(int n, double a[], int ipiv[], double work[], int lwork)
     {
         int i;
         for(i = 0; i < n; ++i ){
@@ -155,7 +191,7 @@ extern "C" {
         }
     }
 
-    DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork)
+    DLLEXPORT void c_lu_inverse_factored(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork)
     {
         int i;
         for(i = 0; i < n; ++i ){
@@ -170,7 +206,7 @@ extern "C" {
         }
     }
 
-    DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork)
+    DLLEXPORT void z_lu_inverse_factored(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork)
     {
         int i;
         for(i = 0; i < n; ++i ){

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

@@ -23,10 +23,10 @@
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength);
             }
 
-			if (!ReferenceEquals(y, result))
-			{
-				Array.Copy(y, 0, result, 0, y.Length);
-			}
+            if (!ReferenceEquals(y, result))
+            {
+                Array.Copy(y, 0, result, 0, y.Length);
+            }
 
             if (alpha == <#=zero#>)
             {
@@ -50,10 +50,10 @@
                 throw new ArgumentNullException("x");
             } 
             
-			if (!ReferenceEquals(x, result))
-			{
-				Array.Copy(x, 0, result, 0, x.Length);
-			}
+            if (!ReferenceEquals(x, result))
+            {
+                Array.Copy(x, 0, result, 0, x.Length);
+            }
 
             if (alpha == <#=one#>)
             {
@@ -197,7 +197,18 @@
         /// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
         public override void LUInverse(<#=dataType#>[] a, int order)
         {
-            throw new NotImplementedException();
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
+            if (a.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
+            }
+
+            var work = new <#=dataType#>[order];
+            SafeNativeMethods.<#=prefix#>_lu_inverse(order, a, work, order);        
         }
 
         /// <summary>
@@ -209,7 +220,28 @@
         /// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
         public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv)
         {
-            throw new NotImplementedException();
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
+            if (ipiv == null)
+            {
+                throw new ArgumentNullException("ipiv");
+            }
+
+            if (a.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
+            }
+
+            if (ipiv.Length != order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
+            }
+
+            var work = new <#=dataType#>[order];
+            SafeNativeMethods.<#=prefix#>_lu_inverse_factored(order, a, ipiv, work, order);        
         }
 
         /// <summary>
@@ -223,7 +255,27 @@
         /// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
         public override void LUInverse(<#=dataType#>[] a, int order, <#=dataType#>[] work)
         {
-            throw new NotImplementedException();
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
+            if (a.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
+            }
+
+            if (work == null)
+            {
+                throw new ArgumentNullException("work");
+            }
+
+            if (work.Length < order)
+            {
+                throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
+            }
+
+            SafeNativeMethods.<#=prefix#>_lu_inverse(order, a, work, work.Length);        
         }
 
         /// <summary>
@@ -238,7 +290,37 @@
         /// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
         public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] work)
         {
-            throw new NotImplementedException();
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
+
+            if (ipiv == null)
+            {
+                throw new ArgumentNullException("ipiv");
+            }
+
+            if (a.Length != order * order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
+            }
+
+            if (ipiv.Length != order)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
+            }
+
+            if (work == null)
+            {
+                throw new ArgumentNullException("work");
+            }
+
+            if (work.Length < order)
+            {
+                throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
+            }
+
+            SafeNativeMethods.<#=prefix#>_lu_inverse_factored(order, a, ipiv, work, order);        
         }
 
         /// <summary>

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

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2010 Math.NET
+// Copyright (c) 2009-2011 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation

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

@@ -2,7 +2,7 @@
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
 //
-// Copyright (c) 2009 Math.NET
+// Copyright (c) 2009-2010 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -137,4 +137,28 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
         internal static extern void z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv);
 
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork);
+    
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork);
+        
         #endregion LAPACK

src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

@@ -194,7 +194,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
                 Assert.AreEqual(_x[i] / _y[i], result[i]);
             }
         }
-
         
         /// <summary>
         /// Can compute L1 norm.
@@ -444,5 +443,110 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
             Assert.AreEqual(ipiv[1], 2);
             Assert.AreEqual(ipiv[2], 2);
         }
+
+        /// <summary>
+        /// Can compute the inverse of a matrix using LU factorization.
+        /// </summary>
+        [Test]
+        public void CanComputeLuInverse()
+        {
+            var matrix = _matrices["Square3x3"];
+            var a = new double[matrix.RowCount * matrix.RowCount];
+            Array.Copy(matrix.Data, a, a.Length);
+            
+            Provider.LUInverse(a, matrix.RowCount);
+
+            AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[1],  -0.909090909090908, 14);
+            AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 14);
+            AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 14);
+            AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 14);
+            AssertHelpers.AlmostEqual(a[6],  -0.113636363636364, 14);
+            AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 14);
+            AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14);
+        }
+
+        /// <summary>
+        /// Can compute the inverse of a matrix using LU factorization
+        /// using a previously factored matrix.
+        /// </summary>
+        [Test]
+        public void CanComputeLuInverseOnFactoredMatrix()
+        {
+            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);
+            Provider.LUInverseFactored(a, matrix.RowCount, ipiv);
+
+            AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[1], -0.909090909090908, 14);
+            AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 14);
+            AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 14);
+            AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 14);
+            AssertHelpers.AlmostEqual(a[6], -0.113636363636364, 14);
+            AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 14);
+            AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14);
+        }
+
+        /// <summary>
+        /// Can compute the inverse of a matrix using LU factorization
+        /// with a work array.
+        /// </summary>
+        [Test]
+        public void CanComputeLuInverseWithWorkArray()
+        {
+            var matrix = _matrices["Square3x3"];
+            var a = new double[matrix.RowCount * matrix.RowCount];
+            Array.Copy(matrix.Data, a, a.Length);
+
+            var work = new double[matrix.RowCount];
+            Provider.LUInverse(a, matrix.RowCount, work);
+
+            AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[1], -0.909090909090908, 14);
+            AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 14);
+            AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 14);
+            AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 14);
+            AssertHelpers.AlmostEqual(a[6], -0.113636363636364, 14);
+            AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 14);
+            AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14);
+        }
+
+        /// <summary>
+        /// Can compute the inverse of a matrix using LU factorization
+        /// using a previously factored matrix with a work array.
+        /// </summary>
+        [Test]
+        public void CanComputeLuInverseOnFactoredMatrixWithWorkArray()
+        {
+            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);
+
+            var work = new double[matrix.RowCount];
+            Provider.LUInverseFactored(a, matrix.RowCount, ipiv, work);
+
+            AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[1], -0.909090909090908, 14);
+            AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 14);
+            AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 14);
+            AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 14);
+            AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 14);
+            AssertHelpers.AlmostEqual(a[6], -0.113636363636364, 14);
+            AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 14);
+            AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14);
+        }
+
     }
 }