diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 529d0949..cd094073 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -79,7 +79,10 @@ Matlab &lt &gt Wikipedia -kronecker +kronecker +blocksize +ipiv +dll diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 2159c0aa..ba2078df 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/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 ){ diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index 9406942d..9a9d1f8e 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/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 @@ /// This is equivalent to the GETRF and GETRI LAPACK routines. 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); } /// @@ -209,7 +220,28 @@ /// This is equivalent to the GETRI LAPACK routine. 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); } /// @@ -223,7 +255,27 @@ /// This is equivalent to the GETRF and GETRI LAPACK routines. 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); } /// @@ -238,7 +290,37 @@ /// This is equivalent to the GETRI LAPACK routine. 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); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.header.include b/src/Numerics/Algorithms/LinearAlgebra/native.header.include index ca4dcae8..3827a9aa 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.header.include +++ b/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 diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index a5e2f389..64005896 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/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 diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index d98152c3..39f48c5c 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -194,7 +194,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double Assert.AreEqual(_x[i] / _y[i], result[i]); } } - /// /// Can compute L1 norm. @@ -444,5 +443,110 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double Assert.AreEqual(ipiv[1], 2); Assert.AreEqual(ipiv[2], 2); } + + /// + /// Can compute the inverse of a matrix using LU factorization. + /// + [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); + } + + /// + /// Can compute the inverse of a matrix using LU factorization + /// using a previously factored matrix. + /// + [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); + } + + /// + /// Can compute the inverse of a matrix using LU factorization + /// with a work array. + /// + [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); + } + + /// + /// Can compute the inverse of a matrix using LU factorization + /// using a previously factored matrix with a work array. + /// + [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); + } + } }