Browse Source

native: added QR solve, tweaking the provider interface as needed

la-knuth
Marcus Cuda 16 years ago
parent
commit
6beb3b6ace
  1. 3
      src/MathNet.Numerics.5.1.ReSharper
  2. 225
      src/NativeWrappers/MKL/lapack.cpp
  3. 53
      src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs
  4. 134
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  5. 134
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  6. 134
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  7. 135
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  8. 251
      src/Numerics/Algorithms/LinearAlgebra/native.generic.include
  9. 24
      src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include
  10. 4
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs
  11. 4
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs
  12. 4
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs
  13. 4
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs
  14. 4
      src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs
  15. 4
      src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs
  16. 4
      src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs
  17. 4
      src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs
  18. 9
      src/Numerics/Properties/Resources.Designer.cs
  19. 3
      src/Numerics/Properties/Resources.resx
  20. 258
      src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

3
src/MathNet.Numerics.5.1.ReSharper

@ -85,7 +85,8 @@ ipiv
dll
Silverlight
namespace
da</UserWords>
da
Dont</UserWords>
</CustomDictionary>
</Dictionaries>
</CustomDictionaries>

225
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;
}

53
src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs

@ -346,46 +346,65 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
void QRSolve(T[] r, int rowsR, int columnsR, T[] q, T[] b, int columnsB, T[] x);
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x);
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
void QRSolve(T[] r, int rowsR, int columnsR, T[] q, T[] b, int columnsB, T[] x, T[] work);
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, T[] work);
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(T[],int,int,T[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[])"/>. </param>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[],T[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] b, int columnsB, T[] x);
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x);
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[],T[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolveFactored(T[] q, T[] r, int rowsR, int columnsR, T[] tau, T[] b, int columnsB, T[] x, T[] work);
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>

134
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@ -1568,128 +1568,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolve(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] b, int columnsB, Complex[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
public virtual void QRSolve(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] b, int columnsB, Complex[] x, Complex[] work)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
@ -1699,10 +1715,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] b, int columnsB, Complex[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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

134
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@ -1568,128 +1568,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolve(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] b, int columnsB, Complex32[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
public virtual void QRSolve(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
@ -1699,10 +1715,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] b, int columnsB, Complex32[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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

134
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@ -1563,128 +1563,144 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolve(double[] r, int rowsR, int columnsR, double[] q, double[] b, int columnsB, double[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
public virtual void QRSolve(double[] r, int rowsR, int columnsR, double[] q, double[] b, int columnsB, double[] x, double[] work)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
@ -1694,10 +1710,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] b, int columnsB, double[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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

135
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@ -1564,130 +1564,145 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolve(float[] r, int rowsR, int columnsR, float[] q, float[] b, int columnsB, float[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
public virtual void QRSolve(float[] r, int rowsR, int columnsR, float[] q, float[] b, int columnsB, float[] x, float[] work)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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;
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
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);
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
@ -1695,10 +1710,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
public virtual void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] b, int columnsB, float[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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

251
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 @@
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
[SecuritySafeCritical]
public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
[SecuritySafeCritical]
public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work)
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(<#=dataType#>[],int,int,<#=dataType#>[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(<#=dataType#>[],int,int,<#=dataType#>[])"/>. </param>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(<#=dataType#>[],int,int,<#=dataType#>[],<#=dataType#>[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">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.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
/// <summary>

24
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

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

4
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);
}
/// <summary>
@ -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);
}
}
}

9
src/Numerics/Properties/Resources.Designer.cs

@ -699,6 +699,15 @@ namespace MathNet.Numerics.Properties {
}
}
/// <summary>
/// Looks up a localized string similar to The number of rows must greater than or equal to the number of columns..
/// </summary>
internal static string RowsLessThanColumns {
get {
return ResourceManager.GetString("RowsLessThanColumns", resourceCulture);
}
}
/// <summary>
/// Looks up a localized string similar to The singular vectors were not computed..
/// </summary>

3
src/Numerics/Properties/Resources.resx

@ -354,4 +354,7 @@
<data name="ArgumentArrayWrongLength" xml:space="preserve">
<value>The given array is the wrong length. Should be {0}.</value>
</data>
<data name="RowsLessThanColumns" xml:space="preserve">
<value>The number of rows must greater than or equal to the number of columns.</value>
</data>
</root>

258
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);
}
/// <summary>
/// Can compute QR factorization of a square matrix.
/// </summary>
[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
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix.
/// </summary>
[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
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix.
/// </summary>
[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
}
}
/// <summary>
/// Can compute QR factorization of a square matrix using a work array.
/// </summary>
[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
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix using a work matrix.
/// </summary>
[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
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix using a work matrix.
/// </summary>
[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
}
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a work array.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[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);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[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);
}
/// <summary>
/// Checks to see if a matrix and array contain the same values.
/// </summary>

Loading…
Cancel
Save