diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 5b2bfd06..341716e1 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -85,7 +85,8 @@ ipiv dll Silverlight namespace -da +da +Dont diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 3f9e62c0..5f54a9f4 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -279,7 +279,7 @@ extern "C" { return info; } - DLLEXPORT int d_lu_solve(int n, int nrhs, double a[], double b[]) + DLLEXPORT int d_lu_solve(int n, int nrhs, double a[], double b[]) { double* clone = new double[n*n]; std::memcpy(clone, a, n*n*sizeof(double)); @@ -523,7 +523,6 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; - r[j * m + i] = 0.0f; } } } @@ -553,7 +552,6 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; - r[j * m + i] = 0.0; } } } @@ -576,10 +574,6 @@ 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) { for (int j = 0; j < m && j < n; ++j) @@ -587,7 +581,6 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; - r[j * m + i] = zero; } } } @@ -610,10 +603,6 @@ 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) @@ -621,7 +610,6 @@ extern "C" { if (i > j) { q[j * m + i] = r[j * m + i]; - r[j * m + i] = zero; } } } @@ -639,76 +627,253 @@ extern "C" { return info; } - DLLEXPORT int s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + DLLEXPORT int s_qr_solve(int m, int n, int bn, float r[], float b[], float x[], float work[], int len) { + int info = 0; + float* clone_r = new float[m*n]; + std::memcpy(clone_r, r, m*n*sizeof(float)); + + float* tau = new float[std::max(1, std::min(m,n))]; + SGEQRF(&m, &n, clone_r, &m, tau, work, &len, &info); + + if (info != 0) + { + delete[] clone_r; + delete[] tau; + return info; + } + + float* clone_b = new float[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(float)); + char side ='L'; char tran = 'T'; + SORMQR(&side, &tran, &m, &bn, &n, clone_r, &m, tau, clone_b, &m, work, &len, &info); + cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, clone_r, m, clone_b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } + } + + delete[] clone_r; + delete[] tau; + delete[] clone_b; + return info; + } + + DLLEXPORT int d_qr_solve(int m, int n, int bn, double r[], double b[], double x[], double work[], int len) + { 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); + double* clone_r = new double[m*n]; + std::memcpy(clone_r, r, m*n*sizeof(double)); + + double* tau = new double[std::max(1, std::min(m,n))]; + DGEQRF(&m, &n, clone_r, &m, tau, work, &len, &info); + + if (info != 0) + { + delete[] clone_r; + delete[] tau; + return info; + } + + double* clone_b = new double[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(double)); + + char side ='L'; + char tran = 'T'; + + DORMQR(&side, &tran, &m, &bn, &n, clone_r, &m, tau, clone_b, &m, work, &len, &info); + cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, clone_r, m, clone_b, m); for (int i = 0; i < n; ++i) { for (int j = 0; j < bn; ++j) { - x[j * n + i] = b[j * m + i]; + x[j * n + i] = clone_b[j * m + i]; } } - + + delete[] clone_b; + delete[] tau; + delete[] clone_r; + return info; + } + + DLLEXPORT int c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 x[], MKL_Complex8 work[], int len) + { + int info = 0; + MKL_Complex8* clone_r = new MKL_Complex8[m*n]; + std::memcpy(clone_r, r, m*n*sizeof(MKL_Complex8)); + + MKL_Complex8* tau = new MKL_Complex8[std::max(1, std::min(m,n))]; + CGEQRF(&m, &n, clone_r, &m, tau, work, &len, &info); + + if (info != 0) + { + delete[] clone_r; + delete[] tau; + return info; + } + + MKL_Complex8* clone_b = new MKL_Complex8[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex8)); + + char side ='L'; + char tran = 'T'; + CUNMQR(&side, &tran, &m, &bn, &n, clone_r, &m, tau, clone_b, &m, work, &len, &info); + MKL_Complex8 one; + one.real = 1.0; + cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, clone_r, m, clone_b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } + } + + delete[] clone_r; + delete[] tau; + delete[] clone_b; + return info; + } + + DLLEXPORT int z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 x[], MKL_Complex16 work[], int len) + { + int info = 0; + MKL_Complex16* clone_r = new MKL_Complex16[m*n]; + std::memcpy(clone_r, r, m*n*sizeof(MKL_Complex16)); + + MKL_Complex16* tau = new MKL_Complex16[std::max(1, std::min(m,n))]; + ZGEQRF(&m, &n, clone_r, &m, tau, work, &len, &info); + + if (info != 0) + { + delete[] clone_r; + delete[] tau; + return info; + } + + MKL_Complex16* clone_b = new MKL_Complex16[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex16)); + + char side ='L'; + char tran = 'T'; + ZUNMQR(&side, &tran, &m, &bn, &n, clone_r, &m, tau, clone_b, &m, work, &len, &info); + MKL_Complex16 one; + one.real = 1.0; + cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, clone_r, m, clone_b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } + } + + delete[] clone_r; + delete[] tau; + delete[] clone_b; return info; } - DLLEXPORT int d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + DLLEXPORT int s_qr_solve_factored(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; - 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); + + float* clone_b = new float[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(float)); + + SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); + cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = clone_b[j * m + i]; + } + } + + delete[] clone_b; + return info; + } + + DLLEXPORT int d_qr_solve_factored(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; + + double* clone_b = new double[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(double)); + + DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); + cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); for (int i = 0; i < n; ++i) { for (int j = 0; j < bn; ++j) { - x[j * n + i] = b[j * m + i]; + x[j * n + i] = clone_b[j * m + i]; } } + + delete[] clone_b; return info; } - DLLEXPORT int 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) + DLLEXPORT int c_qr_solve_factored(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* clone_b = new MKL_Complex8[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex8)); + + CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, clone_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); + cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); for (int i = 0; i < n; ++i) { for (int j = 0; j < bn; ++j) { - x[j * n + i] = b[j * m + i]; + x[j * n + i] = clone_b[j * m + i]; } } + + delete[] clone_b; return info; } - DLLEXPORT int 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) + DLLEXPORT int z_qr_solve_factored(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* clone_b = new MKL_Complex16[m*bn]; + std::memcpy(clone_b, b, m*bn*sizeof(MKL_Complex16)); + + ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, clone_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); + cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); + for (int i = 0; i < n; ++i) { for (int j = 0; j < bn; ++j) { - x[j * n + i] = b[j * m + i]; + x[j * n + i] = clone_b[j * m + i]; } } + + delete[] clone_b; return info; } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 7ff780e7..5b62552e 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -346,46 +346,65 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - void QRSolve(T[] r, int rowsR, int columnsR, T[] q, T[] b, int columnsB, T[] x); + /// Rows must be greater or equal to columns. + void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x); /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - void QRSolve(T[] r, int rowsR, int columnsR, T[] q, T[] b, int columnsB, T[] x, T[] work); + /// Rows must be greater or equal to columns. + void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, T[] work); /// /// 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 QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. /// The number of columns of B. /// On exit, the solution matrix. - void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] b, int columnsB, T[] x); + /// Rows must be greater or equal to columns. + void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x); + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + /// Rows must be greater or equal to columns. + void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x, T[] work); + /// /// Computes the singular value decomposition of A. /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 5ca428dd..bf600e9a 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -1568,128 +1568,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolve(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] b, int columnsB, Complex[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (q.Length != rowsR * rowsR) + if (b.Length != rows * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (b.Length != rowsR * columnsB) + if (x.Length != columns * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (x.Length != columnsR * columnsB) + if (rows < columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(Resources.RowsLessThanColumns); } - var work = new Complex[rowsR * rowsR]; - QRSolve(r, rowsR, columnsR, q, b, columnsB, x, work); + var work = new Complex[rows * rows]; + QRSolve(a, rows, columns, b, columnsB, x, work); } /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - public virtual void QRSolve(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] b, int columnsB, Complex[] x, Complex[] work) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, Complex[] work) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (work == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentNullException("work"); } - if (q.Length != rowsR * rowsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rowsR * columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsR * columnsB) + if (rows < columns) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (work.Length < rowsR * rowsR) + if (work.Length < rows * rows) { - work[0] = rowsR * rowsR; + work[0] = rows * rows; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - QRFactor(r, rowsR, columnsR, q, work); - QRSolveFactored(q, r, rowsR, columnsR, b, columnsB, x); + var clone = new Complex[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfComplex); + var q = new Complex[rows * rows]; + QRFactor(clone, rows, columns, q, work); + QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); - work[0] = rowsR * rowsR; + work[0] = rows * rows; + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work) + { + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x); } /// @@ -1699,10 +1715,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] b, int columnsB, Complex[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x) { if (r == null) { @@ -1731,7 +1750,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (b.Length != rowsR * columnsB) @@ -1744,6 +1763,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + var sol = new Complex[b.Length]; // Copy B matrix to "sol", so B data will not be changed diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 89e4b98f..406b6bdf 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -1568,128 +1568,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolve(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] b, int columnsB, Complex32[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (q.Length != rowsR * rowsR) + if (b.Length != rows * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (b.Length != rowsR * columnsB) + if (x.Length != columns * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (x.Length != columnsR * columnsB) + if (rows < columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(Resources.RowsLessThanColumns); } - var work = new Complex32[rowsR * rowsR]; - QRSolve(r, rowsR, columnsR, q, b, columnsB, x, work); + var work = new Complex32[rows * rows]; + QRSolve(a, rows, columns, b, columnsB, x, work); } /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - public virtual void QRSolve(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (work == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentNullException("work"); } - if (q.Length != rowsR * rowsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rowsR * columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsR * columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (work.Length < rowsR * rowsR) + if (rows < columns) { - work[0] = rowsR * rowsR; + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (work.Length < rows * rows) + { + work[0] = rows * rows; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - QRFactor(r, rowsR, columnsR, q, work); - QRSolveFactored(q, r, rowsR, columnsR, b, columnsB, x); + var clone = new Complex32[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfComplex32); + var q = new Complex32[rows * rows]; + QRFactor(clone, rows, columns, q, work); + QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); - work[0] = rowsR * rowsR; + work[0] = rows * rows; + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work) + { + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x); } /// @@ -1699,10 +1715,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] b, int columnsB, Complex32[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x) { if (r == null) { @@ -1731,7 +1750,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (b.Length != rowsR * columnsB) @@ -1744,6 +1763,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + var sol = new Complex32[b.Length]; // Copy B matrix to "sol", so B data will not be changed diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 25cc4b04..b12c56b9 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -1563,128 +1563,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolve(double[] r, int rowsR, int columnsR, double[] q, double[] b, int columnsB, double[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (q.Length != rowsR * rowsR) + if (b.Length != rows * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (b.Length != rowsR * columnsB) + if (x.Length != columns * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (x.Length != columnsR * columnsB) + if (rows < columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(Resources.RowsLessThanColumns); } - var work = new double[rowsR * rowsR]; - QRSolve(r, rowsR, columnsR, q, b, columnsB, x, work); + var work = new double[rows * rows]; + QRSolve(a, rows, columns, b, columnsB, x, work); } /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - public virtual void QRSolve(double[] r, int rowsR, int columnsR, double[] q, double[] b, int columnsB, double[] x, double[] work) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, double[] work) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (work == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentNullException("work"); } - if (q.Length != rowsR * rowsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rowsR * columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsR * columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (work.Length < rowsR * rowsR) + if (rows < columns) { - work[0] = rowsR * rowsR; + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (work.Length < rows * rows) + { + work[0] = rows * rows; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - QRFactor(r, rowsR, columnsR, q, work); - QRSolveFactored(q, r, rowsR, columnsR, b, columnsB, x); + var clone = new double[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + var q = new double[rows * rows]; + QRFactor(clone, rows, columns, q, work); + QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); - work[0] = rowsR * rowsR; + work[0] = rows * rows; + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, double[] work) + { + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x); } /// @@ -1694,10 +1710,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] b, int columnsB, double[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x) { if (r == null) { @@ -1726,7 +1745,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (b.Length != rowsR * columnsB) @@ -1739,6 +1758,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + var sol = new double[b.Length]; // Copy B matrix to "sol", so B data will not be changed diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 317beeac..e1c626a0 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -1564,130 +1564,145 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolve(float[] r, int rowsR, int columnsR, float[] q, float[] b, int columnsB, float[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (q.Length != rowsR * rowsR) + if (b.Length != rows * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (b.Length != rowsR * columnsB) + if (x.Length != columns * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (x.Length != columnsR * columnsB) + if (rows < columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(Resources.RowsLessThanColumns); } - var work = new float[rowsR * rowsR]; - QRSolve(r, rowsR, columnsR, q, b, columnsB, x, work); + var work = new float[rows * rows]; + QRSolve(a, rows, columns, b, columnsB, x, work); } /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - public virtual void QRSolve(float[] r, int rowsR, int columnsR, float[] q, float[] b, int columnsB, float[] x, float[] work) + /// Rows must be greater or equal to columns. + public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, float[] work) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) + if (a == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("a"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (r.Length != rowsR * columnsR) + if (work == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + throw new ArgumentNullException("work"); } - if (q.Length != rowsR * rowsR) + if (a.Length != rows * columns) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rowsR * columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsR * columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } - if (work.Length < rowsR * rowsR) + if (rows < columns) { - work[0] = rowsR * rowsR; + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (work.Length < rows * rows) + { + work[0] = rows * rows; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - QRFactor(r, rowsR, columnsR, q, work); - QRSolveFactored(q, r, rowsR, columnsR, b, columnsB, x); + var clone = new float[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + var q = new float[rows * rows]; + QRFactor(clone, rows, columns, q, work); + QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x); + + work[0] = rows * rows; + } - work[0] = rowsR * rowsR; + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + public virtual void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, float[] work) + { + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x); } - + /// /// Solves A*X=B for X using a previously QR factored matrix. /// @@ -1695,10 +1710,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] b, int columnsB, float[] x) + /// Rows must be greater or equal to columns. + public virtual void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x) { if (r == null) { @@ -1727,7 +1745,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (q.Length != rowsR * rowsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } if (b.Length != rowsR * columnsB) @@ -1740,6 +1758,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + var sol = new float[b.Length]; // Copy B matrix to "sol", so B data will not be changed diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index e3bc980e..ea4f4f65 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include @@ -582,7 +582,7 @@ if (work == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("work"); } if (r.Length != rowsR * columnsR) @@ -612,56 +612,259 @@ /// /// Solves A*X=B for X using 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. - /// 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. - /// On entry the B matrix; on exit the X matrix. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - [SecuritySafeCritical] - public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) + /// Rows must be greater or equal to columns. + public override void QRSolve(<#=dataType#>[] a, int rows, int columns, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (a.Length != rows * columns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != rows * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columns * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + if (rows < columns) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + var work = new <#=dataType#>[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work); } /// /// Solves A*X=B for X using 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. - /// 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. + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. /// 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. - [SecuritySafeCritical] - public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work) + /// Rows must be greater or equal to columns. + public override void QRSolve(<#=dataType#>[] a, int rows, int columns, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (work == null) + { + throw new ArgumentNullException("work"); + } + + if (a.Length != rows * columns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != rows * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columns * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + if (rows < columns) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (work.Length < 1) + { + work[0] = rows * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + SafeNativeMethods.<#=prefix#>_qr_solve(rows, columns, columnsB, a, b, x, work, work.Length); } /// /// 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. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(<#=dataType#>[] q, <#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) + public override void QRSolveFactored(<#=dataType#>[] q, <#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] tau, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { - throw new NotImplementedException(); + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (b == null) + { + throw new ArgumentNullException("q"); + } + + if (x == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rowsR * columnsR) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rowsR * rowsR) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (b.Length != rowsR * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsR * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + var work = new <#=dataType#>[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be + /// null for the native provider. The native provider uses the Q portion stored in the R matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The work array - only used in the native provider. 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. + /// Rows must be greater or equal to columns. + public override void QRSolveFactored(<#=dataType#>[] q, <#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] tau, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (b == null) + { + throw new ArgumentNullException("q"); + } + + if (x == null) + { + throw new ArgumentNullException("q"); + } + + if (work == null) + { + throw new ArgumentNullException("work"); + } + + if (r.Length != rowsR * columnsR) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rowsR * rowsR) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (b.Length != rowsR * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsR * columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + if (rowsR < columnsR) + { + throw new ArgumentException(Resources.RowsLessThanColumns); + } + + if (work.Length < 1) + { + work[0] = rowsR * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + SafeNativeMethods.<#=prefix#>_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); } /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index 59af2294..8e2695f3 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include @@ -222,4 +222,28 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> [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); + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + #endregion LAPACK diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs index 0e850de6..940e4bcc 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs @@ -167,7 +167,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index a5dd9321..0a5687a4 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs @@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs index 4df3e34e..89ec82b6 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs @@ -167,7 +167,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 403c3840..91aa662b 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs @@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs index 140bd8e8..40877b8f 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs @@ -166,7 +166,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -211,7 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index f7cfc7e4..fc05c96f 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -127,7 +127,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -172,7 +172,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs index 5a02ae42..74a23138 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs @@ -166,7 +166,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -211,7 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixQ.RowCount, MatrixQ.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index 9aaa543b..7bf7a3b0 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs @@ -127,7 +127,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -172,7 +172,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment."); } - Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, null, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 93877d3d..373ba4d3 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -699,6 +699,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The number of rows must greater than or equal to the number of columns.. + /// + internal static string RowsLessThanColumns { + get { + return ResourceManager.GetString("RowsLessThanColumns", resourceCulture); + } + } + /// /// Looks up a localized string similar to The singular vectors were not computed.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 15a681a0..4f40d7d6 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -354,4 +354,7 @@ The given array is the wrong length. Should be {0}. + + The number of rows must greater than or equal to the number of columns. + \ No newline at end of file diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index ca7a999d..e1a82d1d 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -31,7 +31,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double using Algorithms.LinearAlgebra; using LinearAlgebra.Double; - using MathNet.Numerics.LinearAlgebra.Generic; + using LinearAlgebra.Generic; using NUnit.Framework; @@ -641,6 +641,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqual(b[5], 0, 14); } + /// + /// Can compute QR factorization of a square matrix. + /// [Test] public void CanComputeQRFactorSquareMatrix() { @@ -653,7 +656,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var a = mq * mr; for (var row = 0; row < matrix.RowCount; row++) @@ -665,6 +668,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can compute QR factorization of a tall matrix. + /// [Test] public void CanComputeQRFactorTallMatrix() { @@ -676,7 +682,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); var a = mq * mr; @@ -689,6 +695,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can compute QR factorization of a wide matrix. + /// [Test] public void CanComputeQRFactorWideMatrix() { @@ -700,7 +709,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); var a = mq * mr; @@ -713,6 +722,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can compute QR factorization of a square matrix using a work array. + /// [Test] public void CanComputeQRFactorSquareMatrixWithWorkArray() { @@ -726,7 +738,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var a = mq * mr; for (var row = 0; row < matrix.RowCount; row++) @@ -738,6 +750,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can compute QR factorization of a tall matrix using a work matrix. + /// [Test] public void CanComputeQRFactorTallMatrixWithWorkArray() { @@ -750,7 +765,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); var a = mq * mr; @@ -763,6 +778,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can compute QR factorization of a wide matrix using a work matrix. + /// [Test] public void CanComputeQRFactorWideMatrixWithWorkArray() { @@ -775,7 +793,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double 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 mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q); var a = mq * mr; @@ -788,6 +806,232 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } + /// + /// Can solve Ax=b using QR factorization with a square A matrix. + /// + [Test] + public void CanSolveUsingQRSquareMatrix() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x); + + NotModified(3, 3, a, matrix); + + var mx = new DenseMatrix(matrix.ColumnCount, 2, x); + var mb = matrix * mx; + + AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14); + AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14); + AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14); + AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14); + AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14); + AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a tall A matrix. + /// + [Test] + public void CanSolveUsingQRTallMatrix() + { + var matrix = _matrices["Tall3x2"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x); + + NotModified(3, 2, a, matrix); + + var mb = new DenseMatrix(matrix.RowCount, 2, b); + var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb; + + AssertHelpers.AlmostEqual(test[0, 0], x[0], 14); + AssertHelpers.AlmostEqual(test[1, 0], x[1], 14); + AssertHelpers.AlmostEqual(test[0, 1], x[2], 14); + AssertHelpers.AlmostEqual(test[1, 1], x[3], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a square A matrix + /// using a work array. + /// + [Test] + public void CanSolveUsingQRSquareMatrixUsingWorkArray() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + var work = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work); + + NotModified(3, 3, a, matrix); + + var mx = new DenseMatrix(matrix.ColumnCount, 2, x); + var mb = matrix * mx; + + AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14); + AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14); + AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14); + AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14); + AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14); + AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a tall A matrix + /// using a work array. + /// + [Test] + public void CanSolveUsingQRTallMatrixUsingWorkArray() + { + var matrix = _matrices["Tall3x2"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + var work = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work); + + NotModified(3, 2, a, matrix); + + var mb = new DenseMatrix(matrix.RowCount, 2, b); + var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb; + + AssertHelpers.AlmostEqual(test[0, 0], x[0], 14); + AssertHelpers.AlmostEqual(test[1, 0], x[1], 14); + AssertHelpers.AlmostEqual(test[0, 1], x[2], 14); + AssertHelpers.AlmostEqual(test[1, 1], x[3], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a square A matrix + /// using a factored A matrix. + /// + [Test] + public void CanSolveUsingQRSquareMatrixOnFactoredMatrix() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.RowCount]; + Array.Copy(matrix.Data, a, a.Length); + + var tau = new double[matrix.ColumnCount]; + var q = new double[matrix.ColumnCount * matrix.ColumnCount]; + Provider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x); + + var mx = new DenseMatrix(matrix.ColumnCount, 2, x); + var mb = matrix * mx; + + AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14); + AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14); + AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14); + AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14); + AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14); + AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a tall A matrix + /// using a factored A matrix. + /// + [Test] + public void CanSolveUsingQRTallMatrixOnFactoredMatrix() + { + var matrix = _matrices["Tall3x2"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var tau = new double[matrix.ColumnCount]; + var q = new double[matrix.RowCount * matrix.RowCount]; + Provider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x); + + var mb = new DenseMatrix(matrix.RowCount, 2, b); + var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb; + + AssertHelpers.AlmostEqual(test[0, 0], x[0], 14); + AssertHelpers.AlmostEqual(test[1, 0], x[1], 14); + AssertHelpers.AlmostEqual(test[0, 1], x[2], 14); + AssertHelpers.AlmostEqual(test[1, 1], x[3], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a square A matrix + /// using a factored A matrix with a work array. + /// + [Test] + public void CanSolveUsingQRSquareMatrixOnFactoredMatrixWithWorkArray() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.RowCount]; + Array.Copy(matrix.Data, a, a.Length); + + var tau = new double[matrix.ColumnCount]; + var q = new double[matrix.ColumnCount * matrix.ColumnCount]; + var work = new double[2048]; + Provider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work); + + var mx = new DenseMatrix(matrix.ColumnCount, 2, x); + var mb = matrix * mx; + + AssertHelpers.AlmostEqual(mb[0, 0], b[0], 14); + AssertHelpers.AlmostEqual(mb[1, 0], b[1], 14); + AssertHelpers.AlmostEqual(mb[2, 0], b[2], 14); + AssertHelpers.AlmostEqual(mb[0, 1], b[3], 14); + AssertHelpers.AlmostEqual(mb[1, 1], b[4], 14); + AssertHelpers.AlmostEqual(mb[2, 1], b[5], 14); + } + + /// + /// Can solve Ax=b using QR factorization with a tall A matrix + /// using a factored A matrix with a work array. + /// + [Test] + public void CanSolveUsingQRTallMatrixOnFactoredMatrixWithWorkArray() + { + var matrix = _matrices["Tall3x2"]; + var a = new double[matrix.RowCount * matrix.ColumnCount]; + Array.Copy(matrix.Data, a, a.Length); + + var tau = new double[matrix.ColumnCount]; + var q = new double[matrix.RowCount * matrix.RowCount]; + var work = new double[2048]; + Provider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + var x = new double[matrix.ColumnCount * 2]; + Provider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work); + + var mb = new DenseMatrix(matrix.RowCount, 2, b); + var test = (matrix.Transpose() * matrix).Inverse() * matrix.Transpose() * mb; + + AssertHelpers.AlmostEqual(test[0, 0], x[0], 14); + AssertHelpers.AlmostEqual(test[1, 0], x[1], 14); + AssertHelpers.AlmostEqual(test[0, 1], x[2], 14); + AssertHelpers.AlmostEqual(test[1, 1], x[3], 14); + } + /// /// Checks to see if a matrix and array contain the same values. ///