mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-31 19:00:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

port the qr module to ei_solve_xxx.

Browse Source
This commit is contained in:
Benoit Jacob 2009-11-08 10:21:26 -05:00
parent 68210b03c1
commit ba7bfe110c
8 changed files with 195 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/LU/PartialPivLU.h
View File

@ -198,8 +198,6 @@ PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
compute(matrix);
}
/** This is the blocked version of ei_fullpivlu_unblocked() */
template<typename Scalar, int StorageOrder>
struct ei_partial_lu_impl

139
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -42,17 +42,17 @@
*
* \sa MatrixBase::colPivHouseholderQr()
*/
template<typename MatrixType> class ColPivHouseholderQR
template<typename _MatrixType> class ColPivHouseholderQR
{
public:
typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@ -83,22 +83,27 @@ template<typename MatrixType> class ColPivHouseholderQR
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
* *this is the QR decomposition, if any exists.
*
* \returns \c true if a solution exists, \c false if no solution exists.
*
* \param b the right-hand-side of the equation to solve.
*
* \param result a pointer to the vector/matrix in which to store the solution, if any exists.
* Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
* If no solution exists, *result is left with undefined coefficients.
* \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
*
* Example: \include ColPivHouseholderQR_solve.cpp
* Output: \verbinclude ColPivHouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
template<typename Rhs>
inline const ei_solve_return_value<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return ei_solve_return_value<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
HouseholderSequenceType matrixQ(void) const;
@ -204,36 +209,24 @@ template<typename MatrixType> class ColPivHouseholderQR
return isInjective() && isSurjective();
}
/** Computes the inverse of the matrix of which *this is the QR decomposition.
*
* \param result a pointer to the matrix into which to store the inverse. Resized if needed.
*
* \note If this matrix is not invertible, *result is left with undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*
* \sa inverse()
*/
inline void computeInverse(MatrixType *result) const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
}
/** \returns the inverse of the matrix of which *this is the QR decomposition.
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*
* \sa computeInverse()
*/
inline MatrixType inverse() const
inline const
ei_solve_return_value<ColPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return ei_solve_return_value<ColPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
(*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
}
inline int rows() const { return m_qr.rows(); }
inline int cols() const { return m_qr.cols(); }
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@ -331,50 +324,56 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool ColPivHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
template<typename MatrixType, typename Rhs, typename Dest>
struct ei_solve_impl<ColPivHouseholderQR<MatrixType>, Rhs, Dest>
: ei_solve_return_value<ColPivHouseholderQR<MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
result->resize(m_qr.cols(), b.cols());
if(m_rank==0)
void evalTo(Dest& dst) const
{
if(b.squaredNorm() == RealScalar(0))
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
const ColPivHouseholderQR<MatrixType>& dec = this->m_dec;
const Rhs& rhs = this->m_rhs;
const int rows = dec.rows(), cols = dec.cols();
dst.resize(cols, rhs.cols());
ei_assert(rhs.rows() == rows);
// FIXME introduce nonzeroPivots() and use it here. and more generally,
// make the same improvements in this dec as in FullPivLU.
if(dec.rank()==0)
{
result->setZero();
return true;
dst.setZero();
return;
}
else return false;
typename Rhs::PlainMatrixType c(rhs);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(
dec.matrixQR().corner(TopLeft,rows,dec.rank()),
dec.hCoeffs().start(dec.rank())).transpose()
);
if(!dec.isSurjective())
{
// is c is in the image of R ?
RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
// FIXME brain dead
const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
return;
}
dec.matrixQR()
.corner(TopLeft, dec.rank(), dec.rank())
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
}
const int rows = m_qr.rows();
ei_assert(b.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose());
if(!isSurjective())
{
// is c is in the image of R ?
RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
return false;
}
m_qr.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
return true;
}
};
/** \returns the matrix Q as a sequence of householder transformations */
template<typename MatrixType>

149
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -42,17 +42,17 @@
*
* \sa MatrixBase::fullPivHouseholderQr()
*/
template<typename MatrixType> class FullPivHouseholderQR
template<typename _MatrixType> class FullPivHouseholderQR
{
public:
typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@ -78,22 +78,27 @@ template<typename MatrixType> class FullPivHouseholderQR
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
* *this is the QR decomposition, if any exists.
*
* \returns \c true if a solution exists, \c false if no solution exists.
*
* \param b the right-hand-side of the equation to solve.
*
* \param result a pointer to the vector/matrix in which to store the solution, if any exists.
* Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
* If no solution exists, *result is left with undefined coefficients.
* \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
*
* Example: \include FullPivHouseholderQR_solve.cpp
* Output: \verbinclude FullPivHouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
template<typename Rhs>
inline const ei_solve_return_value<FullPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return ei_solve_return_value<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
MatrixQType matrixQ(void) const;
@ -205,36 +210,23 @@ template<typename MatrixType> class FullPivHouseholderQR
return isInjective() && isSurjective();
}
/** Computes the inverse of the matrix of which *this is the QR decomposition.
*
* \param result a pointer to the matrix into which to store the inverse. Resized if needed.
*
* \note If this matrix is not invertible, *result is left with undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*
* \sa inverse()
*/
inline void computeInverse(MatrixType *result) const
{
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
}
/** \returns the inverse of the matrix of which *this is the QR decomposition.
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*
* \sa computeInverse()
*/
inline MatrixType inverse() const
*/ inline const
ei_solve_return_value<FullPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return ei_solve_return_value<FullPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
(*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
}
inline int rows() const { return m_qr.rows(); }
inline int cols() const { return m_qr.cols(); }
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@ -340,56 +332,59 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool FullPivHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
template<typename MatrixType, typename Rhs, typename Dest>
struct ei_solve_impl<FullPivHouseholderQR<MatrixType>, Rhs, Dest>
: ei_solve_return_value<FullPivHouseholderQR<MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
result->resize(m_qr.cols(), b.cols());
if(m_rank==0)
void evalTo(Dest& dst) const
{
if(b.squaredNorm() == RealScalar(0))
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
const FullPivHouseholderQR<MatrixType>& dec = this->m_dec;
const Rhs& rhs = this->m_rhs;
const int rows = dec.rows(), cols = dec.cols();
dst.resize(cols, rhs.cols());
ei_assert(rhs.rows() == rows);
// FIXME introduce nonzeroPivots() and use it here. and more generally,
// make the same improvements in this dec as in FullPivLU.
if(dec.rank()==0)
{
result->setZero();
return true;
dst.setZero();
return;
}
else return false;
typename Rhs::PlainMatrixType c(rhs);
Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs.cols());
for (int k = 0; k < dec.rank(); ++k)
{
int remainingSize = rows-k;
c.row(k).swap(c.row(dec.rowsTranspositions().coeff(k)));
c.corner(BottomRight, remainingSize, rhs.cols())
.applyHouseholderOnTheLeft(dec.matrixQR().col(k).end(remainingSize-1),
dec.hCoeffs().coeff(k), &temp.coeffRef(0));
}
if(!dec.isSurjective())
{
// is c is in the image of R ?
RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
// FIXME brain dead
const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
return;
}
dec.matrixQR()
.corner(TopLeft, dec.rank(), dec.rank())
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
}
const int rows = m_qr.rows();
const int cols = b.cols();
ei_assert(b.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);
Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
for (int k = 0; k < m_rank; ++k)
{
int remainingSize = rows-k;
c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
c.corner(BottomRight, remainingSize, cols)
.applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
if(!isSurjective())
{
// is c is in the image of R ?
RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
return false;
}
m_qr.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
return true;
}
};
/** \returns the matrix Q */
template<typename MatrixType>

72
Eigen/src/QR/HouseholderQR.h
View File

@ -46,17 +46,17 @@
*
* \sa MatrixBase::householderQr()
*/
template<typename MatrixType> class HouseholderQR
template<typename _MatrixType> class HouseholderQR
{
public:
typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
@ -85,19 +85,26 @@ template<typename MatrixType> class HouseholderQR
*
* \param b the right-hand-side of the equation to solve.
*
* \param result a pointer to the vector/matrix in which to store the solution, if any exists.
* Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
* If no solution exists, *result is left with undefined coefficients.
* \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
*
* Example: \include HouseholderQR_solve.cpp
* Output: \verbinclude HouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
template<typename Rhs>
inline const ei_solve_return_value<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return ei_solve_return_value<HouseholderQR, Rhs>(*this, b.derived());
}
MatrixQType matrixQ() const;
HouseholderSequenceType matrixQAsHouseholderSequence() const
@ -145,6 +152,10 @@ template<typename MatrixType> class HouseholderQR
*/
typename MatrixType::RealScalar logAbsDeterminant() const;
inline int rows() const { return m_qr.rows(); }
inline int cols() const { return m_qr.cols(); }
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@ -198,31 +209,36 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
void HouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
template<typename MatrixType, typename Rhs, typename Dest>
struct ei_solve_impl<HouseholderQR<MatrixType>, Rhs, Dest>
: ei_solve_return_value<HouseholderQR<MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
result->derived().resize(m_qr.cols(), b.cols());
const int rows = m_qr.rows();
const int rank = std::min(m_qr.rows(), m_qr.cols());
ei_assert(b.rows() == rows);
void evalTo(Dest& dst) const
{
const HouseholderQR<MatrixType>& dec = this->m_dec;
const Rhs& rhs = this->m_rhs;
const int rows = dec.rows(), cols = dec.cols();
dst.resize(cols, rhs.cols());
const int rank = std::min(rows, cols);
ei_assert(rhs.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);
typename Rhs::PlainMatrixType c(rhs);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(
dec.matrixQR().corner(TopLeft,rows,rank),
dec.hCoeffs().start(rank)).transpose()
);
m_qr.corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, rank, c.cols()));
dec.matrixQR()
.corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, rank, c.cols()));
result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
}
dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols());
dst.corner(BottomLeft, cols-rank, c.cols()).setZero();
}
};
/** \returns the matrix Q */
template<typename MatrixType>

1
Eigen/src/SVD/SVD.h
View File

@ -86,7 +86,6 @@ template<typename _MatrixType> class SVD
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
* \note_about_using_kernel_to_study_multiple_solutions
*
* \sa MatrixBase::svd(),
*/

6
test/qr.cpp
View File

@ -63,7 +63,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
qr.solve(m3, &m2);
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
}
@ -86,7 +86,7 @@ template<typename MatrixType> void qr_invertible()
HouseholderQR<MatrixType> qr(m1);
m3 = MatrixType::Random(size,size);
qr.solve(m3, &m2);
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
@ -106,7 +106,7 @@ template<typename MatrixType> void qr_verify_assert()
HouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
VERIFY_RAISES_ASSERT(qr.absDeterminant())
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())

13
test/qr_colpivoting.cpp
View File

@ -61,10 +61,8 @@ template<typename MatrixType> void qr()
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
VERIFY(qr.solve(m3, &m2));
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
m3 = MatrixType::Random(rows,cols2);
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType, int Cols2> void qr_fixedsize()
@ -95,10 +93,8 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
VERIFY(qr.solve(m3, &m2));
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
m3 = Matrix<Scalar,Rows,Cols2>::Random(Rows,Cols2);
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType> void qr_invertible()
@ -120,7 +116,7 @@ template<typename MatrixType> void qr_invertible()
ColPivHouseholderQR<MatrixType> qr(m1);
m3 = MatrixType::Random(size,size);
qr.solve(m3, &m2);
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
@ -140,13 +136,12 @@ template<typename MatrixType> void qr_verify_assert()
ColPivHouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
VERIFY_RAISES_ASSERT(qr.isInjective())
VERIFY_RAISES_ASSERT(qr.isSurjective())
VERIFY_RAISES_ASSERT(qr.isInvertible())
VERIFY_RAISES_ASSERT(qr.computeInverse(&tmp))
VERIFY_RAISES_ASSERT(qr.inverse())
VERIFY_RAISES_ASSERT(qr.absDeterminant())
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())

9
test/qr_fullpivoting.cpp
View File

@ -61,10 +61,8 @@ template<typename MatrixType> void qr()
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
VERIFY(qr.solve(m3, &m2));
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
m3 = MatrixType::Random(rows,cols2);
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType> void qr_invertible()
@ -90,7 +88,7 @@ template<typename MatrixType> void qr_invertible()
VERIFY(qr.isSurjective());
m3 = MatrixType::Random(size,size);
VERIFY(qr.solve(m3, &m2));
m2 = qr.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
@ -110,13 +108,12 @@ template<typename MatrixType> void qr_verify_assert()
FullPivHouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
VERIFY_RAISES_ASSERT(qr.isInjective())
VERIFY_RAISES_ASSERT(qr.isSurjective())
VERIFY_RAISES_ASSERT(qr.isInvertible())
VERIFY_RAISES_ASSERT(qr.computeInverse(&tmp))
VERIFY_RAISES_ASSERT(qr.inverse())
VERIFY_RAISES_ASSERT(qr.absDeterminant())
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())