diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 33ab37e3..3f9e62c0 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -523,6 +523,7 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; + r[j * m + i] = 0.0f; } } } @@ -552,6 +553,7 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; + r[j * m + i] = 0.0; } } } @@ -573,6 +575,10 @@ extern "C" { { int info = 0; CGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + MKL_Complex8 zero; + zero.real = 0.0f; + zero.imag = 0.0f; for (int i = 0; i < m; ++i) { @@ -581,6 +587,7 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; + r[j * m + i] = zero; } } } @@ -603,6 +610,10 @@ extern "C" { int info = 0; ZGEQRF(&m, &n, r, &m, tau, work, &len, &info); + MKL_Complex16 zero; + zero.real = 0.0; + zero.imag = 0.0; + for (int i = 0; i < m; ++i) { for (int j = 0; j < m && j < n; ++j) @@ -610,6 +621,7 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; + r[j * m + i] = zero; } } } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index e0a9b4c6..7ff780e7 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -316,13 +316,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// Computes the QR factorization of A. /// /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. + /// it is overwritten with the R matrix of the QR factorization. /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - void QRFactor(T[] r, int rowsR, int columnsR, T[] q); + void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] tau); /// /// Computes the QR factorization of A. @@ -333,11 +335,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. - void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] work); + void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] tau, T[] work); /// /// Solves A*X=B for X using QR factorization of A. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 06e497a9..5ca428dd 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -1359,8 +1359,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q) + public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) { if (r == null) { @@ -1374,16 +1376,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } var work = new Complex[rowsR * rowsR]; - QRFactor(r, rowsR, columnsR, q, work); + QRFactor(r, rowsR, columnsR, q, tau, work); } /// @@ -1395,11 +1402,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] work) + public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau, Complex[] work) { if (r == null) { @@ -1418,12 +1427,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (work.Length < rowsR * rowsR) @@ -1681,8 +1695,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// The B matrix. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 8a8f6b99..89e4b98f 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -1359,8 +1359,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q) + public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) { if (r == null) { @@ -1374,16 +1376,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } var work = new Complex32[rowsR * rowsR]; - QRFactor(r, rowsR, columnsR, q, work); + QRFactor(r, rowsR, columnsR, q, tau, work); } /// @@ -1395,11 +1402,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] work) + public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau, Complex32[] work) { if (r == null) { @@ -1418,12 +1427,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (work.Length < rowsR * rowsR) @@ -1681,8 +1695,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// The B matrix. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 10ec9bd3..25cc4b04 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -1353,8 +1353,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q) + public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) { if (r == null) { @@ -1368,16 +1370,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } var work = new double[rowsR * rowsR]; - QRFactor(r, rowsR, columnsR, q, work); + QRFactor(r, rowsR, columnsR, q, tau, work); } /// @@ -1389,11 +1396,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] work) + public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau, double[] work) { if (r == null) { @@ -1412,12 +1421,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (work.Length < rowsR * rowsR) @@ -1676,8 +1690,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// The B matrix. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index d9f99f7d..317beeac 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -1354,8 +1354,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q) + public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) { if (r == null) { @@ -1369,16 +1371,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } var work = new float[rowsR * rowsR]; - QRFactor(r, rowsR, columnsR, q, work); + QRFactor(r, rowsR, columnsR, q, tau, work); } /// @@ -1390,11 +1397,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] work) + public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau, float[] work) { if (r == null) { @@ -1413,12 +1422,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (work.Length < rowsR * rowsR) @@ -1677,8 +1691,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// The B matrix. diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index 7da5fb3f..e3bc980e 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include @@ -517,11 +517,39 @@ /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q) + public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] tau) { - throw new NotImplementedException(); + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rowsR * columnsR) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (q.Length != rowsR * rowsR) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); + } + + var work = new <#=dataType#>[columnsR * Control.BlockSize]; + SafeNativeMethods.<#=prefix#>_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// @@ -533,14 +561,52 @@ /// The number of columns in the A matrix. /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] work) + public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] tau, <#=dataType#>[] work) { - throw new NotImplementedException(); + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (work == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rowsR * columnsR) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); + } + + if (tau.Length < Math.Min(rowsR, columnsR)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (q.Length != rowsR * rowsR) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); + } + + if (work.Length < columnsR * Control.BlockSize) + { + work[0] = columnsR * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + SafeNativeMethods.<#=prefix#>_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index adae155b..59af2294 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include @@ -187,27 +187,39 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, float[] b); + internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, double[] b); + internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, Complex32[] b); + internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, Complex[] b); + internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); #endregion LAPACK diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index 13dceba7..abd698ac 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -46,6 +46,11 @@ namespace MathNet.Numerics private static int _numberOfThreads = Environment.ProcessorCount; #endif + /// + /// Initial block size for the native linear algebra provider. + /// + private static int _blockSize = 512; + /// /// Initializes static members of the Control class. /// @@ -60,7 +65,7 @@ namespace MathNet.Numerics /// /// Gets or sets a value indicating whether the distribution classes check validate each parameter. /// For the multivariate distributions this could involve an expensive matrix factorization. - /// The default setting of this property is true. + /// The default setting of this property is true. /// public static bool CheckDistributionParameters { get; set; } @@ -109,5 +114,26 @@ namespace MathNet.Numerics } #endif } + + /// + /// Gets or sets the the block size to use for the native linear + /// algebra provider. + /// + /// The block size. Must be at least 32. + public static int BlockSize + { + get + { + return _blockSize; + } + + set + { + if (_blockSize > 31) + { + _blockSize = value; + } + } + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index 22dcd05f..a5dd9321 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs @@ -46,6 +46,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// public class DenseQR : QR { + /// + /// Gets or sets Tau vector. Contains additional information on Q - used for native solver. + /// + public Complex[] Tau + { + get; + set; + } + /// /// Initializes a new instance of the class. This object will compute the /// QR factorization when the constructor is called and cache it's factorization. @@ -67,7 +76,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization MatrixR = matrix.Clone(); MatrixQ = new DenseMatrix(matrix.RowCount); - Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data); + Tau = new Complex[Math.Min(matrix.RowCount, matrix.ColumnCount)]; + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau); } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 64ef6526..403c3840 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs @@ -46,6 +46,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// public class DenseQR : QR { + /// + /// Gets or sets Tau vector. Contains additional information on Q - used for native solver. + /// + public Complex32[] Tau + { + get; + set; + } + /// /// Initializes a new instance of the class. This object will compute the /// QR factorization when the constructor is called and cache it's factorization. @@ -67,7 +76,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization MatrixR = matrix.Clone(); MatrixQ = new DenseMatrix(matrix.RowCount); - Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data); + Tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)]; + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 48017348..1e626da7 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -45,6 +45,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// public class DenseQR : QR { + /// + /// Gets or sets Tau vector. Contains additional information on Q - used for native solver. + /// + public double[] Tau + { + get; + set; + } + /// /// Initializes a new instance of the class. This object will compute the /// QR factorization when the constructor is called and cache it's factorization. @@ -54,11 +63,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// If row count is less then column count public DenseQR(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount < matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); @@ -66,7 +70,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization MatrixR = matrix.Clone(); MatrixQ = new DenseMatrix(matrix.RowCount); - Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data); + Tau = new double[Math.Min(matrix.RowCount, matrix.ColumnCount)]; + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs index 44d3c8d1..f9fcdc6a 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs @@ -91,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Generate column from initial matrix to work array /// /// Initial matrix - /// The firts row + /// The first row /// Column index /// Generated vector private static double[] GenerateColumn(Matrix a, int row, int column) diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index a0fabfc2..9aaa543b 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs @@ -45,6 +45,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// public class DenseQR : QR { + /// + /// Gets or sets Tau vector. Contains additional information on Q - used for native solver. + /// + internal float[] Tau + { + get; + set; + } + /// /// Initializes a new instance of the class. This object will compute the /// QR factorization when the constructor is called and cache it's factorization. @@ -66,7 +75,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization MatrixR = matrix.Clone(); MatrixQ = new DenseMatrix(matrix.RowCount); - Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data); + Tau = new float[Math.Min(matrix.RowCount, matrix.ColumnCount)]; + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau); } /// diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 5c7a637b..93877d3d 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -69,6 +69,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The given array is the wrong length. Should be {0}.. + /// + internal static string ArgumentArrayWrongLength { + get { + return ResourceManager.GetString("ArgumentArrayWrongLength", resourceCulture); + } + } + /// /// Looks up a localized string similar to The argument must be between 0 and 1.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index f653a30b..15a681a0 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -351,4 +351,7 @@ The given array is too small. It must be at least {0} long. + + The given array is the wrong length. Should be {0}. + \ No newline at end of file diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index 7aa38802..ca7a999d 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -641,7 +641,152 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqual(b[5], 0, 14); } + [Test] + public void CanComputeQRFactorSquareMatrix() + { + var matrix = _matrices["Square3x3"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau); + + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } + + [Test] + public void CanComputeQRFactorTallMatrix() + { + var matrix = _matrices["Tall3x2"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau); + + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } + + [Test] + public void CanComputeQRFactorWideMatrix() + { + var matrix = _matrices["Wide2x3"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau); + + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } + [Test] + public void CanComputeQRFactorSquareMatrixWithWorkArray() + { + var matrix = _matrices["Square3x3"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + var work = new double[matrix.ColumnCount * Control.BlockSize]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work); + + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } + + [Test] + public void CanComputeQRFactorTallMatrixWithWorkArray() + { + var matrix = _matrices["Tall3x2"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + var work = new double[matrix.ColumnCount * Control.BlockSize]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work); + + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } + + [Test] + public void CanComputeQRFactorWideMatrixWithWorkArray() + { + var matrix = _matrices["Wide2x3"]; + var r = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, r, r.Length); + + var tau = new double[3]; + var q = new double[matrix.RowCount * matrix.RowCount]; + var work = new double[matrix.ColumnCount * Control.BlockSize]; + Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work); + + var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r); + var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); + var a = mq * mr; + + for (var row = 0; row < matrix.RowCount; row++) + { + for (var col = 0; col < matrix.ColumnCount; col++) + { + AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14); + } + } + } /// /// Checks to see if a matrix and array contain the same values.