GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-15 09:18:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge with head

Browse Source
This commit is contained in:
Thomas Capricelli 2009-08-24 08:55:27 +02:00
commit dff5135026
29 changed files with 1270 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/QR
View File

@ -36,6 +36,8 @@ namespace Eigen {
*/
#include "src/QR/QR.h"
#include "src/QR/FullPivotingHouseholderQR.h"
#include "src/QR/ColPivotingHouseholderQR.h"
#include "src/QR/Tridiagonalization.h"
#include "src/QR/EigenSolver.h"
#include "src/QR/SelfAdjointEigenSolver.h"

9
Eigen/StdVector
View File

@ -136,10 +136,8 @@ class vector<T,Eigen::aligned_allocator<T> >
{ return vector_base::insert(position,x); }
void insert(const_iterator position, size_type new_size, const value_type& x)
{ vector_base::insert(position, new_size, x); }
#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,1)
#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2)
// workaround GCC std::vector implementation
// Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
// no no need to workaround !
void resize(size_type new_size, const value_type& x)
{
if (new_size < vector_base::size())
@ -147,9 +145,12 @@ class vector<T,Eigen::aligned_allocator<T> >
else
vector_base::insert(vector_base::end(), new_size - vector_base::size(), x);
}
#elif defined(_GLIBCXX_VECTOR)
#elif defined(_GLIBCXX_VECTOR) && (!EIGEN_GNUC_AT_LEAST(4,1))
// Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
// no no need to workaround !
using vector_base::resize;
#else
// either GCC 4.1 or non-GCC
// default implementation which should always work.
void resize(size_type new_size, const value_type& x)
{

2
Eigen/src/Core/MatrixBase.h
View File

@ -745,6 +745,8 @@ template<typename Derived> class MatrixBase
/////////// QR module ///////////
const HouseholderQR<PlainMatrixType> householderQr() const;
const ColPivotingHouseholderQR<PlainMatrixType> colPivotingHouseholderQr() const;
const FullPivotingHouseholderQR<PlainMatrixType> fullPivotingHouseholderQr() const;
EigenvaluesReturnType eigenvalues() const;
RealScalar operatorNorm() const;

21
Eigen/src/Core/Product.h
View File

@ -61,11 +61,22 @@ template<typename Lhs, typename Rhs> struct ei_product_type
enum {
Rows = Lhs::RowsAtCompileTime,
Cols = Rhs::ColsAtCompileTime,
Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime),
value = ei_product_type_selector<(Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small)),
(Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small)),
(Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small))>::ret
Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime)
};
// the splitting into different lines of code here, introducing the _select enums and the typedef below,
// is to work around an internal compiler error with gcc 4.1 and 4.2.
private:
enum {
rows_select = Rows >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Rows==1 ? 1 : Small),
cols_select = Cols >=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Cols==1 ? 1 : Small),
depth_select = Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small)
};
typedef ei_product_type_selector<rows_select, cols_select, depth_select> product_type_selector;
public:
enum {
value = product_type_selector::ret
};
};

2
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -117,6 +117,8 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse;
template<typename MatrixType> class LU;
template<typename MatrixType> class PartialLU;
template<typename MatrixType> class HouseholderQR;
template<typename MatrixType> class ColPivotingHouseholderQR;
template<typename MatrixType> class FullPivotingHouseholderQR;
template<typename MatrixType> class SVD;
template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
template<typename MatrixType> class LDLT;

5
Eigen/src/Householder/Householder.h
View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -73,10 +74,10 @@ void MatrixBase<Derived>::makeHouseholder(
EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
RealScalar tailSqNorm;
RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm();
Scalar c0 = coeff(0);
if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0))
if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0))
{
*tau = 0;
*beta = ei_real(c0);

38
Eigen/src/Jacobi/Jacobi.h
View File

@ -26,17 +26,32 @@
#ifndef EIGEN_JACOBI_H
#define EIGEN_JACOBI_H
/** Applies the counter clock wise 2D rotation of angle \c theta given by its
* cosine \a c and sine \a s to the set of 2D vectors of cordinates \a x and \a y:
* \f$ x = c x - s' y \f$
* \f$ y = s x + c y \f$
*
* \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight()
*/
template<typename VectorX, typename VectorY>
void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s);
/** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this.
* More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ]
* \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
inline void MatrixBase<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s)
{
RowXpr x(row(p));
RowXpr y(row(q));
ei_apply_rotation_in_the_plane(x, y, c, s);
ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s));
}
/** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this.
* More precisely, it computes B = B * J, with J = [c s ; -s' c] and B = [ *this.col(p) ; *this.col(q) ]
* \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane()
*/
template<typename Derived>
inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s)
{
@ -45,6 +60,12 @@ inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, S
ei_apply_rotation_in_the_plane(x, y, c, s);
}
/** Computes the cosine-sine pair (\a c, \a s) such that its associated
* rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$
* applied to both the right and left of the 2x2 matrix
* \f$ B = ( \begin{array}{cc} x & y \\ _ & z \end{array} )\f$ yields
* a diagonal matrix A: \f$ A = J' B J \f$
*/
template<typename Scalar>
bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s)
{
@ -128,12 +149,13 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
const Packet pc = ei_pset1(c);
const Packet ps = ei_pset1(s);
ei_conj_helper<true,false> cj;
for(int i=0; i<alignedStart; ++i)
{
Scalar xi = x[i];
Scalar yi = y[i];
x[i] = c * xi - s * yi;
x[i] = c * xi - ei_conj(s) * yi;
y[i] = s * xi + c * yi;
}
@ -146,7 +168,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Packet xi = ei_pload(px);
Packet yi = ei_pload(py);
ei_pstore(px, ei_psub(ei_pmul(pc,xi),ei_pmul(ps,yi)));
ei_pstore(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore(py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
px += PacketSize;
py += PacketSize;
@ -161,8 +183,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
Packet xi1 = ei_ploadu(px+PacketSize);
Packet yi = ei_pload (py);
Packet yi1 = ei_pload (py+PacketSize);
ei_pstoreu(px, ei_psub(ei_pmul(pc,xi),ei_pmul(ps,yi)));
ei_pstoreu(px+PacketSize, ei_psub(ei_pmul(pc,xi1),ei_pmul(ps,yi1)));
ei_pstoreu(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstoreu(px+PacketSize, ei_psub(ei_pmul(pc,xi1),cj.pmul(ps,yi1)));
ei_pstore (py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
ei_pstore (py+PacketSize, ei_padd(ei_pmul(ps,xi1),ei_pmul(pc,yi1)));
px += Peeling*PacketSize;
@ -172,7 +194,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Packet xi = ei_ploadu(x+peelingEnd);
Packet yi = ei_pload (y+peelingEnd);
ei_pstoreu(x+peelingEnd, ei_psub(ei_pmul(pc,xi),ei_pmul(ps,yi)));
ei_pstoreu(x+peelingEnd, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi)));
ei_pstore (y+peelingEnd, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi)));
}
}
@ -181,7 +203,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Scalar xi = x[i];
Scalar yi = y[i];
x[i] = c * xi - s * yi;
x[i] = c * xi - ei_conj(s) * yi;
y[i] = s * xi + c * yi;
}
}
@ -191,7 +213,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY&
{
Scalar xi = *x;
Scalar yi = *y;
*x = c * xi - s * yi;
*x = c * xi - ei_conj(s) * yi;
*y = s * xi + c * yi;
x += incrx;
y += incry;

26
Eigen/src/LU/LU.h
View File

@ -42,11 +42,10 @@
* This decomposition provides the generic approach to solving systems of linear equations, computing
* the rank, invertibility, inverse, kernel, and determinant.
*
* This LU decomposition is very stable and well tested with large matrices. Even exact rank computation
* works at sizes larger than 1000x1000. However there are use cases where the SVD decomposition is inherently
* more stable when dealing with numerically damaged input. For example, computing the kernel is more stable with
* SVD because the SVD can determine which singular values are negligible while LU has to work at the level of matrix
* coefficients that are less meaningful in this respect.
* This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD
* decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix,
* working with the SVD allows to select the smallest singular values of the matrix, something that
* the LU decomposition doesn't see.
*
* The data of the LU decomposition can be directly accessed through the methods matrixLU(),
* permutationP(), permutationQ().
@ -326,6 +325,7 @@ template<typename MatrixType> class LU
inline void computeInverse(MatrixType *result) const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
}
@ -456,6 +456,7 @@ template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
return Scalar(m_det_pq) * m_lu.diagonal().prod();
}
@ -533,7 +534,8 @@ bool LU<MatrixType>::solve(
) const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
if(m_rank==0) return false;
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
* So we proceed as follows:
* Step 1: compute c = Pb.
@ -552,7 +554,8 @@ bool LU<MatrixType>::solve(
for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i);
// Step 2
m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template triangularView<UnitLowerTriangular>()
m_lu.corner(Eigen::TopLeft,smalldim,smalldim)
.template triangularView<UnitLowerTriangular>()
.solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols)
{
@ -564,11 +567,10 @@ bool LU<MatrixType>::solve(
if(!isSurjective())
{
// is c is in the image of U ?
RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : 0;
for(int col = 0; col < c.cols(); ++col)
for(int row = m_rank; row < c.rows(); ++row)
if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
return false;
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_lu.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>()

266
Eigen/src/QR/ColPivotingHouseholderQR.h Normal file
View File

@ -0,0 +1,266 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
/** \ingroup QR_Module
* \nonstableyet
*
* \class ColPivotingHouseholderQR
*
* \brief Householder rank-revealing QR decomposition of a matrix
*
* \param MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a rank-revealing QR decomposition using Householder transformations.
*
* This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal
* numerical stability.
*
* \sa MatrixBase::colPivotingHouseholderQr()
*/
template<typename MatrixType> class ColPivotingHouseholderQR
{
public:
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;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ColPivotingHouseholderQR::compute(const MatrixType&).
*/
ColPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
ColPivotingHouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
m_isInitialized(false)
{
compute(matrix);
}
/** 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.
*
* \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.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* Example: \include ColPivotingHouseholderQR_solve.cpp
* Output: \verbinclude ColPivotingHouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
const MatrixType& matrixQR() const { return m_qr; }
ColPivotingHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
return m_cols_permutation;
}
inline int rank() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
return m_rank;
}
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
IntRowVectorType m_cols_permutation;
bool m_isInitialized;
RealScalar m_precision;
int m_rank;
int m_det_pq;
};
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
int rows = matrix.rows();
int cols = matrix.cols();
int size = std::min(rows,cols);
m_rank = size;
m_qr = matrix;
m_hCoeffs.resize(size);
RowVectorType temp(cols);
m_precision = epsilon<Scalar>() * size;
IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
RealRowVectorType colSqNorms(cols);
for(int k = 0; k < cols; ++k)
colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
RealScalar biggestColSqNorm = colSqNorms.maxCoeff();
for (int k = 0; k < size; ++k)
{
int biggest_col_in_corner;
RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner);
biggest_col_in_corner += k;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision))
{
m_rank = k;
for(int i = k; i < size; i++)
{
cols_transpositions.coeffRef(i) = i;
m_hCoeffs.coeffRef(i) = Scalar(0);
}
break;
}
cols_transpositions.coeffRef(k) = biggest_col_in_corner;
if(k != biggest_col_in_corner) {
m_qr.col(k).swap(m_qr.col(biggest_col_in_corner));
++number_of_transpositions;
}
RealScalar beta;
m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta);
m_qr.coeffRef(k,k) = beta;
m_qr.corner(BottomRight, rows-k, cols-k-1)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2();
}
for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k;
for(int k = 0; k < size; ++k)
std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
void ColPivotingHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
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.corner(BottomRight, remainingSize, cols)
.applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
m_qr.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
result->resize(m_qr.cols(), b.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();
}
/** \returns the matrix Q */
template<typename MatrixType>
MatrixType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
// compute the product H'_0 H'_1 ... H'_n-1,
// where H_k is the k-th Householder transformation I - h_k v_k v_k'
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
int rows = m_qr.rows();
int cols = m_qr.cols();
int size = std::min(rows,cols);
MatrixType res = MatrixType::Identity(rows, rows);
Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
for (int k = size-1; k >= 0; k--)
{
res.block(k, k, rows-k, rows-k)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
}
return res;
}
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the column-pivoting Householder QR decomposition of \c *this.
*
* \sa class ColPivotingHouseholderQR
*/
template<typename Derived>
const ColPivotingHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::colPivotingHouseholderQr() const
{
return ColPivotingHouseholderQR<PlainMatrixType>(eval());
}
#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H

395
Eigen/src/QR/FullPivotingHouseholderQR.h Normal file
View File

@ -0,0 +1,395 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
/** \ingroup QR_Module
* \nonstableyet
*
* \class FullPivotingHouseholderQR
*
* \brief Householder rank-revealing QR decomposition of a matrix
*
* \param MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a rank-revealing QR decomposition using Householder transformations.
*
* This decomposition performs full-pivoting in order to be rank-revealing and achieve optimal
* numerical stability.
*
* \sa MatrixBase::householderRrqr()
*/
template<typename MatrixType> class FullPivotingHouseholderQR
{
public:
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;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&).
*/
FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
FullPivotingHouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
m_isInitialized(false)
{
compute(matrix);
}
/** 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.
*
* \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.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* Example: \include FullPivotingHouseholderQR_solve.cpp
* Output: \verbinclude FullPivotingHouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
const MatrixType& matrixQR() const { return m_qr; }
FullPivotingHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_cols_permutation;
}
const IntColVectorType& rowsTranspositions() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_rows_transpositions;
}
/** \returns the absolute value of the determinant of the matrix of which
* *this is the QR decomposition. It has only linear complexity
* (that is, O(n) where n is the dimension of the square matrix)
* as the QR decomposition has already been computed.
*
* \note This is only for square matrices.
*
* \warning a determinant can be very big or small, so for matrices
* of large enough dimension, there is a risk of overflow/underflow.
*
* \sa MatrixBase::determinant()
*/
typename MatrixType::RealScalar absDeterminant() const;
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note This is computed at the time of the construction of the QR decomposition. This
* method does not perform any further computation.
*/
inline int rank() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_rank;
}
/** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
*
* \note Since the rank is computed at the time of the construction of the QR decomposition, this
* method almost does not perform any further computation.
*/
inline int dimensionOfKernel() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_qr.cols() - m_rank;
}
/** \returns true if the matrix of which *this is the QR decomposition represents an injective
* linear map, i.e. has trivial kernel; false otherwise.
*
* \note Since the rank is computed at the time of the construction of the QR decomposition, this
* method almost does not perform any further computation.
*/
inline bool isInjective() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_rank == m_qr.cols();
}
/** \returns true if the matrix of which *this is the QR decomposition represents a surjective
* linear map; false otherwise.
*
* \note Since the rank is computed at the time of the construction of the QR decomposition, this
* method almost does not perform any further computation.
*/
inline bool isSurjective() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_rank == m_qr.rows();
}
/** \returns true if the matrix of which *this is the QR decomposition is invertible.
*
* \note Since the rank is computed at the time of the construction of the QR decomposition, this
* method almost does not perform any further computation.
*/
inline bool isInvertible() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
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 && "FullPivotingHouseholderQR 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
{
MatrixType result;
computeInverse(&result);
return result;
}
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
IntColVectorType m_rows_transpositions;
IntRowVectorType m_cols_permutation;
bool m_isInitialized;
RealScalar m_precision;
int m_rank;
int m_det_pq;
};
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::absDeterminant() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
return ei_abs(m_qr.diagonal().prod());
}
template<typename MatrixType>
FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
int rows = matrix.rows();
int cols = matrix.cols();
int size = std::min(rows,cols);
m_rank = size;
m_qr = matrix;
m_hCoeffs.resize(size);
RowVectorType temp(cols);
m_precision = epsilon<Scalar>() * size;
m_rows_transpositions.resize(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
RealScalar biggest;
for (int k = 0; k < size; ++k)
{
int row_of_biggest_in_corner, col_of_biggest_in_corner;
RealScalar biggest_in_corner;
biggest_in_corner = m_qr.corner(Eigen::BottomRight, rows-k, cols-k)
.cwise().abs()
.maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
row_of_biggest_in_corner += k;
col_of_biggest_in_corner += k;
if(k==0) biggest = biggest_in_corner;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
{
m_rank = k;
for(int i = k; i < size; i++)
{
m_rows_transpositions.coeffRef(i) = i;
cols_transpositions.coeffRef(i) = i;
m_hCoeffs.coeffRef(i) = Scalar(0);
}
break;
}
m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) {
m_qr.row(k).end(cols-k).swap(m_qr.row(row_of_biggest_in_corner).end(cols-k));
++number_of_transpositions;
}
if(k != col_of_biggest_in_corner) {
m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
++number_of_transpositions;
}
RealScalar beta;
m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta);
m_qr.coeffRef(k,k) = beta;
m_qr.corner(BottomRight, rows-k, cols-k-1)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
}
for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k;
for(int k = 0; k < size; ++k)
std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool FullPivotingHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
if(m_rank==0) return false;
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()));
result->resize(m_qr.cols(), b.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>
MatrixType FullPivotingHouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
// compute the product H'_0 H'_1 ... H'_n-1,
// where H_k is the k-th Householder transformation I - h_k v_k v_k'
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
int rows = m_qr.rows();
int cols = m_qr.cols();
int size = std::min(rows,cols);
MatrixType res = MatrixType::Identity(rows, rows);
Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
for (int k = size-1; k >= 0; k--)
{
res.block(k, k, rows-k, rows-k)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
}
return res;
}
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the full-pivoting Householder QR decomposition of \c *this.
*
* \sa class FullPivotingHouseholderQR
*/
template<typename Derived>
const FullPivotingHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::fullPivotingHouseholderQr() const
{
return FullPivotingHouseholderQR<PlainMatrixType>(eval());
}
#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H

9
Eigen/src/QR/QR.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -39,7 +39,7 @@
*
* Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
*
* \sa MatrixBase::qr()
* \sa MatrixBase::householderQr()
*/
template<typename MatrixType> class HouseholderQR
{
@ -54,6 +54,7 @@ template<typename MatrixType> class HouseholderQR
typedef Block<MatrixType, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixRBlockType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
typedef Matrix<Scalar, MinSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
/**
* \brief Default Constructor.
@ -125,12 +126,12 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
m_qr = matrix;
m_hCoeffs.resize(size);
Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
RowVectorType temp(cols);
for (int k = 0; k < size; ++k)
{
int remainingRows = rows - k;
int remainingCols = cols - k -1;
int remainingCols = cols - k - 1;
RealScalar beta;
m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta);

2
cmake/EigenTesting.cmake
View File

@ -74,7 +74,7 @@ macro(ei_add_test testname)
string(STRIP "${ARGV2}" ARGV2_stripped)
string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length)
if(${ARGV2_stripped_length} GREATER 0)
target_link_libraries(${targetname} "${ARGV2}")
target_link_libraries(${targetname} ${ARGV2})
endif(${ARGV2_stripped_length} GREATER 0)
endif(${ARGC} GREATER 2)

2
doc/C01_QuickStartGuide.dox
View File

@ -1,6 +1,6 @@
namespace Eigen {
/** \page TutorialCore Tutorial 1/3 - Core features
/** \page TutorialCore Tutorial 1/4 - Core features
\ingroup Tutorial
<div class="eimainmenu">\ref index "Overview"

2
doc/C03_TutorialGeometry.dox
View File

@ -1,6 +1,6 @@
namespace Eigen {
/** \page TutorialGeometry Tutorial 2/3 - Geometry
/** \page TutorialGeometry Tutorial 2/4 - Geometry
\ingroup Tutorial
<div class="eimainmenu">\ref index "Overview"

262
doc/C05_TutorialLinearAlgebra.dox
View File

@ -1,7 +1,6 @@
namespace Eigen {
/** \page TutorialAdvancedLinearAlgebra Tutorial 3/3 - Advanced linear algebra
/** \page TutorialAdvancedLinearAlgebra Tutorial 3/4 - Advanced linear algebra
\ingroup Tutorial
<div class="eimainmenu">\ref index "Overview"
@ -11,6 +10,9 @@ namespace Eigen {
| \ref TutorialSparse "Sparse matrix"
</div>
This tutorial chapter explains how you can use Eigen to tackle various problems involving matrices:
solving systems of linear equations, finding eigenvalues and eigenvectors, and so on.
\b Table \b of \b contents
- \ref TutorialAdvSolvers
- \ref TutorialAdvLU
@ -18,53 +20,129 @@ namespace Eigen {
- \ref TutorialAdvQR
- \ref TutorialAdvEigenProblems
\section TutorialAdvSolvers Solving linear problems
This part of the tutorial focuses on solving linear problem of the form \f$ A \mathbf{x} = b \f$,
where both \f$ A \f$ and \f$ b \f$ are known, and \f$ x \f$ is the unknown. Moreover, \f$ A \f$
assumed to be a squared matrix. Of course, the methods described here can be used whenever an expression
involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} B \f$.
Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of
the matrix \f$ A \f$, such as its shape, size and numerical properties.
This part of the tutorial focuses on solving systems of linear equations. Such statems can be
written in the form \f$ A \mathbf{x} = \mathbf{b} \f$, where both \f$ A \f$ and \f$ \mathbf{b} \f$
are known, and \f$ \mathbf{x} \f$ is the unknown. Moreover, \f$ A \f$ is assumed to be a square
matrix.
The equation \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution if \f$ A \f$ is invertible. If
the matrix is not invertible, then the equation may have no or infinitely many solutions. All
solvers assume that \f$ A \f$ is invertible, unless noted otherwise.
Eigen offers various algorithms to this problem. The choice of algorithm mainly depends on the
nature of the matrix \f$ A \f$, such as its shape, size and numerical properties.
- The \ref TutorialAdvSolvers_LU "LU decomposition" (with partial pivoting) is a general-purpose
algorithm which works for most problems.
- Use the \ref TutorialAdvSolvers_Cholesky "Cholesky decomposition" if the matrix \f$ A \f$ is
positive definite.
- Use a special \ref TutorialAdvSolvers_Triangular "triangular solver" if the matrix \f$ A \f$ is
upper or lower triangular.
- Use of the \ref TutorialAdvSolvers_Inverse "matrix inverse" is not recommended in general, but
may be appropriate in special cases, for instance if you want to solve several systems with the
same matrix \f$ A \f$ and that matrix is small.
- \ref TutorialAdvSolvers_Misc "Other solvers" (%LU decomposition with full pivoting, the singular
value decomposition) are provided for special cases, such as when \f$ A \f$ is not invertible.
The methods described here can be used whenever an expression involve the product of an inverse
matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1} B \f$.
\subsection TutorialAdvSolvers_LU LU decomposition (with partial pivoting)
This is a general-purpose algorithm which performs well in most cases (provided the matrix \f$ A \f$
is invertible), so if you are unsure about which algorithm to pick, choose this. The method proceeds
in two steps. First, the %LU decomposition with partial pivoting is computed using the
MatrixBase::partialLu() function. This yields an object of the class PartialLU. Then, the
PartialLU::solve() method is called to compute a solution.
As an example, suppose we want to solve the following system of linear equations:
\f[ \begin{aligned}
x + 2y + 3z &= 3 \\
4x + 5y + 6z &= 3 \\
7x + 8y + 10z &= 4.
\end{aligned} \f]
The following program solves this system:
\subsection TutorialAdvSolvers_Triangular Triangular solver
If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the diagonal
are all not zero), then the problem can be solved directly using MatrixBase::solveTriangular(), or better,
MatrixBase::solveTriangularInPlace(). Here is an example:
<table class="tutorial_code"><tr><td>
\include MatrixBase_marked.cpp
</td>
<td>
output:
\include MatrixBase_marked.out
\include Tutorial_PartialLU_solve.cpp
</td><td>
output: \include Tutorial_PartialLU_solve.out
</td></tr></table>
See MatrixBase::solveTriangular() for more details.
There are many situations in which we want to solve the same system of equations with different
right-hand sides. One possibility is to put the right-hand sides as columns in a matrix \f$ B \f$
and then solve the equation \f$ A X = B \f$. For instance, suppose that we want to solve the same
system as before, but now we also need the solution of the same equations with 1 on the right-hand
side. The following code computes the required solutions:
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_multiple_rhs.cpp
</td><td>
output: \include Tutorial_solve_multiple_rhs.out
</td></tr></table>
However, this is not always possible. Often, you only know the right-hand side of the second
problem, and whether you want to solve it at all, after you solved the first problem. In such a
case, it's best to save the %LU decomposition and reuse it to solve the second problem. This is
worth the effort because computing the %LU decomposition is much more expensive than using it to
solve the equation. Here is some code to illustrate the procedure. It uses the constructor
PartialLU::PartialLU(const MatrixType&) to compute the %LU decomposition.
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_reuse_decomposition.cpp
</td><td>
output: \include Tutorial_solve_reuse_decomposition.out
</td></tr></table>
\b Warning: All this code presumes that the matrix \f$ A \f$ is invertible, so that the system
\f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution. If the matrix \f$ A \f$ is not invertible,
then the system \f$ A \mathbf{x} = \mathbf{b} \f$ has either zero or infinitely many solutions. In
both cases, PartialLU::solve() will give nonsense results. For example, suppose that we want to
solve the same system as above, but with the 10 in the last equation replaced by 9. Then the system
of equations is inconsistent: adding the first and the third equation gives \f$ 8x + 10y + 12z = 7 \f$,
which implies \f$ 4x + 5y + 6z = 3\frac12 \f$, in contradiction with the seocond equation. If we try
to solve this inconsistent system with Eigen, we find:
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_singular.cpp
</td><td>
output: \include Tutorial_solve_singular.out
</td></tr></table>
The %LU decomposition with \b full pivoting (class LU) and the singular value decomposition (class
SVD) may be helpful in this case, as explained in the section \ref TutorialAdvSolvers_Misc below.
\sa LU_Module, MatrixBase::partialLu(), PartialLU::solve(), class PartialLU.
\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices)
If the matrix \f$ A \f$ is small (\f$ \leq 4 \f$) and invertible, then a good approach is to directly compute
the inverse of the matrix \f$ A \f$, and then obtain the solution \f$ x \f$ by \f$ \mathbf{x} = A^{-1} b \f$. With Eigen,
this can be implemented like this:
\subsection TutorialAdvSolvers_Cholesky Cholesky decomposition
\code
#include <Eigen/LU>
Matrix4f A = Matrix4f::Random();
Vector4f b = Vector4f::Random();
Vector4f x = A.inverse() * b;
\endcode
If the matrix \f$ A \f$ is \b symmetric \b positive \b definite, then the best method is to use a
Cholesky decomposition: it is both faster and more accurate than the %LU decomposition. Such
positive definite matrices often arise when solving overdetermined problems. These are linear
systems \f$ A \mathbf{x} = \mathbf{b} \f$ in which the matrix \f$ A \f$ has more rows than columns,
so that there are more equations than unknowns. Typically, there is no vector \f$ \mathbf{x} \f$
which satisfies all the equation. Instead, we look for the least-square solution, that is, the
vector \f$ \mathbf{x} \f$ for which \f$ \| A \mathbf{x} - \mathbf{b} \|_2 \f$ is minimal. You can
find this vector by solving the equation \f$ A^T \! A \mathbf{x} = A^T \mathbf{b} \f$. If the matrix
\f$ A \f$ has full rank, then \f$ A^T \! A \f$ is positive definite and thus you can use the
Cholesky decomposition to solve this system and find the least-square solution to the original
system \f$ A \mathbf{x} = \mathbf{b} \f$.
Note that the function inverse() is defined in the LU module.
See MatrixBase::inverse() for more details.
Eigen offers two different Cholesky decompositions: the LLT class provides a \f$ LL^T \f$
decomposition where L is a lower triangular matrix, and the LDLT class provides a \f$ LDL^T \f$
decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. The latter
includes pivoting and avoids square roots; this makes the %LDLT decomposition slightly more stable
than the %LLT decomposition. The LDLT class is able to handle both positive- and negative-definite
matrices, but not indefinite matrices.
The API is the same as when using the %LU decomposition.
\subsection TutorialAdvSolvers_Symmetric Cholesky (for positive definite matrices)
If the matrix \f$ A \f$ is \b positive \b definite, then
the best method is to use a Cholesky decomposition.
Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below).
Eigen offers two different Cholesky decompositions: a \f$ LL^T \f$ decomposition where L is a lower triangular matrix,
and a \f$ LDL^T \f$ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix.
The latter avoids square roots and is therefore slightly more stable than the former one.
\code
#include <Eigen/Cholesky>
MatrixXf D = MatrixXf::Random(8,4);
@ -74,15 +152,19 @@ VectorXf x;
A.llt().solve(b,&x); // using a LLT factorization
A.ldlt().solve(b,&x); // using a LDLT factorization
\endcode
when the value of the right hand side \f$ b \f$ is not needed anymore, then it is faster to use
the \em in \em place API, e.g.:
The LLT and LDLT classes also provide an \em in \em place API for the case where the value of the
right hand-side \f$ b \f$ is not needed anymore.
\code
A.llt().solveInPlace(b);
\endcode
In that case the value of \f$ b \f$ is replaced by the result \f$ x \f$.
If the linear problem has to solved for various right hand side \f$ b_i \f$, but same matrix \f$ A \f$,
then it is highly recommended to perform the decomposition of \$ A \$ only once, e.g.:
This code replaces the vector \f$ b \f$ by the result \f$ x \f$.
As before, you can reuse the factorization if you have to solve the same linear problem with
different right-hand sides, e.g.:
\code
// ...
LLT<MatrixXf> lltOfA(A);
@ -91,40 +173,69 @@ lltOfA.solveInPlace(b1);
// ...
\endcode
\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT.
\sa Cholesky_Module, MatrixBase::llt(), MatrixBase::ldlt(), LLT::solve(), LLT::solveInPlace(),
LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT.
\subsection TutorialAdvSolvers_LU LU decomposition (for most cases)
If the matrix \f$ A \f$ does not fit in any of the previous categories, or if you are unsure about the numerical
stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section \ref TutorialAdvLU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition
with full pivoting and rank update which has much better numerical stability.
The API of the LU solver is the same than the Cholesky one, except that there is no \em in \em place variant:
\code
#include <Eigen/LU>
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x;
A.lu().solve(b, &x);
\endcode
\subsection TutorialAdvSolvers_Triangular Triangular solver
Again, the LU decomposition can be stored to be reused or to perform other kernel operations:
\code
// ...
LU<MatrixXf> luOfA(A);
luOfA.solve(b, &x);
// ...
\endcode
If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the
diagonal are all not zero), then the problem can be solved directly using the TriangularView
class. This is much faster than using an %LU or Cholesky decomposition (in fact, the triangular
solver is used when you solve a system using the %LU or Cholesky decomposition). Here is an example:
See the section \ref TutorialAdvLU below.
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_triangular.cpp
</td><td>
output: \include Tutorial_solve_triangular.out
</td></tr></table>
\sa class LU, LU::solve(), LU_Module
The MatrixBase::triangularView() function constructs an object of the class TriangularView, and
TriangularView::solve() then solves the system. There is also an \e in \e place variant:
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_triangular_inplace.cpp
</td><td>
output: \include Tutorial_solve_triangular_inplace.out
</td></tr></table>
\sa MatrixBase::triangularView(), TriangularView::solve(), TriangularView::solveInPlace(),
TriangularView class.
\subsection TutorialAdvSolvers_SVD SVD solver (for singular matrices and special cases)
Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the
same than with Cholesky or LU:
\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices)
The solution of the system \f$ A \mathbf{x} = \mathbf{b} \f$ is given by \f$ \mathbf{x} = A^{-1}
\mathbf{b} \f$. This suggests the following approach for solving the system: compute the matrix
inverse and multiply that with the right-hand side. This is often not a good approach: using the %LU
decomposition with partial pivoting yields a more accurate algorithm that is usually just as fast or
even faster. However, using the matrix inverse can be faster if the matrix \f$ A \f$ is small
(&le;4) and fixed size, though numerical stability problems may still remain. Here is an example of
how you would write this in Eigen:
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_matrix_inverse.cpp
</td><td>
output: \include Tutorial_solve_matrix_inverse.out
</td></tr></table>
Note that the function inverse() is defined in the \ref LU_Module.
\sa MatrixBase::inverse().
\subsection TutorialAdvSolvers_Misc Other solvers (for singular matrices and special cases)
Finally, Eigen also offer solvers based on a singular value decomposition (%SVD) or the %LU
decomposition with full pivoting. These have the same API as the solvers based on the %LU
decomposition with partial pivoting (PartialLU).
The solver based on the %SVD uses the class SVD. It can handle singular matrices. Here is an example
of its use:
\code
#include <Eigen/SVD>
// ...
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x;
@ -133,8 +244,23 @@ SVD<MatrixXf> svdOfA(A);
svdOfA.solve(b, &x);
\endcode
\sa class SVD, SVD::solve(), SVD_Module
%LU decomposition with full pivoting has better numerical stability than %LU decomposition with
partial pivoting. It is defined in the class LU. The solver can also handle singular matrices.
\code
#include <Eigen/LU>
// ...
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x;
A.lu().solve(b, &x);
LU<MatrixXf> luOfA(A);
luOfA.solve(b, &x);
\endcode
See the section \ref TutorialAdvLU below.
\sa class SVD, SVD::solve(), SVD_Module, class LU, LU::solve(), LU_Module.

2
doc/C07_TutorialSparse.dox
View File

@ -1,6 +1,6 @@
namespace Eigen {
/** \page TutorialSparse Tutorial - Getting started with the sparse module
/** \page TutorialSparse Tutorial 4/4 - Getting started with the sparse module
\ingroup Tutorial
<div class="eimainmenu">\ref index "Overview"

18
doc/examples/Tutorial_PartialLU_solve.cpp Normal file
View File

@ -0,0 +1,18 @@
#include <Eigen/Core>
#include <Eigen/LU>
using namespace std;
using namespace Eigen;
int main(int, char *[])
{
Matrix3f A;
Vector3f b;
A << 1,2,3, 4,5,6, 7,8,10;
b << 3, 3, 4;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the vector b:" << endl << b << endl;
Vector3f x;
A.partialLu().solve(b, &x);
cout << "The solution is:" << endl << x << endl;
}

6
doc/snippets/Tutorial_solve_matrix_inverse.cpp Normal file
View File

@ -0,0 +1,6 @@
Matrix3f A;
Vector3f b;
A << 1,2,3, 4,5,6, 7,8,10;
b << 3, 3, 4;
Vector3f x = A.inverse() * b;
cout << "The solution is:" << endl << x << endl;

10
doc/snippets/Tutorial_solve_multiple_rhs.cpp Normal file
View File

@ -0,0 +1,10 @@
Matrix3f A(3,3);
A << 1,2,3, 4,5,6, 7,8,10;
Matrix<float,3,2> B;
B << 3,1, 3,1, 4,1;
Matrix<float,3,2> X;
A.partialLu().solve(B, &X);
cout << "The solution with right-hand side (3,3,4) is:" << endl;
cout << X.col(0) << endl;
cout << "The solution with right-hand side (1,1,1) is:" << endl;
cout << X.col(1) << endl;

13
doc/snippets/Tutorial_solve_reuse_decomposition.cpp Normal file
View File

@ -0,0 +1,13 @@
Matrix3f A(3,3);
A << 1,2,3, 4,5,6, 7,8,10;
PartialLU<Matrix3f> luOfA(A); // compute LU decomposition of A
Vector3f b;
b << 3,3,4;
Vector3f x;
luOfA.solve(b, &x);
cout << "The solution with right-hand side (3,3,4) is:" << endl;
cout << x << endl;
b << 1,1,1;
luOfA.solve(b, &x);
cout << "The solution with right-hand side (1,1,1) is:" << endl;
cout << x << endl;

9
doc/snippets/Tutorial_solve_singular.cpp Normal file
View File

@ -0,0 +1,9 @@
Matrix3f A;
Vector3f b;
A << 1,2,3, 4,5,6, 7,8,9;
b << 3, 3, 4;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the vector b:" << endl << b << endl;
Vector3f x;
A.partialLu().solve(b, &x);
cout << "The solution is:" << endl << x << endl;

8
doc/snippets/Tutorial_solve_triangular.cpp Normal file
View File

@ -0,0 +1,8 @@
Matrix3f A;
Vector3f b;
A << 1,2,3, 0,5,6, 0,0,10;
b << 3, 3, 4;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the vector b:" << endl << b << endl;
Vector3f x = A.triangularView<UpperTriangular>().solve(b);
cout << "The solution is:" << endl << x << endl;

6
doc/snippets/Tutorial_solve_triangular_inplace.cpp Normal file
View File

@ -0,0 +1,6 @@
Matrix3f A;
Vector3f b;
A << 1,2,3, 0,5,6, 0,0,10;
b << 3, 3, 4;
A.triangularView<UpperTriangular>().solveInPlace(b);
cout << "The solution is:" << endl << b << endl;

8
test/CMakeLists.txt
View File

@ -78,7 +78,6 @@ else(QT4_FOUND)
ei_add_property(EIGEN_MISSING_BACKENDS "Qt4 support, ")
endif(QT4_FOUND)
if(TEST_LIB)
add_definitions("-DEIGEN_EXTERN_INSTANTIATIONS=1")
endif(TEST_LIB)
@ -122,6 +121,8 @@ ei_add_test(lu ${EI_OFLAG})
ei_add_test(determinant)
ei_add_test(inverse ${EI_OFLAG})
ei_add_test(qr)
ei_add_test(qr_colpivoting)
ei_add_test(qr_fullpivoting)
ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}")
ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}")
ei_add_test(svd)
@ -137,11 +138,6 @@ ei_add_test(regression)
ei_add_test(stdvector)
ei_add_test(resize)
if(QT4_FOUND)
if(QT_QTCORE_LIBRARY_DEBUG)
set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_DEBUG})
else(QT_QTCORE_LIBRARY_DEBUG)
set(QT_QTCORE_LIBRARY ${QT_QTCORE_LIBRARY_RELEASE})
endif(QT_QTCORE_LIBRARY_DEBUG)
ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}")
endif(QT4_FOUND)
ei_add_test(sparse_vector)

4
test/lu.cpp
View File

@ -44,7 +44,7 @@ template<typename MatrixType> void lu_non_invertible()
VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
VERIFY(!lu.isInjective());
VERIFY(!lu.isInvertible());
VERIFY(lu.isSurjective() == (lu.rank() == rows));
VERIFY(!lu.isSurjective());
VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
VERIFY(m1image.lu().rank() == rank);
MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
@ -53,7 +53,7 @@ template<typename MatrixType> void lu_non_invertible()
m2 = MatrixType::Random(cols,cols2);
m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
lu.solve(m3, &m2);
VERIFY(lu.solve(m3, &m2));
VERIFY_IS_APPROX(m3, m1*m2);
m3 = MatrixType::Random(rows,cols2);
VERIFY(!lu.solve(m3, &m2));

9
test/main.h
View File

@ -227,12 +227,11 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m,
return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>());
}
template<typename Derived>
void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::MatrixBase<Derived>& m)
template<typename MatrixType>
void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
{
typedef Derived MatrixType;
typedef typename ei_traits<MatrixType>::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
typedef Matrix<Scalar, Dynamic, 1> VectorType;
MatrixType a = MatrixType::Random(rows,rows);
MatrixType d = MatrixType::Identity(rows,cols);
@ -244,7 +243,7 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri
HouseholderQR<MatrixType> qra(a);
HouseholderQR<MatrixType> qrb(b);
m = (qra.matrixQ() * d * qrb.matrixQ()).lazy();
m = qra.matrixQ() * d * qrb.matrixQ();
}
} // end namespace Eigen

116
test/qr_colpivoting.cpp Normal file
View File

@ -0,0 +1,116 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/QR>
template<typename MatrixType> void qr()
{
/* this test covers the following files: QR.h */
int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
int rank = ei_random<int>(1, std::min(rows, cols)-1);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
createRandomMatrixOfRank(rank,rows,cols,m1);
ColPivotingHouseholderQR<MatrixType> qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
MatrixType r = qr.matrixQR();
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
MatrixType b = qr.matrixQ() * r;
MatrixType c = MatrixType::Zero(rows,cols);
for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c);
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
qr.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2);
}
template<typename MatrixType> void qr_invertible()
{
/* this test covers the following files: RRQR.h */
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
int size = ei_random<int>(10,50);
MatrixType m1(size, size), m2(size, size), m3(size, size);
m1 = MatrixType::Random(size,size);
if (ei_is_same_type<RealScalar,float>::ret)
{
// let's build a matrix more stable to inverse
MatrixType a = MatrixType::Random(size,size*2);
m1 += a * a.adjoint();
}
ColPivotingHouseholderQR<MatrixType> qr(m1);
m3 = MatrixType::Random(size,size);
qr.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2);
}
template<typename MatrixType> void qr_verify_assert()
{
MatrixType tmp;
ColPivotingHouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixR())
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
}
void test_qr_colpivoting()
{
for(int i = 0; i < 1; i++) {
// FIXME : very weird bug here
// CALL_SUBTEST( qr(Matrix2f()) );
CALL_SUBTEST( qr<MatrixXf>() );
CALL_SUBTEST( qr<MatrixXd>() );
CALL_SUBTEST( qr<MatrixXcd>() );
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( qr_invertible<MatrixXf>() );
CALL_SUBTEST( qr_invertible<MatrixXd>() );
CALL_SUBTEST( qr_invertible<MatrixXcf>() );
CALL_SUBTEST( qr_invertible<MatrixXcd>() );
}
CALL_SUBTEST(qr_verify_assert<Matrix3f>());
CALL_SUBTEST(qr_verify_assert<Matrix3d>());
CALL_SUBTEST(qr_verify_assert<MatrixXf>());
CALL_SUBTEST(qr_verify_assert<MatrixXd>());
CALL_SUBTEST(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST(qr_verify_assert<MatrixXcd>());
}

137
test/qr_fullpivoting.cpp Normal file
View File

@ -0,0 +1,137 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
#include <Eigen/QR>
template<typename MatrixType> void qr()
{
/* this test covers the following files: QR.h */
int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
int rank = ei_random<int>(1, std::min(rows, cols)-1);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
createRandomMatrixOfRank(rank,rows,cols,m1);
FullPivotingHouseholderQR<MatrixType> qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
VERIFY(!qr.isInvertible());
VERIFY(!qr.isSurjective());
MatrixType r = qr.matrixQR();
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
MatrixType b = qr.matrixQ() * r;
MatrixType c = MatrixType::Zero(rows,cols);
for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c);
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
VERIFY(qr.solve(m3, &m2));
VERIFY_IS_APPROX(m3, m1*m2);
m3 = MatrixType::Random(rows,cols2);
VERIFY(!qr.solve(m3, &m2));
}
template<typename MatrixType> void qr_invertible()
{
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename MatrixType::Scalar Scalar;
int size = ei_random<int>(10,50);
MatrixType m1(size, size), m2(size, size), m3(size, size);
m1 = MatrixType::Random(size,size);
if (ei_is_same_type<RealScalar,float>::ret)
{
// let's build a matrix more stable to inverse
MatrixType a = MatrixType::Random(size,size*2);
m1 += a * a.adjoint();
}
FullPivotingHouseholderQR<MatrixType> qr(m1);
VERIFY(qr.isInjective());
VERIFY(qr.isInvertible());
VERIFY(qr.isSurjective());
m3 = MatrixType::Random(size,size);
VERIFY(qr.solve(m3, &m2));
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
RealScalar absdet = ei_abs(m1.diagonal().prod());
m3 = qr.matrixQ(); // get a unitary
m1 = m3 * m1 * m3;
qr.compute(m1);
VERIFY_IS_APPROX(absdet, qr.absDeterminant());
}
template<typename MatrixType> void qr_verify_assert()
{
MatrixType tmp;
FullPivotingHouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixR())
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
}
void test_qr_fullpivoting()
{
for(int i = 0; i < 1; i++) {
// FIXME : very weird bug here
// CALL_SUBTEST( qr(Matrix2f()) );
CALL_SUBTEST( qr<MatrixXf>() );
CALL_SUBTEST( qr<MatrixXd>() );
CALL_SUBTEST( qr<MatrixXcd>() );
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( qr_invertible<MatrixXf>() );
CALL_SUBTEST( qr_invertible<MatrixXd>() );
CALL_SUBTEST( qr_invertible<MatrixXcf>() );
CALL_SUBTEST( qr_invertible<MatrixXcd>() );
}
CALL_SUBTEST(qr_verify_assert<Matrix3f>());
CALL_SUBTEST(qr_verify_assert<Matrix3d>());
CALL_SUBTEST(qr_verify_assert<MatrixXf>());
CALL_SUBTEST(qr_verify_assert<MatrixXd>());
CALL_SUBTEST(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST(qr_verify_assert<MatrixXcd>());
}

2
test/stdvector.cpp
View File

@ -159,5 +159,5 @@ void test_stdvector()
// some Quaternion
CALL_SUBTEST(check_stdvector_quaternion(Quaternionf()));
CALL_SUBTEST(check_stdvector_quaternion(Quaternionf()));
CALL_SUBTEST(check_stdvector_quaternion(Quaterniond()));
}