merge with main repository

This commit is contained in:
Thomas Capricelli 2009-11-09 19:02:52 +01:00
commit 42b92c2022
123 changed files with 3010 additions and 2213 deletions

View File

@ -3,12 +3,12 @@ project(Eigen)
cmake_minimum_required(VERSION 2.6.2)
# automatically parse the version number
file(READ "${CMAKE_SOURCE_DIR}/Eigen/src/Core/util/Macros.h" _eigen2_version_header LIMIT 5000 OFFSET 1000)
string(REGEX MATCH "define *EIGEN_WORLD_VERSION ([0-9]*)" _eigen2_world_version_match "${_eigen2_version_header}")
file(READ "${CMAKE_SOURCE_DIR}/Eigen/src/Core/util/Macros.h" _eigen2_version_header)
string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen2_world_version_match "${_eigen2_version_header}")
set(EIGEN2_WORLD_VERSION "${CMAKE_MATCH_1}")
string(REGEX MATCH "define *EIGEN_MAJOR_VERSION ([0-9]*)" _eigen2_major_version_match "${_eigen2_version_header}")
string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen2_major_version_match "${_eigen2_version_header}")
set(EIGEN2_MAJOR_VERSION "${CMAKE_MATCH_1}")
string(REGEX MATCH "define *EIGEN_MINOR_VERSION ([0-9]*)" _eigen2_minor_version_match "${_eigen2_version_header}")
string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen2_minor_version_match "${_eigen2_version_header}")
set(EIGEN2_MINOR_VERSION "${CMAKE_MATCH_1}")
set(EIGEN_VERSION_NUMBER ${EIGEN2_WORLD_VERSION}.${EIGEN2_MAJOR_VERSION}.${EIGEN2_MINOR_VERSION})
@ -142,21 +142,23 @@ if(EIGEN_BUILD_BTL)
endif(EIGEN_BUILD_BTL)
ei_testing_print_summary()
if(NOT MSVC_IDE)
message("")
message("Configured Eigen ${EIGEN_VERSION_NUMBER}")
message("You can now do the following:")
message("--------------+----------------------------------------------------------------")
message("Command | Description")
message("--------------+----------------------------------------------------------------")
message("make install | Install to ${CMAKE_INSTALL_PREFIX}")
message(" | * To change that: cmake . -DCMAKE_INSTALL_PREFIX=yourpath")
message("make btest | Build the unit tests")
message(" | * That takes lots of memory! Easy on the -j option")
message("make test | Build and run the unit tests (using CTest)")
message("make test_qr | Build a specific test, here test_qr. To run it: test/test_qr")
message("make debug_qr | Build a test with full debug info. To run it: test/debug_qr")
message("make blas | Build BLAS library (not the same thing as Eigen)")
message("make doc | Generate the API documentation, requires Doxygen & LaTeX")
message("--------------+----------------------------------------------------------------")
endif(NOT MSVC_IDE)
string(TOLOWER "${CMAKE_GENERATOR}" cmake_generator_tolower)
if(cmake_generator_tolower MATCHES "makefile")
message("You can now do the following:")
message("--------------+----------------------------------------------------------------")
message("Command | Description")
message("--------------+----------------------------------------------------------------")
message("make install | Install to ${CMAKE_INSTALL_PREFIX}")
message(" | * To change that: cmake . -DCMAKE_INSTALL_PREFIX=yourpath")
message("make btest | Build the unit tests")
message("make test | Build and run the unit tests (using CTest)")
message("make test_qr | Build a specific test, here test_qr. To run it: test/test_qr")
message("make debug_qr | Build a test with full debug info. To run it: test/debug_qr")
message("make blas | Build BLAS library (not the same thing as Eigen)")
message("make doc | Generate the API documentation, requires Doxygen & LaTeX")
message("--------------+----------------------------------------------------------------")
endif()

View File

@ -1,25 +1,5 @@
set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi Eigenvalues)
if(EIGEN_BUILD_LIB)
set(Eigen_SRCS
src/Core/CoreInstantiations.cpp
src/Cholesky/CholeskyInstantiations.cpp
src/QR/QrInstantiations.cpp
)
add_library(Eigen2 SHARED ${Eigen_SRCS})
install(TARGETS Eigen2
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
endif(EIGEN_BUILD_LIB)
if(CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g1 -O2")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -g1 -O2")
endif(CMAKE_COMPILER_IS_GNUCXX)
install(FILES
${Eigen_HEADERS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen COMPONENT Devel

View File

@ -30,6 +30,7 @@ namespace Eigen {
* \endcode
*/
#include "src/misc/Solve.h"
#include "src/Array/CwiseOperators.h"
#include "src/Array/Functors.h"
#include "src/Cholesky/LLT.h"

View File

@ -18,8 +18,11 @@ namespace Eigen {
* \endcode
*/
#include "src/LU/LU.h"
#include "src/LU/PartialLU.h"
#include "src/misc/Solve.h"
#include "src/misc/Kernel.h"
#include "src/misc/Image.h"
#include "src/LU/FullPivLU.h"
#include "src/LU/PartialPivLU.h"
#include "src/LU/Determinant.h"
#include "src/LU/Inverse.h"

View File

@ -33,9 +33,10 @@ namespace Eigen {
* \endcode
*/
#include "src/misc/Solve.h"
#include "src/QR/HouseholderQR.h"
#include "src/QR/FullPivotingHouseholderQR.h"
#include "src/QR/ColPivotingHouseholderQR.h"
#include "src/QR/FullPivHouseholderQR.h"
#include "src/QR/ColPivHouseholderQR.h"
// declare all classes for a given matrix type
#define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \

View File

@ -22,6 +22,7 @@ namespace Eigen {
* \endcode
*/
#include "src/misc/Solve.h"
#include "src/SVD/SVD.h"
#include "src/SVD/JacobiSVD.h"

View File

@ -10,3 +10,4 @@ ADD_SUBDIRECTORY(Sparse)
ADD_SUBDIRECTORY(Jacobi)
ADD_SUBDIRECTORY(Householder)
ADD_SUBDIRECTORY(Eigenvalues)
ADD_SUBDIRECTORY(misc)

View File

@ -43,8 +43,8 @@
* zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
* on D also stabilizes the computation.
*
* Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky decomposition to determine
* whether a system of equations has a solution.
* Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
* decomposition to determine whether a system of equations has a solution.
*
* \sa MatrixBase::ldlt(), class LLT
*/
@ -52,10 +52,10 @@
* Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
* the strict lower part does not have to store correct values.
*/
template<typename MatrixType> class LDLT
template<typename _MatrixType> class LDLT
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
@ -88,7 +88,7 @@ template<typename MatrixType> class LDLT
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class LU.
* see the examples given in the documentation of class FullPivLU.
*/
inline const IntColVectorType& permutationP() const
{
@ -117,14 +117,40 @@ template<typename MatrixType> class LDLT
return m_sign == -1;
}
template<typename RhsDerived, typename ResultType>
bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
/** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* \note_about_checking_solutions
*
* \sa solveInPlace(), MatrixBase::ldlt()
*/
template<typename Rhs>
inline const ei_solve_retval<LDLT, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
ei_assert(m_matrix.rows()==b.rows()
&& "LDLT::solve(): invalid number of rows of the right hand side matrix b");
return ei_solve_retval<LDLT, Rhs>(*this, b.derived());
}
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const;
LDLT& compute(const MatrixType& matrix);
/** \returns the LDLT decomposition matrix
*
* TODO: document the storage layout
*/
inline const MatrixType& matrixLDLT() const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix;
}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
protected:
/** \internal
* Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
@ -134,7 +160,7 @@ template<typename MatrixType> class LDLT
*/
MatrixType m_matrix;
IntColVectorType m_p;
IntColVectorType m_transpositions;
IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
int m_sign;
bool m_isInitialized;
};
@ -238,27 +264,18 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
return *this;
}
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
* The result is stored in \a result
*
* \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
*
* In other words, it computes \f$ b = A^{-1} b \f$ with
* \f$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left.
*
* \sa LDLT::solveInPlace(), MatrixBase::ldlt()
*/
template<typename MatrixType>
template<typename RhsDerived, typename ResultType>
bool LDLT<MatrixType>
::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<LDLT<_MatrixType>, Rhs>
: ei_solve_retval_base<LDLT<_MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
*result = b;
return solveInPlace(*result);
}
EIGEN_MAKE_SOLVE_HELPERS(LDLT<_MatrixType>,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
dst = rhs();
dec().solveInPlace(dst);
}
};
/** This is the \em in-place version of solve().
*

View File

@ -53,9 +53,10 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
* Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
* the strict lower part does not have to store correct values.
*/
template<typename MatrixType, int _UpLo> class LLT
template<typename _MatrixType, int _UpLo> class LLT
{
private:
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
@ -68,8 +69,6 @@ template<typename MatrixType, int _UpLo> class LLT
typedef LLT_Traits<MatrixType,UpLo> Traits;
public:
/**
* \brief Default Constructor.
*
@ -99,14 +98,44 @@ template<typename MatrixType, int _UpLo> class LLT
return Traits::getL(m_matrix);
}
template<typename RhsDerived, typename ResultType>
bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* Since this LLT class assumes anyway that the matrix A is invertible, the solution
* theoretically exists and is unique regardless of b.
*
* Example: \include LLT_solve.cpp
* Output: \verbinclude LLT_solve.out
*
* \sa solveInPlace(), MatrixBase::llt()
*/
template<typename Rhs>
inline const ei_solve_retval<LLT, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
ei_assert(m_matrix.rows()==b.rows()
&& "LLT::solve(): invalid number of rows of the right hand side matrix b");
return ei_solve_retval<LLT, Rhs>(*this, b.derived());
}
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const;
LLT& compute(const MatrixType& matrix);
/** \returns the LLT decomposition matrix
*
* TODO: document the storage layout
*/
inline const MatrixType& matrixLLT() const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
return m_matrix;
}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
protected:
/** \internal
* Used to compute and store L
@ -229,28 +258,19 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
return *this;
}
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
* The result is stored in \a result
*
* \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
*
* In other words, it computes \f$ b = A^{-1} b \f$ with
* \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
*
* Example: \include LLT_solve.cpp
* Output: \verbinclude LLT_solve.out
*
* \sa LLT::solveInPlace(), MatrixBase::llt()
*/
template<typename MatrixType, int _UpLo>
template<typename RhsDerived, typename ResultType>
bool LLT<MatrixType,_UpLo>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
template<typename _MatrixType, int UpLo, typename Rhs>
struct ei_solve_retval<LLT<_MatrixType, UpLo>, Rhs>
: ei_solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
{
ei_assert(m_isInitialized && "LLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
return solveInPlace((*result) = b);
}
typedef LLT<_MatrixType,UpLo> LLTType;
EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
dst = rhs();
dec().solveInPlace(dst);
}
};
/** This is the \em in-place version of solve().
*

View File

@ -206,7 +206,7 @@ template<typename Derived> class MatrixBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal the plain matrix type corresponding to this expression. Note that is not necessarily
* exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const
* reference to a matrix, not a matrix! It guaranteed however, that the return type of eval() is either
* reference to a matrix, not a matrix! It is however guaranteed that the return type of eval() is either
* PlainMatrixType or const PlainMatrixType&.
*/
typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
@ -713,13 +713,23 @@ template<typename Derived> class MatrixBase
/////////// LU module ///////////
const LU<PlainMatrixType> lu() const;
const PartialLU<PlainMatrixType> partialLu() const;
const PlainMatrixType inverse() const;
const FullPivLU<PlainMatrixType> fullPivLu() const;
const PartialPivLU<PlainMatrixType> partialPivLu() const;
const PartialPivLU<PlainMatrixType> lu() const;
const ei_inverse_impl<Derived> inverse() const;
template<typename ResultType>
void computeInverse(ResultType *result) const;
void computeInverseAndDetWithCheck(
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible,
const RealScalar& absDeterminantThreshold = precision<Scalar>()
) const;
template<typename ResultType>
bool computeInverseWithCheck(ResultType *result ) const;
void computeInverseWithCheck(
ResultType& inverse,
bool& invertible,
const RealScalar& absDeterminantThreshold = precision<Scalar>()
) const;
Scalar determinant() const;
/////////// Cholesky module ///////////
@ -730,8 +740,8 @@ template<typename Derived> class MatrixBase
/////////// QR module ///////////
const HouseholderQR<PlainMatrixType> householderQr() const;
const ColPivotingHouseholderQR<PlainMatrixType> colPivotingHouseholderQr() const;
const FullPivotingHouseholderQR<PlainMatrixType> fullPivotingHouseholderQr() const;
const ColPivHouseholderQR<PlainMatrixType> colPivHouseholderQr() const;
const FullPivHouseholderQR<PlainMatrixType> fullPivHouseholderQr() const;
EigenvaluesReturnType eigenvalues() const;
RealScalar operatorNorm() const;

View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 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
@ -33,10 +34,23 @@ struct ei_traits<ReturnByValue<Derived> >
: public ei_traits<typename ei_traits<Derived>::ReturnMatrixType>
{
enum {
Flags = ei_traits<typename ei_traits<Derived>::ReturnMatrixType>::Flags | EvalBeforeNestingBit
// FIXME had to remove the DirectAccessBit for usage like
// matrix.inverse().block(...)
// because the Block ctor with direct access
// wants to call coeffRef() to get an address, and that fails (infinite recursion) as ReturnByValue
// doesnt implement coeffRef().
// The fact that I had to do that shows that when doing xpr.block() with a non-direct-access xpr,
// even if xpr has the EvalBeforeNestingBit, the block() doesn't use direct access on the evaluated
// xpr.
Flags = (ei_traits<typename ei_traits<Derived>::ReturnMatrixType>::Flags
| EvalBeforeNestingBit) & ~DirectAccessBit
};
};
/* The ReturnByValue object doesn't even have a coeff() method.
* So the only way that nesting it in an expression can work, is by evaluating it into a plain matrix.
* So ei_nested always gives the plain return matrix type.
*/
template<typename Derived,int n,typename PlainMatrixType>
struct ei_nested<ReturnByValue<Derived>, n, PlainMatrixType>
{
@ -46,9 +60,9 @@ struct ei_nested<ReturnByValue<Derived>, n, PlainMatrixType>
template<typename Derived>
class ReturnByValue : public MatrixBase<ReturnByValue<Derived> >
{
typedef typename ei_traits<Derived>::ReturnMatrixType ReturnMatrixType;
public:
EIGEN_GENERIC_PUBLIC_INTERFACE(ReturnByValue)
typedef typename ei_traits<Derived>::ReturnMatrixType ReturnMatrixType;
template<typename Dest>
inline void evalTo(Dest& dst) const
{ static_cast<const Derived* const>(this)->evalTo(dst); }

View File

@ -220,13 +220,13 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const Packet2
template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float* to, const Packet4f& from) { ei_pstoreu((double*)to, _mm_castps_pd(from)); }
template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int* to, const Packet4i& from) { ei_pstoreu((double*)to, _mm_castsi128_pd(from)); }
#if (_MSC_VER <= 1500) && defined(_WIN64)
#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64)
// The temporary variable fixes an internal compilation error.
// Direct of the struct members fixed bug #62.
template<> EIGEN_STRONG_INLINE float ei_pfirst<Packet4f>(const Packet4f& a) { return a.m128_f32[0]; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<Packet2d>(const Packet2d& a) { return a.m128d_f64[0]; }
template<> EIGEN_STRONG_INLINE int ei_pfirst<Packet4i>(const Packet4i& a) { int x = _mm_cvtsi128_si32(a); return x; }
#elif (_MSC_VER <= 1500)
#elif defined(_MSC_VER) && (_MSC_VER <= 1500)
// The temporary variable fixes an internal compilation error.
template<> EIGEN_STRONG_INLINE float ei_pfirst<Packet4f>(const Packet4f& a) { float x = _mm_cvtss_f32(a); return x; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<Packet2d>(const Packet2d& a) { double x = _mm_cvtsd_f64(a); return x; }

View File

@ -258,6 +258,11 @@ namespace {
EIGEN_UNUSED NoChange_t NoChange;
}
struct Default_t {};
namespace {
EIGEN_UNUSED Default_t Default;
}
enum {
IsDense = 0,
IsSparse = SparseBit,

View File

@ -66,6 +66,13 @@ template<typename ExpressionType> class WithFormat;
template<typename MatrixType> struct CommaInitializer;
template<typename Derived> class ReturnByValue;
template<typename DecompositionType, typename Rhs> struct ei_solve_retval_base;
template<typename DecompositionType, typename Rhs> struct ei_solve_retval;
template<typename DecompositionType> struct ei_kernel_retval_base;
template<typename DecompositionType> struct ei_kernel_retval;
template<typename DecompositionType> struct ei_image_retval_base;
template<typename DecompositionType> struct ei_image_retval;
template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int Supers=Dynamic, int Subs=Dynamic, int Options=0> class BandMatrix;
@ -114,11 +121,12 @@ template<typename ExpressionType, int Direction> class VectorwiseOp;
template<typename MatrixType,int RowFactor,int ColFactor> class Replicate;
template<typename MatrixType, int Direction = BothDirections> class Reverse;
template<typename MatrixType> class LU;
template<typename MatrixType> class PartialLU;
template<typename MatrixType> class FullPivLU;
template<typename MatrixType> class PartialPivLU;
template<typename MatrixType> struct ei_inverse_impl;
template<typename MatrixType> class HouseholderQR;
template<typename MatrixType> class ColPivotingHouseholderQR;
template<typename MatrixType> class FullPivotingHouseholderQR;
template<typename MatrixType> class ColPivHouseholderQR;
template<typename MatrixType> class FullPivHouseholderQR;
template<typename MatrixType> class SVD;
template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
template<typename MatrixType, int UpLo = LowerTriangular> class LLT;

View File

@ -30,7 +30,7 @@
#define EIGEN_WORLD_VERSION 2
#define EIGEN_MAJOR_VERSION 90
#define EIGEN_MINOR_VERSION 0
#define EIGEN_MINOR_VERSION 1
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \

View File

@ -26,6 +26,155 @@
#ifndef EIGEN_QUATERNION_H
#define EIGEN_QUATERNION_H
/***************************************************************************
* Definition of QuaternionBase<Derived>
* The implementation is at the end of the file
***************************************************************************/
template<typename Other,
int OtherRows=Other::RowsAtCompileTime,
int OtherCols=Other::ColsAtCompileTime>
struct ei_quaternionbase_assign_impl;
template<class Derived>
class QuaternionBase : public RotationBase<Derived, 3>
{
typedef RotationBase<Derived, 3> Base;
public:
using Base::operator*;
using Base::derived;
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename ei_traits<Derived>::Coefficients Coefficients;
// typedef typename Matrix<Scalar,4,1> Coefficients;
/** the type of a 3D vector */
typedef Matrix<Scalar,3,1> Vector3;
/** the equivalent rotation matrix type */
typedef Matrix<Scalar,3,3> Matrix3;
/** the equivalent angle-axis type */
typedef AngleAxis<Scalar> AngleAxisType;
/** \returns the \c x coefficient */
inline Scalar x() const { return this->derived().coeffs().coeff(0); }
/** \returns the \c y coefficient */
inline Scalar y() const { return this->derived().coeffs().coeff(1); }
/** \returns the \c z coefficient */
inline Scalar z() const { return this->derived().coeffs().coeff(2); }
/** \returns the \c w coefficient */
inline Scalar w() const { return this->derived().coeffs().coeff(3); }
/** \returns a reference to the \c x coefficient */
inline Scalar& x() { return this->derived().coeffs().coeffRef(0); }
/** \returns a reference to the \c y coefficient */
inline Scalar& y() { return this->derived().coeffs().coeffRef(1); }
/** \returns a reference to the \c z coefficient */
inline Scalar& z() { return this->derived().coeffs().coeffRef(2); }
/** \returns a reference to the \c w coefficient */
inline Scalar& w() { return this->derived().coeffs().coeffRef(3); }
/** \returns a read-only vector expression of the imaginary part (x,y,z) */
inline const VectorBlock<Coefficients,3> vec() const { return coeffs().template start<3>(); }
/** \returns a vector expression of the imaginary part (x,y,z) */
inline VectorBlock<Coefficients,3> vec() { return coeffs().template start<3>(); }
/** \returns a read-only vector expression of the coefficients (x,y,z,w) */
inline const typename ei_traits<Derived>::Coefficients& coeffs() const { return derived().coeffs(); }
/** \returns a vector expression of the coefficients (x,y,z,w) */
inline typename ei_traits<Derived>::Coefficients& coeffs() { return derived().coeffs(); }
template<class OtherDerived> Derived& operator=(const QuaternionBase<OtherDerived>& other);
// disabled this copy operator as it is giving very strange compilation errors when compiling
// test_stdvector with GCC 4.4.2. This looks like a GCC bug though, so feel free to re-enable it if it's
// useful; however notice that we already have the templated operator= above and e.g. in MatrixBase
// we didn't have to add, in addition to templated operator=, such a non-templated copy operator.
// Derived& operator=(const QuaternionBase& other)
// { return operator=<Derived>(other); }
Derived& operator=(const AngleAxisType& aa);
template<class OtherDerived> Derived& operator=(const MatrixBase<OtherDerived>& m);
/** \returns a quaternion representing an identity rotation
* \sa MatrixBase::Identity()
*/
inline static Quaternion<Scalar> Identity() { return Quaternion<Scalar>(1, 0, 0, 0); }
/** \sa QuaternionBase::Identity(), MatrixBase::setIdentity()
*/
inline QuaternionBase& setIdentity() { coeffs() << 0, 0, 0, 1; return *this; }
/** \returns the squared norm of the quaternion's coefficients
* \sa QuaternionBase::norm(), MatrixBase::squaredNorm()
*/
inline Scalar squaredNorm() const { return coeffs().squaredNorm(); }
/** \returns the norm of the quaternion's coefficients
* \sa QuaternionBase::squaredNorm(), MatrixBase::norm()
*/
inline Scalar norm() const { return coeffs().norm(); }
/** Normalizes the quaternion \c *this
* \sa normalized(), MatrixBase::normalize() */
inline void normalize() { coeffs().normalize(); }
/** \returns a normalized copy of \c *this
* \sa normalize(), MatrixBase::normalized() */
inline Quaternion<Scalar> normalized() const { return Quaternion<Scalar>(coeffs().normalized()); }
/** \returns the dot product of \c *this and \a other
* Geometrically speaking, the dot product of two unit quaternions
* corresponds to the cosine of half the angle between the two rotations.
* \sa angularDistance()
*/
template<class OtherDerived> inline Scalar dot(const QuaternionBase<OtherDerived>& other) const { return coeffs().dot(other.coeffs()); }
template<class OtherDerived> inline Scalar angularDistance(const QuaternionBase<OtherDerived>& other) const;
Matrix3 toRotationMatrix() const;
template<typename Derived1, typename Derived2>
Derived& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b);
template<class OtherDerived> inline Quaternion<Scalar> operator* (const QuaternionBase<OtherDerived>& q) const;
template<class OtherDerived> inline Derived& operator*= (const QuaternionBase<OtherDerived>& q);
Quaternion<Scalar> inverse() const;
Quaternion<Scalar> conjugate() const;
template<class OtherDerived> Quaternion<Scalar> slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const;
/** \returns \c true if \c *this is approximately equal to \a other, within the precision
* determined by \a prec.
*
* \sa MatrixBase::isApprox() */
template<class OtherDerived>
bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = precision<Scalar>()) const
{ return coeffs().isApprox(other.coeffs(), prec); }
Vector3 _transformVector(Vector3 v) const;
/** \returns \c *this with scalar type casted to \a NewScalarType
*
* Note that if \a NewScalarType is equal to the current scalar type of \c *this
* then this function smartly returns a const reference to \c *this.
*/
template<typename NewScalarType>
inline typename ei_cast_return_type<Derived,Quaternion<NewScalarType> >::type cast() const
{
return typename ei_cast_return_type<Derived,Quaternion<NewScalarType> >::type(
coeffs().template cast<NewScalarType>());
}
};
/***************************************************************************
* Definition/implementation of Quaternion<Scalar>
***************************************************************************/
/** \geometry_module \ingroup Geometry_Module
*
* \class Quaternion
@ -48,152 +197,13 @@
* \sa class AngleAxis, class Transform
*/
template<typename Other,
int OtherRows=Other::RowsAtCompileTime,
int OtherCols=Other::ColsAtCompileTime>
struct ei_quaternionbase_assign_impl;
template<typename Scalar> class Quaternion; // [XXX] => remove when Quaternion becomes Quaternion
template<typename Derived>
struct ei_traits<QuaternionBase<Derived> >
{
typedef typename ei_traits<Derived>::Scalar Scalar;
enum {
PacketAccess = ei_traits<Derived>::PacketAccess
};
};
template<class Derived>
class QuaternionBase : public RotationBase<Derived, 3>
{
typedef RotationBase<Derived, 3> Base;
public:
using Base::operator*;
typedef typename ei_traits<QuaternionBase<Derived> >::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
// typedef typename Matrix<Scalar,4,1> Coefficients;
/** the type of a 3D vector */
typedef Matrix<Scalar,3,1> Vector3;
/** the equivalent rotation matrix type */
typedef Matrix<Scalar,3,3> Matrix3;
/** the equivalent angle-axis type */
typedef AngleAxis<Scalar> AngleAxisType;
/** \returns the \c x coefficient */
inline Scalar x() const { return this->derived().coeffs().coeff(0); }
/** \returns the \c y coefficient */
inline Scalar y() const { return this->derived().coeffs().coeff(1); }
/** \returns the \c z coefficient */
inline Scalar z() const { return this->derived().coeffs().coeff(2); }
/** \returns the \c w coefficient */
inline Scalar w() const { return this->derived().coeffs().coeff(3); }
/** \returns a reference to the \c x coefficient */
inline Scalar& x() { return this->derived().coeffs().coeffRef(0); }
/** \returns a reference to the \c y coefficient */
inline Scalar& y() { return this->derived().coeffs().coeffRef(1); }
/** \returns a reference to the \c z coefficient */
inline Scalar& z() { return this->derived().coeffs().coeffRef(2); }
/** \returns a reference to the \c w coefficient */
inline Scalar& w() { return this->derived().coeffs().coeffRef(3); }
/** \returns a read-only vector expression of the imaginary part (x,y,z) */
inline const VectorBlock<typename ei_traits<Derived>::Coefficients,3> vec() const { return this->derived().coeffs().template start<3>(); }
/** \returns a vector expression of the imaginary part (x,y,z) */
inline VectorBlock<typename ei_traits<Derived>::Coefficients,3> vec() { return this->derived().coeffs().template start<3>(); }
/** \returns a read-only vector expression of the coefficients (x,y,z,w) */
inline const typename ei_traits<Derived>::Coefficients& coeffs() const { return this->derived().coeffs(); }
/** \returns a vector expression of the coefficients (x,y,z,w) */
inline typename ei_traits<Derived>::Coefficients& coeffs() { return this->derived().coeffs(); }
template<class OtherDerived> QuaternionBase& operator=(const QuaternionBase<OtherDerived>& other);
QuaternionBase& operator=(const AngleAxisType& aa);
template<class OtherDerived>
QuaternionBase& operator=(const MatrixBase<OtherDerived>& m);
/** \returns a quaternion representing an identity rotation
* \sa MatrixBase::Identity()
*/
inline static Quaternion<Scalar> Identity() { return Quaternion<Scalar>(1, 0, 0, 0); }
/** \sa Quaternion2::Identity(), MatrixBase::setIdentity()
*/
inline QuaternionBase& setIdentity() { coeffs() << 0, 0, 0, 1; return *this; }
/** \returns the squared norm of the quaternion's coefficients
* \sa Quaternion2::norm(), MatrixBase::squaredNorm()
*/
inline Scalar squaredNorm() const { return coeffs().squaredNorm(); }
/** \returns the norm of the quaternion's coefficients
* \sa Quaternion2::squaredNorm(), MatrixBase::norm()
*/
inline Scalar norm() const { return coeffs().norm(); }
/** Normalizes the quaternion \c *this
* \sa normalized(), MatrixBase::normalize() */
inline void normalize() { coeffs().normalize(); }
/** \returns a normalized version of \c *this
* \sa normalize(), MatrixBase::normalized() */
inline Quaternion<Scalar> normalized() const { return Quaternion<Scalar>(coeffs().normalized()); }
/** \returns the dot product of \c *this and \a other
* Geometrically speaking, the dot product of two unit quaternions
* corresponds to the cosine of half the angle between the two rotations.
* \sa angularDistance()
*/
template<class OtherDerived> inline Scalar dot(const QuaternionBase<OtherDerived>& other) const { return coeffs().dot(other.coeffs()); }
template<class OtherDerived> inline Scalar angularDistance(const QuaternionBase<OtherDerived>& other) const;
Matrix3 toRotationMatrix(void) const;
template<typename Derived1, typename Derived2>
QuaternionBase& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b);
template<class OtherDerived> inline Quaternion<Scalar> operator* (const QuaternionBase<OtherDerived>& q) const;
template<class OtherDerived> inline QuaternionBase& operator*= (const QuaternionBase<OtherDerived>& q);
Quaternion<Scalar> inverse(void) const;
Quaternion<Scalar> conjugate(void) const;
template<class OtherDerived> Quaternion<Scalar> slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const;
/** \returns \c true if \c *this is approximately equal to \a other, within the precision
* determined by \a prec.
*
* \sa MatrixBase::isApprox() */
bool isApprox(const QuaternionBase& other, RealScalar prec = precision<Scalar>()) const
{ return coeffs().isApprox(other.coeffs(), prec); }
Vector3 _transformVector(Vector3 v) const;
/** \returns \c *this with scalar type casted to \a NewScalarType
*
* Note that if \a NewScalarType is equal to the current scalar type of \c *this
* then this function smartly returns a const reference to \c *this.
*/
template<typename NewScalarType>
inline typename ei_cast_return_type<Derived,Quaternion<NewScalarType> >::type cast() const
{
return typename ei_cast_return_type<Derived,Quaternion<NewScalarType> >::type(
coeffs().template cast<NewScalarType>());
}
};
template<typename _Scalar>
struct ei_traits<Quaternion<_Scalar> >
{
typedef _Scalar Scalar;
typedef Matrix<_Scalar,4,1> Coefficients;
enum{
PacketAccess = Aligned
PacketAccess = Aligned
};
};
@ -239,7 +249,7 @@ public:
explicit inline Quaternion(const MatrixBase<Derived>& other) { *this = other; }
/** Copy constructor with scalar type conversion */
template<class Derived>
template<typename Derived>
inline explicit Quaternion(const QuaternionBase<Derived>& other)
{ m_coeffs = other.coeffs().template cast<Scalar>(); }
@ -250,16 +260,29 @@ protected:
Coefficients m_coeffs;
};
/* ########### Map<Quaternion> */
/** \ingroup Geometry_Module
* single precision quaternion type */
typedef Quaternion<float> Quaternionf;
/** \ingroup Geometry_Module
* double precision quaternion type */
typedef Quaternion<double> Quaterniond;
/***************************************************************************
* Specialization of Map<Quaternion<Scalar>>
***************************************************************************/
/** \class Map<Quaternion>
* \nonstableyet
*
* \brief Expression of a quaternion
* \brief Expression of a quaternion from a memory buffer
*
* \param Scalar the type of the vector of diagonal coefficients
* \param _Scalar the type of the Quaternion coefficients
* \param PacketAccess see class Map
*
* \sa class Quaternion, class QuaternionBase
* This is a specialization of class Map for Quaternion. This class allows to view
* a 4 scalar memory buffer as an Eigen's Quaternion object.
*
* \sa class Map, class Quaternion, class QuaternionBase
*/
template<typename _Scalar, int _PacketAccess>
struct ei_traits<Map<Quaternion<_Scalar>, _PacketAccess> >:
@ -273,15 +296,23 @@ ei_traits<Quaternion<_Scalar> >
};
template<typename _Scalar, int PacketAccess>
class Map<Quaternion<_Scalar>, PacketAccess > : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >, ei_no_assignment_operator {
class Map<Quaternion<_Scalar>, PacketAccess >
: public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >,
ei_no_assignment_operator
{
public:
typedef _Scalar Scalar;
typedef typename ei_traits<Map>::Coefficients Coefficients;
typedef typename ei_traits<Map<Quaternion<Scalar>, PacketAccess> >::Coefficients Coefficients;
/** Constructs a Mapped Quaternion object from the pointer \a coeffs
*
* The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
* \code *coeffs == {x, y, z, w} \endcode
*
* If the template paramter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */
inline Map(const Scalar* coeffs) : m_coeffs(coeffs) {}
inline Map<Quaternion<Scalar>, PacketAccess >(const Scalar* coeffs) : m_coeffs(coeffs) {}
inline Coefficients& coeffs() { return m_coeffs;}
inline const Coefficients& coeffs() const { return m_coeffs;}
@ -289,15 +320,20 @@ class Map<Quaternion<_Scalar>, PacketAccess > : public QuaternionBase<Map<Quater
Coefficients m_coeffs;
};
typedef Map<Quaternion<double> > QuaternionMapd;
typedef Map<Quaternion<float> > QuaternionMapf;
typedef Map<Quaternion<double>, Aligned> QuaternionMapAlignedd;
typedef Map<Quaternion<float>, Aligned> QuaternionMapAlignedf;
typedef Map<Quaternion<double> > QuaternionMapd;
typedef Map<Quaternion<float> > QuaternionMapf;
typedef Map<Quaternion<double>, Aligned> QuaternionMapAlignedd;
typedef Map<Quaternion<float>, Aligned> QuaternionMapAlignedf;
/***************************************************************************
* Implementation of QuaternionBase methods
***************************************************************************/
// Generic Quaternion * Quaternion product
template<int Arch, class Derived, class OtherDerived, typename Scalar, int PacketAccess> struct ei_quat_product
// This product can be specialized for a given architecture via the Arch template argument.
template<int Arch, class Derived1, class Derived2, typename Scalar, int PacketAccess> struct ei_quat_product
{
inline static Quaternion<Scalar> run(const QuaternionBase<Derived>& a, const QuaternionBase<OtherDerived>& b){
inline static Quaternion<Scalar> run(const QuaternionBase<Derived1>& a, const QuaternionBase<Derived2>& b){
return Quaternion<Scalar>
(
a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(),
@ -311,21 +347,22 @@ template<int Arch, class Derived, class OtherDerived, typename Scalar, int Packe
/** \returns the concatenation of two rotations as a quaternion-quaternion product */
template <class Derived>
template <class OtherDerived>
inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) const
inline Quaternion<typename ei_traits<Derived>::Scalar>
QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) const
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename OtherDerived::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
return ei_quat_product<EiArch, Derived, OtherDerived,
typename ei_traits<Derived>::Scalar,
ei_traits<Derived>::PacketAccess && ei_traits<OtherDerived>::PacketAccess>::run(*this, other);
return ei_quat_product<EiArch, Derived, OtherDerived,
typename ei_traits<Derived>::Scalar,
ei_traits<Derived>::PacketAccess && ei_traits<OtherDerived>::PacketAccess>::run(*this, other);
}
/** \sa operator*(Quaternion) */
template <class Derived>
template <class OtherDerived>
inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator*= (const QuaternionBase<OtherDerived>& other)
inline Derived& QuaternionBase<Derived>::operator*= (const QuaternionBase<OtherDerived>& other)
{
return (*this = *this * other);
return (derived() = derived() * other.derived());
}
/** Rotation of a vector by a quaternion.
@ -350,21 +387,21 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const
template<class Derived>
template<class OtherDerived>
inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const QuaternionBase<OtherDerived>& other)
inline Derived& QuaternionBase<Derived>::operator=(const QuaternionBase<OtherDerived>& other)
{
coeffs() = other.coeffs();
return *this;
return derived();
}
/** Set \c *this from an angle-axis \a aa and returns a reference to \c *this
*/
template<class Derived>
inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const AngleAxisType& aa)
inline Derived& QuaternionBase<Derived>::operator=(const AngleAxisType& aa)
{
Scalar ha = Scalar(0.5)*aa.angle(); // Scalar(0.5) to suppress precision loss warnings
this->w() = ei_cos(ha);
this->vec() = ei_sin(ha) * aa.axis();
return *this;
return derived();
}
/** Set \c *this from the expression \a xpr:
@ -375,12 +412,12 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const AngleAx
template<class Derived>
template<class MatrixDerived>
inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const MatrixBase<MatrixDerived>& xpr)
inline Derived& QuaternionBase<Derived>::operator=(const MatrixBase<MatrixDerived>& xpr)
{
EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename MatrixDerived::Scalar>::ret),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
ei_quaternionbase_assign_impl<MatrixDerived>::run(*this, xpr.derived());
return *this;
return derived();
}
/** Convert the quaternion to a 3x3 rotation matrix. The quaternion is required to
@ -434,7 +471,7 @@ QuaternionBase<Derived>::toRotationMatrix(void) const
*/
template<class Derived>
template<typename Derived1, typename Derived2>
inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b)
inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b)
{
Vector3 v0 = a.normalized();
Vector3 v1 = b.normalized();
@ -458,7 +495,7 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const
Scalar w2 = (Scalar(1)+c)*Scalar(0.5);
this->w() = ei_sqrt(w2);
this->vec() = axis * ei_sqrt(Scalar(1) - w2);
return *this;
return derived();
}
Vector3 axis = v0.cross(v1);
Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2));
@ -466,17 +503,17 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const
this->vec() = axis * invs;
this->w() = s * Scalar(0.5);
return *this;
return derived();
}
/** \returns the multiplicative inverse of \c *this
* Note that in most cases, i.e., if you simply want the opposite rotation,
* and/or the quaternion is normalized, then it is enough to use the conjugate.
*
* \sa Quaternion2::conjugate()
* \sa QuaternionBase::conjugate()
*/
template <class Derived>
inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::inverse() const
inline Quaternion<typename ei_traits<Derived>::Scalar> QuaternionBase<Derived>::inverse() const
{
// FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ??
Scalar n2 = this->squaredNorm();
@ -485,7 +522,7 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni
else
{
// return an invalid result to flag the error
return Quaternion<Scalar>(ei_traits<Derived>::Coefficients::Zero());
return Quaternion<Scalar>(Coefficients::Zero());
}
}
@ -496,7 +533,8 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni
* \sa Quaternion2::inverse()
*/
template <class Derived>
inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::conjugate() const
inline Quaternion<typename ei_traits<Derived>::Scalar>
QuaternionBase<Derived>::conjugate() const
{
return Quaternion<Scalar>(this->w(),-this->x(),-this->y(),-this->z());
}
@ -506,11 +544,12 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni
*/
template <class Derived>
template <class OtherDerived>
inline typename ei_traits<QuaternionBase<Derived> >::Scalar QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const
inline typename ei_traits<Derived>::Scalar
QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const
{
double d = ei_abs(this->dot(other));
if (d>=1.0)
return 0;
return Scalar(0);
return Scalar(2) * std::acos(d);
}
@ -519,13 +558,14 @@ inline typename ei_traits<QuaternionBase<Derived> >::Scalar QuaternionBase<Deriv
*/
template <class Derived>
template <class OtherDerived>
Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const
Quaternion<typename ei_traits<Derived>::Scalar>
QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const
{
static const Scalar one = Scalar(1) - precision<Scalar>();
Scalar d = this->dot(other);
Scalar absD = ei_abs(d);
if (absD>=one)
return Quaternion<Scalar>(*this);
return Quaternion<Scalar>(derived());
// theta is the angle between the 2 quaternions
Scalar theta = std::acos(absD);
@ -549,7 +589,7 @@ struct ei_quaternionbase_assign_impl<Other,3,3>
// This algorithm comes from "Quaternion Calculus and Fast Animation",
// Ken Shoemake, 1987 SIGGRAPH course notes
Scalar t = mat.trace();
if (t > 0)
if (t > Scalar(0))
{
t = ei_sqrt(t + Scalar(1.0));
q.w() = Scalar(0.5)*t;

View File

@ -885,6 +885,24 @@ Transform<Scalar,Dim,Mode>::fromPositionOrientationScale(const MatrixBase<Positi
return *this;
}
// selector needed to avoid taking the inverse of a 3x4 matrix
template<typename TransformType, int Mode=TransformType::Mode>
struct ei_projective_transform_inverse
{
static inline void run(const TransformType&, TransformType&)
{}
};
template<typename TransformType>
struct ei_projective_transform_inverse<TransformType, Projective>
{
static inline void run(const TransformType& m, TransformType& res)
{
res.matrix() = m.matrix().inverse();
}
};
/** \nonstableyet
*
* \returns the inverse transformation according to some given knowledge
@ -911,7 +929,7 @@ Transform<Scalar,Dim,Mode>::inverse(TransformTraits hint) const
Transform res;
if (hint == Projective)
{
res.matrix() = m_matrix.inverse();
ei_projective_transform_inverse<Transform>::run(*this, res);
}
else
{

View File

@ -53,7 +53,7 @@ template<typename Derived,
{
static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
{
return m.partialLu().determinant();
return m.partialPivLu().determinant();
}
};

684
Eigen/src/LU/FullPivLU.h Normal file
View File

@ -0,0 +1,684 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-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_LU_H
#define EIGEN_LU_H
/** \ingroup LU_Module
*
* \class FullPivLU
*
* \brief LU decomposition of a matrix with complete pivoting, and related features
*
* \param MatrixType the type of the matrix of which we are computing the LU decomposition
*
* This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A
* is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q
* are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal
* coefficients) of U are sorted in such a way that any zeros are at the end.
*
* 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. 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().
*
* As an exemple, here is how the original matrix can be retrieved:
* \include class_FullPivLU.cpp
* Output: \verbinclude class_FullPivLU.out
*
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
*/
template<typename _MatrixType> class FullPivLU
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LU::compute(const MatrixType&).
*/
FullPivLU();
/** Constructor.
*
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*/
FullPivLU(const MatrixType& matrix);
/** Computes the LU decomposition of the given matrix.
*
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*
* \returns a reference to *this
*/
FullPivLU& compute(const MatrixType& matrix);
/** \returns the LU decomposition matrix: the upper-triangular part is U, the
* unit-lower-triangular part is L (at least for square matrices; in the non-square
* case, special care is needed, see the documentation of class FullPivLU).
*
* \sa matrixL(), matrixU()
*/
inline const MatrixType& matrixLU() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_lu;
}
/** \returns the number of nonzero pivots in the LU decomposition.
* Here nonzero is meant in the exact sense, not in a fuzzy sense.
* So that notion isn't really intrinsically interesting, but it is
* still useful when implementing algorithms.
*
* \sa rank()
*/
inline int nonzeroPivots() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_nonzero_pivots;
}
/** \returns the absolute value of the biggest pivot, i.e. the biggest
* diagonal coefficient of U.
*/
RealScalar maxPivot() const { return m_maxpivot; }
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class FullPivLU.
*
* \sa permutationQ()
*/
inline const IntColVectorType& permutationP() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_p;
}
/** \returns a vector of integers, whose size is the number of columns of the matrix being
* decomposed, representing the Q permutation i.e. the permutation of the columns.
* For its precise meaning, see the examples given in the documentation of class FullPivLU.
*
* \sa permutationP()
*/
inline const IntRowVectorType& permutationQ() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_q;
}
/** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
* will form a basis of the kernel.
*
* \note If the kernel has dimension zero, then the returned matrix is a column-vector filled with zeros.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*
* Example: \include FullPivLU_kernel.cpp
* Output: \verbinclude FullPivLU_kernel.out
*
* \sa image()
*/
inline const ei_kernel_retval<FullPivLU> kernel() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return ei_kernel_retval<FullPivLU>(*this);
}
/** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
* will form a basis of the kernel.
*
* \param originalMatrix the original matrix, of which *this is the LU decomposition.
* The reason why it is needed to pass it here, is that this allows
* a large optimization, as otherwise this method would need to reconstruct it
* from the LU decomposition.
*
* \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*
* Example: \include FullPivLU_image.cpp
* Output: \verbinclude FullPivLU_image.out
*
* \sa kernel()
*/
inline const ei_image_retval<FullPivLU>
image(const MatrixType& originalMatrix) const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return ei_image_retval<FullPivLU>(*this, originalMatrix);
}
/** \return a solution x to the equation Ax=b, where A is the matrix of which
* *this is the LU decomposition.
*
* \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
* the only requirement in order for the equation to make sense is that
* b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
*
* \returns a solution.
*
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
* \note_about_using_kernel_to_study_multiple_solutions
*
* Example: \include FullPivLU_solve.cpp
* Output: \verbinclude FullPivLU_solve.out
*
* \sa TriangularView::solve(), kernel(), inverse()
*/
template<typename Rhs>
inline const ei_solve_retval<FullPivLU, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return ei_solve_retval<FullPivLU, Rhs>(*this, b.derived());
}
/** \returns the determinant of the matrix of which
* *this is the LU decomposition. It has only linear complexity
* (that is, O(n) where n is the dimension of the square matrix)
* as the LU decomposition has already been computed.
*
* \note This is only for square matrices.
*
* \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
* optimized paths.
*
* \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 ei_traits<MatrixType>::Scalar determinant() const;
/** Allows to prescribe a threshold to be used by certain methods, such as rank(),
* who need to determine when pivots are to be considered nonzero. This is not used for the
* LU decomposition itself.
*
* When it needs to get the threshold value, Eigen calls threshold(). By default, this calls
* defaultThreshold(). Once you have called the present method setThreshold(const RealScalar&),
* your value is used instead.
*
* \param threshold The new value to use as the threshold.
*
* A pivot will be considered nonzero if its absolute value is strictly greater than
* \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
* where maxpivot is the biggest pivot.
*
* If you want to come back to the default behavior, call setThreshold(Default_t)
*/
FullPivLU& setThreshold(const RealScalar& threshold)
{
m_usePrescribedThreshold = true;
m_prescribedThreshold = threshold;
}
/** Allows to come back to the default behavior, letting Eigen use its default formula for
* determining the threshold.
*
* You should pass the special object Eigen::Default as parameter here.
* \code lu.setThreshold(Eigen::Default); \endcode
*
* See the documentation of setThreshold(const RealScalar&).
*/
FullPivLU& setThreshold(Default_t)
{
m_usePrescribedThreshold = false;
}
/** Returns the threshold that will be used by certain methods such as rank().
*
* See the documentation of setThreshold(const RealScalar&).
*/
RealScalar threshold() const
{
ei_assert(m_isInitialized || m_usePrescribedThreshold);
return m_usePrescribedThreshold ? m_prescribedThreshold
// this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt.
: epsilon<Scalar>() * m_lu.diagonalSize();
}
/** \returns the rank of the matrix of which *this is the LU decomposition.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline int rank() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
RealScalar premultiplied_threshold = ei_abs(m_maxpivot) * threshold();
int result = 0;
for(int i = 0; i < m_nonzero_pivots; ++i)
result += (ei_abs(m_lu.coeff(i,i)) > premultiplied_threshold);
return result;
}
/** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline int dimensionOfKernel() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_lu.cols() - rank();
}
/** \returns true if the matrix of which *this is the LU decomposition represents an injective
* linear map, i.e. has trivial kernel; false otherwise.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline bool isInjective() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return rank() == m_lu.cols();
}
/** \returns true if the matrix of which *this is the LU decomposition represents a surjective
* linear map; false otherwise.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline bool isSurjective() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return rank() == m_lu.rows();
}
/** \returns true if the matrix of which *this is the LU decomposition is invertible.
*
* \note This method has to determine which pivots should be considered nonzero.
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline bool isInvertible() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return isInjective() && (m_lu.rows() == m_lu.cols());
}
/** \returns the inverse of the matrix of which *this is the LU 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 MatrixBase::inverse()
*/
inline const ei_solve_retval<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
return ei_solve_retval<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> >
(*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
}
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
protected:
MatrixType m_lu;
IntColVectorType m_p;
IntRowVectorType m_q;
int m_det_pq, m_nonzero_pivots;
RealScalar m_maxpivot, m_prescribedThreshold;
bool m_isInitialized, m_usePrescribedThreshold;
};
template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU()
: m_isInitialized(false), m_usePrescribedThreshold(false)
{
}
template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
: m_isInitialized(false), m_usePrescribedThreshold(false)
{
compute(matrix);
}
template<typename MatrixType>
FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
{
m_isInitialized = true;
m_lu = matrix;
m_p.resize(matrix.rows());
m_q.resize(matrix.cols());
const int size = matrix.diagonalSize();
const int rows = matrix.rows();
const int cols = matrix.cols();
// will store the transpositions, before we accumulate them at the end.
// can't accumulate on-the-fly because that will be done in reverse order for the rows.
IntColVectorType rows_transpositions(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols());
int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. rows_transpositions[i]!=i
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
for(int k = 0; k < size; ++k)
{
// First, we need to find the pivot.
// biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
int row_of_biggest_in_corner, col_of_biggest_in_corner;
RealScalar biggest_in_corner;
biggest_in_corner = m_lu.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; // correct the values! since they were computed in the corner,
col_of_biggest_in_corner += k; // need to add k to them.
// if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values
if(biggest_in_corner == RealScalar(0))
{
// before exiting, make sure to initialize the still uninitialized row_transpositions
// in a sane state without destroying what we already have.
m_nonzero_pivots = k;
for(int i = k; i < size; i++)
{
rows_transpositions.coeffRef(i) = i;
cols_transpositions.coeffRef(i) = i;
}
break;
}
if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
// Now that we've found the pivot, we need to apply the row/col swaps to
// bring it to the location (k,k).
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_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
++number_of_transpositions;
}
if(k != col_of_biggest_in_corner) {
m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
++number_of_transpositions;
}
// Now that the pivot is at the right location, we update the remaining
// bottom-right corner by Gaussian elimination.
if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= m_lu.coeff(k,k);
if(k<size-1)
m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).end(rows-k-1) * m_lu.row(k).end(cols-k-1);
}
// the main loop is over, we still have to accumulate the transpositions to find the
// permutations P and Q
for(int k = 0; k < matrix.rows(); ++k) m_p.coeffRef(k) = k;
for(int k = size-1; k >= 0; --k)
std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
for(int k = 0; k < matrix.cols(); ++k) m_q.coeffRef(k) = k;
for(int k = 0; k < size; ++k)
std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
return *this;
}
template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
{
ei_assert(m_isInitialized && "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();
}
/********* Implementation of kernel() **************************************************/
template<typename _MatrixType>
struct ei_kernel_retval<FullPivLU<_MatrixType> >
: ei_kernel_retval_base<FullPivLU<_MatrixType> >
{
EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
MatrixType::MaxColsAtCompileTime,
MatrixType::MaxRowsAtCompileTime)
};
template<typename Dest> void evalTo(Dest& dst) const
{
const int cols = dec().matrixLU().cols(), dimker = cols - rank();
if(dimker == 0)
{
// The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
// avoid crashing/asserting as that depends on floating point calculations. Let's
// just return a single column vector filled with zeros.
dst.setZero();
return;
}
/* Let us use the following lemma:
*
* Lemma: If the matrix A has the LU decomposition PAQ = LU,
* then Ker A = Q(Ker U).
*
* Proof: trivial: just keep in mind that P, Q, L are invertible.
*/
/* Thus, all we need to do is to compute Ker U, and then apply Q.
*
* U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
* Thus, the diagonal of U ends with exactly
* dimKer zero's. Let us use that to construct dimKer linearly
* independent vectors in Ker U.
*/
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
int p = 0;
for(int i = 0; i < dec().nonzeroPivots(); ++i)
if(ei_abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
pivots.coeffRef(p++) = i;
ei_internal_assert(p == rank());
// we construct a temporaty trapezoid matrix m, by taking the U matrix and
// permuting the rows and cols to bring the nonnegligible pivots to the top of
// the main diagonal. We need that to be able to apply our triangular solvers.
// FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
m(dec().matrixLU().block(0, 0, rank(), cols));
for(int i = 0; i < rank(); ++i)
{
if(i) m.row(i).start(i).setZero();
m.row(i).end(cols-i) = dec().matrixLU().row(pivots.coeff(i)).end(cols-i);
}
m.block(0, 0, rank(), rank());
m.block(0, 0, rank(), rank()).template triangularView<StrictlyLowerTriangular>().setZero();
for(int i = 0; i < rank(); ++i)
m.col(i).swap(m.col(pivots.coeff(i)));
// ok, we have our trapezoid matrix, we can apply the triangular solver.
// notice that the math behind this suggests that we should apply this to the
// negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
m.corner(TopLeft, rank(), rank())
.template triangularView<UpperTriangular>().solveInPlace(
m.corner(TopRight, rank(), dimker)
);
// now we must undo the column permutation that we had applied!
for(int i = rank()-1; i >= 0; --i)
m.col(i).swap(m.col(pivots.coeff(i)));
// see the negative sign in the next line, that's what we were talking about above.
for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().coeff(i)) = -m.row(i).end(dimker);
for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().coeff(i)).setZero();
for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().coeff(rank()+k), k) = Scalar(1);
}
};
/***** Implementation of image() *****************************************************/
template<typename _MatrixType>
struct ei_image_retval<FullPivLU<_MatrixType> >
: ei_image_retval_base<FullPivLU<_MatrixType> >
{
EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
MatrixType::MaxColsAtCompileTime,
MatrixType::MaxRowsAtCompileTime)
};
template<typename Dest> void evalTo(Dest& dst) const
{
if(rank() == 0)
{
// The Image is just {0}, so it doesn't have a basis properly speaking, but let's
// avoid crashing/asserting as that depends on floating point calculations. Let's
// just return a single column vector filled with zeros.
dst.setZero();
return;
}
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
int p = 0;
for(int i = 0; i < dec().nonzeroPivots(); ++i)
if(ei_abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
pivots.coeffRef(p++) = i;
ei_internal_assert(p == rank());
for(int i = 0; i < rank(); ++i)
dst.col(i) = originalMatrix().col(dec().permutationQ().coeff(pivots.coeff(i)));
}
};
/***** Implementation of solve() *****************************************************/
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
: ei_solve_retval_base<FullPivLU<_MatrixType>, Rhs>
{
EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
* So we proceed as follows:
* Step 1: compute c = P * rhs.
* Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
* Step 3: replace c by the solution x to Ux = c. May or may not exist.
* Step 4: result = Q * c;
*/
const int rows = dec().matrixLU().rows(), cols = dec().matrixLU().cols(),
nonzero_pivots = dec().nonzeroPivots();
ei_assert(rhs().rows() == rows);
const int smalldim = std::min(rows, cols);
if(nonzero_pivots == 0)
{
dst.setZero();
return;
}
typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols());
// Step 1
for(int i = 0; i < rows; ++i)
c.row(dec().permutationP().coeff(i)) = rhs().row(i);
// Step 2
dec().matrixLU()
.corner(Eigen::TopLeft,smalldim,smalldim)
.template triangularView<UnitLowerTriangular>()
.solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols)
{
c.corner(Eigen::BottomLeft, rows-cols, c.cols())
-= dec().matrixLU().corner(Eigen::BottomLeft, rows-cols, cols)
* c.corner(Eigen::TopLeft, cols, c.cols());
}
// Step 3
dec().matrixLU()
.corner(TopLeft, nonzero_pivots, nonzero_pivots)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols()));
// Step 4
for(int i = 0; i < nonzero_pivots; ++i)
dst.row(dec().permutationQ().coeff(i)) = c.row(i);
for(int i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
dst.row(dec().permutationQ().coeff(i)).setZero();
}
};
/******* MatrixBase methods *****************************************************************/
/** \lu_module
*
* \return the full-pivoting LU decomposition of \c *this.
*
* \sa class FullPivLU
*/
template<typename Derived>
inline const FullPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::fullPivLu() const
{
return FullPivLU<PlainMatrixType>(eval());
}
#endif // EIGEN_LU_H

View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-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
@ -25,78 +25,165 @@
#ifndef EIGEN_INVERSE_H
#define EIGEN_INVERSE_H
/********************************************************************
*** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
********************************************************************/
/**********************************
*** General case implementation ***
**********************************/
template<typename XprType, typename MatrixType>
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, ResultType& result)
{
result = matrix.partialPivLu().inverse();
}
};
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse_and_det_with_check { /* nothing! general case not supported. */ };
/****************************
*** Size 1 implementation ***
****************************/
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 1>
{
static inline void run(const MatrixType& matrix, ResultType& result)
{
typedef typename MatrixType::Scalar Scalar;
result.coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 1>
{
static inline void run(
const MatrixType& matrix,
const typename MatrixType::RealScalar& absDeterminantThreshold,
ResultType& result,
typename ResultType::Scalar& determinant,
bool& invertible
)
{
determinant = matrix.coeff(0,0);
invertible = ei_abs(determinant) > absDeterminantThreshold;
if(invertible) result.coeffRef(0,0) = typename ResultType::Scalar(1) / determinant;
}
};
/****************************
*** Size 2 implementation ***
****************************/
template<typename MatrixType, typename ResultType>
inline void ei_compute_inverse_size2_helper(
const XprType& matrix, const typename MatrixType::Scalar& invdet,
MatrixType* result)
const MatrixType& matrix, const typename ResultType::Scalar& invdet,
ResultType& result)
{
result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
}
template<typename MatrixType>
inline void ei_compute_inverse_size2(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar invdet = Scalar(1) / matrix.determinant();
ei_compute_inverse_size2_helper( matrix, invdet, result );
}
template<typename XprType, typename MatrixType>
bool ei_compute_inverse_size2_with_check(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det = matrix.determinant();
if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
ei_compute_inverse_size2_helper( matrix, invdet, result );
return true;
}
template<typename XprType, typename MatrixType>
void ei_compute_inverse_size3_helper(
const XprType& matrix,
const typename MatrixType::Scalar& invdet,
const typename MatrixType::Scalar& det_minor00,
const typename MatrixType::Scalar& det_minor10,
const typename MatrixType::Scalar& det_minor20,
MatrixType* result)
{
result->coeffRef(0, 0) = det_minor00 * invdet;
result->coeffRef(0, 1) = -det_minor10 * invdet;
result->coeffRef(0, 2) = det_minor20 * invdet;
result->coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
result->coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
result->coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
}
template<bool Check, typename XprType, typename MatrixType>
bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar det = ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );
if(Check) if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
ei_compute_inverse_size3_helper( matrix, invdet, det_minor00, det_minor10, det_minor20, result );
return true;
result.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result.coeffRef(1,1) = matrix.coeff(0,0) * invdet;
}
template<typename MatrixType, typename ResultType>
bool ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType* result)
struct ei_compute_inverse<MatrixType, ResultType, 2>
{
static inline void run(const MatrixType& matrix, ResultType& result)
{
typedef typename ResultType::Scalar Scalar;
const Scalar invdet = typename MatrixType::Scalar(1) / matrix.determinant();
ei_compute_inverse_size2_helper(matrix, invdet, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
{
static inline void run(
const MatrixType& matrix,
const typename MatrixType::RealScalar& absDeterminantThreshold,
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible
)
{
typedef typename ResultType::Scalar Scalar;
determinant = matrix.determinant();
invertible = ei_abs(determinant) > absDeterminantThreshold;
if(!invertible) return;
const Scalar invdet = Scalar(1) / determinant;
ei_compute_inverse_size2_helper(matrix, invdet, inverse);
}
};
/****************************
*** Size 3 implementation ***
****************************/
template<typename MatrixType, typename ResultType>
void ei_compute_inverse_size3_helper(
const MatrixType& matrix,
const typename ResultType::Scalar& invdet,
const Matrix<typename ResultType::Scalar,3,1>& cofactors_col0,
ResultType& result)
{
result.row(0) = cofactors_col0 * invdet;
result.coeffRef(1,0) = -matrix.minor(0,1).determinant() * invdet;
result.coeffRef(1,1) = matrix.minor(1,1).determinant() * invdet;
result.coeffRef(1,2) = -matrix.minor(2,1).determinant() * invdet;
result.coeffRef(2,0) = matrix.minor(0,2).determinant() * invdet;
result.coeffRef(2,1) = -matrix.minor(1,2).determinant() * invdet;
result.coeffRef(2,2) = matrix.minor(2,2).determinant() * invdet;
}
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 3>
{
static inline void run(const MatrixType& matrix, ResultType& result)
{
typedef typename ResultType::Scalar Scalar;
Matrix<Scalar,3,1> cofactors_col0;
cofactors_col0.coeffRef(0) = matrix.minor(0,0).determinant();
cofactors_col0.coeffRef(1) = -matrix.minor(1,0).determinant();
cofactors_col0.coeffRef(2) = matrix.minor(2,0).determinant();
const Scalar det = (cofactors_col0.cwise()*matrix.col(0)).sum();
const Scalar invdet = Scalar(1) / det;
ei_compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
{
static inline void run(
const MatrixType& matrix,
const typename MatrixType::RealScalar& absDeterminantThreshold,
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible
)
{
typedef typename ResultType::Scalar Scalar;
Matrix<Scalar,3,1> cofactors_col0;
cofactors_col0.coeffRef(0) = matrix.minor(0,0).determinant();
cofactors_col0.coeffRef(1) = -matrix.minor(1,0).determinant();
cofactors_col0.coeffRef(2) = matrix.minor(2,0).determinant();
determinant = (cofactors_col0.cwise()*matrix.col(0)).sum();
invertible = ei_abs(determinant) > absDeterminantThreshold;
if(!invertible) return;
const Scalar invdet = Scalar(1) / determinant;
ei_compute_inverse_size3_helper(matrix, invdet, cofactors_col0, inverse);
}
};
/****************************
*** Size 4 implementation ***
****************************/
template<typename MatrixType, typename ResultType>
void ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType& result)
{
/* Let's split M into four 2x2 blocks:
* (P Q)
@ -111,257 +198,224 @@ bool ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType* resul
* Q' = -(P_inverse*Q) * S'
* R' = -S' * (R*P_inverse)
*/
typedef Block<MatrixType,2,2> XprBlock22;
typedef Block<ResultType,2,2> XprBlock22;
typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
Block22 P_inverse;
if(ei_compute_inverse_size2_with_check(matrix.template block<2,2>(0,0), &P_inverse))
{
const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q;
const XprBlock22 R = matrix.template block<2,2>(2,0);
const Block22 R_times_P_inverse = R * P_inverse;
const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y;
ei_compute_inverse_size2(X, &Y);
result->template block<2,2>(2,2) = Y;
result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
const Block22 Z = P_inverse_times_Q * Y;
result->template block<2,2>(0,2) = - Z;
result->template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
return true;
}
else
{
return false;
}
ei_compute_inverse<XprBlock22, Block22>::run(matrix.template block<2,2>(0,0), P_inverse);
const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q;
const XprBlock22 R = matrix.template block<2,2>(2,0);
const Block22 R_times_P_inverse = R * P_inverse;
const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y;
ei_compute_inverse<Block22, Block22>::run(X, Y);
result.template block<2,2>(2,2) = Y;
result.template block<2,2>(2,0) = - Y * R_times_P_inverse;
const Block22 Z = P_inverse_times_Q * Y;
result.template block<2,2>(0,2) = - Z;
result.template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
}
template<typename XprType, typename MatrixType>
bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* result)
{
if(ei_compute_inverse_size4_helper(matrix, result))
{
// good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
return true;
}
else
{
// rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
// since this is a rare case, we don't need to optimize it. We just want to handle it with little
// additional code.
MatrixType m(matrix);
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
if(ei_compute_inverse_size4_helper(m, result))
{
// good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
// rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
// the corresponding columns.
result->col(0).swap(result->col(2));
result->col(1).swap(result->col(3));
return true;
}
else
{
// first, undo the swaps previously made
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
// swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs
int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1;
m.row(0).swap(m.row(swap0with));
// swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
m.row(1).swap(m.row(swap1with));
if( ei_compute_inverse_size4_helper(m, result) )
{
result->col(1).swap(result->col(swap1with));
result->col(0).swap(result->col(swap0with));
return true;
}
else
{
// non-invertible matrix
return false;
}
}
}
}
/***********************************************
*** Part 2 : selector and MatrixBase methods ***
***********************************************/
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, ResultType* result)
{
matrix.partialLu().computeInverse(result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 1>
{
static inline void run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 2>
{
static inline void run(const MatrixType& matrix, ResultType* result)
{
ei_compute_inverse_size2(matrix, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 3>
{
static inline void run(const MatrixType& matrix, ResultType* result)
{
ei_compute_inverse_size3<false, MatrixType, ResultType>(matrix, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse<MatrixType, ResultType, 4>
{
static inline void run(const MatrixType& matrix, ResultType* result)
static inline void run(const MatrixType& _matrix, ResultType& result)
{
ei_compute_inverse_size4_with_check(matrix, result);
typedef typename ResultType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
// we will do row permutations on the matrix. This copy should have negligible cost.
// if not, consider working in-place on the matrix (const-cast it, but then undo the permutations
// to nevertheless honor constness)
typename MatrixType::PlainMatrixType matrix(_matrix);
// let's extract from the 2 first colums a 2x2 block whose determinant is as big as possible.
int good_row0, good_row1, good_i;
Matrix<RealScalar,6,1> absdet;
// any 2x2 block with determinant above this threshold will be considered good enough
RealScalar d = (matrix.col(0).squaredNorm()+matrix.col(1).squaredNorm()) * RealScalar(1e-2);
#define ei_inv_size4_helper_macro(i,row0,row1) \
absdet[i] = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1) \
- matrix.coeff(row0,1)*matrix.coeff(row1,0)); \
if(absdet[i] > d) { good_row0=row0; good_row1=row1; goto good; }
ei_inv_size4_helper_macro(0,0,1)
ei_inv_size4_helper_macro(1,0,2)
ei_inv_size4_helper_macro(2,0,3)
ei_inv_size4_helper_macro(3,1,2)
ei_inv_size4_helper_macro(4,1,3)
ei_inv_size4_helper_macro(5,2,3)
// no 2x2 block has determinant bigger than the threshold. So just take the one that
// has the biggest determinant
absdet.maxCoeff(&good_i);
good_row0 = good_i <= 2 ? 0 : good_i <= 4 ? 1 : 2;
good_row1 = good_i <= 2 ? good_i+1 : good_i <= 4 ? good_i-1 : 3;
// now good_row0 and good_row1 are correctly set
good:
// do row permutations to move this 2x2 block to the top
matrix.row(0).swap(matrix.row(good_row0));
matrix.row(1).swap(matrix.row(good_row1));
// now applying our helper function is numerically stable
ei_compute_inverse_size4_helper(matrix, result);
// Since we did row permutations on the original matrix, we need to do column permutations
// in the reverse order on the inverse
result.col(1).swap(result.col(good_row1));
result.col(0).swap(result.col(good_row0));
}
};
/** \lu_module
*
* Computes the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
* computeInverseWithCheck().
*
* \param result Pointer to the matrix in which to store the result.
*
* Example: \include MatrixBase_computeInverse.cpp
* Output: \verbinclude MatrixBase_computeInverse.out
*
* \sa inverse(), computeInverseWithCheck()
*/
template<typename Derived>
template<typename ResultType>
inline void MatrixBase<Derived>::computeInverse(ResultType *result) const
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
{
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
ei_compute_inverse<PlainMatrixType, ResultType>::run(eval(), result);
}
static inline void run(
const MatrixType& matrix,
const typename MatrixType::RealScalar& absDeterminantThreshold,
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible
)
{
determinant = matrix.determinant();
invertible = ei_abs(determinant) > absDeterminantThreshold;
if(invertible) ei_compute_inverse<MatrixType, ResultType>::run(matrix, inverse);
}
};
/*************************
*** MatrixBase methods ***
*************************/
template<typename MatrixType>
struct ei_traits<ei_inverse_impl<MatrixType> >
{
typedef typename MatrixType::PlainMatrixType ReturnMatrixType;
};
template<typename MatrixType>
struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> >
{
// for 2x2, it's worth giving a chance to avoid evaluating.
// for larger sizes, evaluating has negligible cost and limits code size.
typedef typename ei_meta_if<
MatrixType::RowsAtCompileTime == 2,
typename ei_nested<MatrixType,2>::type,
typename ei_eval<MatrixType>::type
>::ret MatrixTypeNested;
typedef typename ei_cleantype<MatrixTypeNested>::type MatrixTypeNestedCleaned;
const MatrixTypeNested m_matrix;
ei_inverse_impl(const MatrixType& matrix)
: m_matrix(matrix)
{}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{
ei_compute_inverse<MatrixTypeNestedCleaned, Dest>::run(m_matrix, dst);
}
};
/** \lu_module
*
* \returns the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
* computeInverseWithCheck().
* For small fixed sizes up to 4x4, this method uses ad-hoc methods (cofactors up to 3x3, Euler's trick for 4x4).
* In the general case, this method uses class PartialPivLU.
*
* \note This method returns a matrix by value, which can be inefficient. To avoid that overhead,
* use computeInverse() instead.
* \note This matrix must be invertible, otherwise the result is undefined. If you need an
* invertibility check, do the following:
* \li for fixed sizes up to 4x4, use computeInverseAndDetWithCheck().
* \li for the general case, use class FullPivLU.
*
* Example: \include MatrixBase_inverse.cpp
* Output: \verbinclude MatrixBase_inverse.out
*
* \sa computeInverse(), computeInverseWithCheck()
* \sa computeInverseAndDetWithCheck()
*/
template<typename Derived>
inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const
inline const ei_inverse_impl<Derived> MatrixBase<Derived>::inverse() const
{
PlainMatrixType result(rows(), cols());
computeInverse(&result);
return result;
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
ei_assert(rows() == cols());
return ei_inverse_impl<Derived>(derived());
}
/********************************************
* Compute inverse with invertibility check *
*******************************************/
template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse_with_check
/** \lu_module
*
* Computation of matrix inverse and determinant, with invertibility check.
*
* This is only for fixed-size square matrices of size up to 4x4.
*
* \param inverse Reference to the matrix in which to store the inverse.
* \param determinant Reference to the variable in which to store the inverse.
* \param invertible Reference to the bool variable in which to store whether the matrix is invertible.
* \param absDeterminantThreshold Optional parameter controlling the invertibility check.
* The matrix will be declared invertible if the absolute value of its
* determinant is greater than this threshold.
*
* Example: \include MatrixBase_computeInverseAndDetWithCheck.cpp
* Output: \verbinclude MatrixBase_computeInverseAndDetWithCheck.out
*
* \sa inverse(), computeInverseWithCheck()
*/
template<typename Derived>
template<typename ResultType>
inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(
ResultType& inverse,
typename ResultType::Scalar& determinant,
bool& invertible,
const RealScalar& absDeterminantThreshold
) const
{
static inline bool run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
LU<MatrixType> lu( matrix );
if( !lu.isInvertible() ) return false;
lu.computeInverse(result);
return true;
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 1>
{
static inline bool run(const MatrixType& matrix, ResultType* result)
{
typedef typename MatrixType::Scalar Scalar;
if( matrix.coeff(0,0) == Scalar(0) ) return false;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
return true;
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 2>
{
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size2_with_check(matrix, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 3>
{
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size3<true, MatrixType, ResultType>(matrix, result);
}
};
template<typename MatrixType, typename ResultType>
struct ei_compute_inverse_with_check<MatrixType, ResultType, 4>
{
static inline bool run(const MatrixType& matrix, ResultType* result)
{
return ei_compute_inverse_size4_with_check(matrix, result);
}
};
// i'd love to put some static assertions there, but SFINAE means that they have no effect...
ei_assert(rows() == cols());
// for 2x2, it's worth giving a chance to avoid evaluating.
// for larger sizes, evaluating has negligible cost and limits code size.
typedef typename ei_meta_if<
RowsAtCompileTime == 2,
typename ei_cleantype<typename ei_nested<Derived, 2>::type>::type,
PlainMatrixType
>::ret MatrixType;
ei_compute_inverse_and_det_with_check<MatrixType, ResultType>::run
(derived(), absDeterminantThreshold, inverse, determinant, invertible);
}
/** \lu_module
*
* Computation of matrix inverse, with invertibility check.
*
* \returns true if the matrix is invertible, false otherwise.
* This is only for fixed-size square matrices of size up to 4x4.
*
* \param result Pointer to the matrix in which to store the result.
* \param inverse Reference to the matrix in which to store the inverse.
* \param invertible Reference to the bool variable in which to store whether the matrix is invertible.
* \param absDeterminantThreshold Optional parameter controlling the invertibility check.
* The matrix will be declared invertible if the absolute value of its
* determinant is greater than this threshold.
*
* \sa inverse(), computeInverse()
* Example: \include MatrixBase_computeInverseWithCheck.cpp
* Output: \verbinclude MatrixBase_computeInverseWithCheck.out
*
* \sa inverse(), computeInverseAndDetWithCheck()
*/
template<typename Derived>
template<typename ResultType>
inline bool MatrixBase<Derived>::computeInverseWithCheck(ResultType *result) const
inline void MatrixBase<Derived>::computeInverseWithCheck(
ResultType& inverse,
bool& invertible,
const RealScalar& absDeterminantThreshold
) const
{
RealScalar determinant;
// i'd love to put some static assertions there, but SFINAE means that they have no effect...
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
return ei_compute_inverse_with_check<PlainMatrixType, ResultType>::run(eval(), result);
computeInverseAndDetWithCheck(inverse,determinant,invertible,absDeterminantThreshold);
}
#endif // EIGEN_INVERSE_H

View File

@ -1,607 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 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_LU_H
#define EIGEN_LU_H
/** \ingroup LU_Module
*
* \class LU
*
* \brief LU decomposition of a matrix with complete pivoting, and related features
*
* \param MatrixType the type of the matrix of which we are computing the LU decomposition
*
* This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A
* is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q
* are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal
* coefficients) of U are sorted in such a way that any zeros are at the end, so that the rank
* of A is the index of the first zero on the diagonal of U (with indices starting at 0) if any.
*
* 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. 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().
*
* As an exemple, here is how the original matrix can be retrieved:
* \include class_LU.cpp
* Output: \verbinclude class_LU.out
*
* \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse()
*/
template<typename MatrixType> class LU
{
public:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
MatrixType::MaxColsAtCompileTime,
MatrixType::MaxRowsAtCompileTime)
};
typedef Matrix<typename MatrixType::Scalar,
MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" is the number of cols of the original matrix
// so that the product "matrix * kernel = zero" makes sense
Dynamic, // we don't know at compile-time the dimension of the kernel
MatrixType::Options,
MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter
MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, whose dimension is the number
// of columns of the original matrix
> KernelResultType;
typedef Matrix<typename MatrixType::Scalar,
MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose dimension is the number
// of rows of the original matrix
Dynamic, // we don't know at compile time the dimension of the image (the rank)
MatrixType::Options,
MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix,
MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns.
> ImageResultType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LU::compute(const MatrixType&).
*/
LU();
/** Constructor.
*
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*/
LU(const MatrixType& matrix);
/** Computes the LU decomposition of the given matrix.
*
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*
* \returns a reference to *this
*/
LU& compute(const MatrixType& matrix);
/** \returns the LU decomposition matrix: the upper-triangular part is U, the
* unit-lower-triangular part is L (at least for square matrices; in the non-square
* case, special care is needed, see the documentation of class LU).
*
* \sa matrixL(), matrixU()
*/
inline const MatrixType& matrixLU() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_lu;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class LU.
*
* \sa permutationQ()
*/
inline const IntColVectorType& permutationP() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_p;
}
/** \returns a vector of integers, whose size is the number of columns of the matrix being
* decomposed, representing the Q permutation i.e. the permutation of the columns.
* For its precise meaning, see the examples given in the documentation of class LU.
*
* \sa permutationP()
*/
inline const IntRowVectorType& permutationQ() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_q;
}
/** Computes a basis of the kernel of the matrix, also called the null-space of the matrix.
*
* \note This method is only allowed on non-invertible matrices, as determined by
* isInvertible(). Calling it on an invertible matrix will make an assertion fail.
*
* \param result a pointer to the matrix in which to store the kernel. The columns of this
* matrix will be set to form a basis of the kernel (it will be resized
* if necessary).
*
* Example: \include LU_computeKernel.cpp
* Output: \verbinclude LU_computeKernel.out
*
* \sa kernel(), computeImage(), image()
*/
template<typename KernelMatrixType>
void computeKernel(KernelMatrixType *result) const;
/** Computes a basis of the image of the matrix, also called the column-space or range of he matrix.
*
* \note Calling this method on the zero matrix will make an assertion fail.
*
* \param result a pointer to the matrix in which to store the image. The columns of this
* matrix will be set to form a basis of the image (it will be resized
* if necessary).
*
* Example: \include LU_computeImage.cpp
* Output: \verbinclude LU_computeImage.out
*
* \sa image(), computeKernel(), kernel()
*/
template<typename ImageMatrixType>
void computeImage(ImageMatrixType *result) const;
/** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
* will form a basis of the kernel.
*
* \note: this method is only allowed on non-invertible matrices, as determined by
* isInvertible(). Calling it on an invertible matrix will make an assertion fail.
*
* \note: this method returns a matrix by value, which induces some inefficiency.
* If you prefer to avoid this overhead, use computeKernel() instead.
*
* Example: \include LU_kernel.cpp
* Output: \verbinclude LU_kernel.out
*
* \sa computeKernel(), image()
*/
const KernelResultType kernel() const;
/** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
* will form a basis of the kernel.
*
* \note: Calling this method on the zero matrix will make an assertion fail.
*
* \note: this method returns a matrix by value, which induces some inefficiency.
* If you prefer to avoid this overhead, use computeImage() instead.
*
* Example: \include LU_image.cpp
* Output: \verbinclude LU_image.out
*
* \sa computeImage(), kernel()
*/
const ImageResultType image() const;
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
* *this is the LU decomposition, if any exists.
*
* \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
* the only requirement in order for the equation to make sense is that
* b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
* \param result a pointer to the vector or 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 true if any solution exists, false if no solution exists.
*
* \note If there exist more than one solution, this method will arbitrarily choose one.
* If you need a complete analysis of the space of solutions, take the one solution obtained
* by this method and add to it elements of the kernel, as determined by kernel().
*
* Example: \include LU_solve.cpp
* Output: \verbinclude LU_solve.out
*
* \sa TriangularView::solve(), kernel(), computeKernel(), inverse(), computeInverse()
*/
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
/** \returns the determinant of the matrix of which
* *this is the LU decomposition. It has only linear complexity
* (that is, O(n) where n is the dimension of the square matrix)
* as the LU decomposition has already been computed.
*
* \note This is only for square matrices.
*
* \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
* optimized paths.
*
* \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 ei_traits<MatrixType>::Scalar determinant() const;
/** \returns the rank of the matrix of which *this is the LU decomposition.
*
* \note This is computed at the time of the construction of the LU decomposition. This
* method does not perform any further computation.
*/
inline int rank() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_rank;
}
/** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition.
*
* \note Since the rank is computed at the time of the construction of the LU decomposition, this
* method almost does not perform any further computation.
*/
inline int dimensionOfKernel() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_lu.cols() - m_rank;
}
/** \returns true if the matrix of which *this is the LU 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 LU decomposition, this
* method almost does not perform any further computation.
*/
inline bool isInjective() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_rank == m_lu.cols();
}
/** \returns true if the matrix of which *this is the LU decomposition represents a surjective
* linear map; false otherwise.
*
* \note Since the rank is computed at the time of the construction of the LU decomposition, this
* method almost does not perform any further computation.
*/
inline bool isSurjective() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return m_rank == m_lu.rows();
}
/** \returns true if the matrix of which *this is the LU decomposition is invertible.
*
* \note Since the rank is computed at the time of the construction of the LU decomposition, this
* method almost does not perform any further computation.
*/
inline bool isInvertible() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
return isInjective() && isSurjective();
}
/** Computes the inverse of the matrix of which *this is the LU 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 MatrixBase::computeInverse(), inverse()
*/
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);
}
/** \returns the inverse of the matrix of which *this is the LU 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(), MatrixBase::inverse()
*/
inline MatrixType inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
}
protected:
const MatrixType* m_originalMatrix;
MatrixType m_lu;
IntColVectorType m_p;
IntRowVectorType m_q;
int m_det_pq;
int m_rank;
RealScalar m_precision;
};
template<typename MatrixType>
LU<MatrixType>::LU()
: m_originalMatrix(0),
m_lu(),
m_p(),
m_q(),
m_det_pq(0),
m_rank(-1),
m_precision(precision<RealScalar>())
{
}
template<typename MatrixType>
LU<MatrixType>::LU(const MatrixType& matrix)
: m_originalMatrix(0),
m_lu(),
m_p(),
m_q(),
m_det_pq(0),
m_rank(-1),
m_precision(precision<RealScalar>())
{
compute(matrix);
}
template<typename MatrixType>
LU<MatrixType>& LU<MatrixType>::compute(const MatrixType& matrix)
{
m_originalMatrix = &matrix;
m_lu = matrix;
m_p.resize(matrix.rows());
m_q.resize(matrix.cols());
const int size = matrix.diagonalSize();
const int rows = matrix.rows();
const int cols = matrix.cols();
// this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt.
m_precision = epsilon<Scalar>() * size;
IntColVectorType rows_transpositions(matrix.rows());
IntRowVectorType cols_transpositions(matrix.cols());
int number_of_transpositions = 0;
RealScalar biggest = RealScalar(0);
m_rank = size;
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_lu.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++)
{
rows_transpositions.coeffRef(i) = i;
cols_transpositions.coeffRef(i) = i;
}
break;
}
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_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
++number_of_transpositions;
}
if(k != col_of_biggest_in_corner) {
m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
++number_of_transpositions;
}
if(k<rows-1)
m_lu.col(k).end(rows-k-1) /= m_lu.coeff(k,k);
if(k<size-1)
m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).end(rows-k-1) * m_lu.row(k).end(cols-k-1);
}
for(int k = 0; k < matrix.rows(); ++k) m_p.coeffRef(k) = k;
for(int k = size-1; k >= 0; --k)
std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
for(int k = 0; k < matrix.cols(); ++k) m_q.coeffRef(k) = k;
for(int k = 0; k < size; ++k)
std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k)));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
return *this;
}
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();
}
template<typename MatrixType>
template<typename KernelMatrixType>
void LU<MatrixType>::computeKernel(KernelMatrixType *result) const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
const int dimker = dimensionOfKernel(), cols = m_lu.cols();
result->resize(cols, dimker);
/* Let us use the following lemma:
*
* Lemma: If the matrix A has the LU decomposition PAQ = LU,
* then Ker A = Q(Ker U).
*
* Proof: trivial: just keep in mind that P, Q, L are invertible.
*/
/* Thus, all we need to do is to compute Ker U, and then apply Q.
*
* U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
* Thus, the diagonal of U ends with exactly
* m_dimKer zero's. Let us use that to construct m_dimKer linearly
* independent vectors in Ker U.
*/
Matrix<Scalar, Dynamic, Dynamic, MatrixType::Options,
MatrixType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime>
y(-m_lu.corner(TopRight, m_rank, dimker));
m_lu.corner(TopLeft, m_rank, m_rank)
.template triangularView<UpperTriangular>().solveInPlace(y);
for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = y.row(i);
for(int i = m_rank; i < cols; ++i) result->row(m_q.coeff(i)).setZero();
for(int k = 0; k < dimker; ++k) result->coeffRef(m_q.coeff(m_rank+k), k) = Scalar(1);
}
template<typename MatrixType>
const typename LU<MatrixType>::KernelResultType
LU<MatrixType>::kernel() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
KernelResultType result(m_lu.cols(), dimensionOfKernel());
computeKernel(&result);
return result;
}
template<typename MatrixType>
template<typename ImageMatrixType>
void LU<MatrixType>::computeImage(ImageMatrixType *result) const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
result->resize(m_originalMatrix->rows(), m_rank);
for(int i = 0; i < m_rank; ++i)
result->col(i) = m_originalMatrix->col(m_q.coeff(i));
}
template<typename MatrixType>
const typename LU<MatrixType>::ImageResultType
LU<MatrixType>::image() const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
ImageResultType result(m_originalMatrix->rows(), m_rank);
computeImage(&result);
return result;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool LU<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
{
ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
result->resize(m_lu.cols(), b.cols());
if(m_rank==0)
{
if(b.squaredNorm() == RealScalar(0))
{
result->setZero();
return true;
}
else 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.
* Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
* Step 3: replace c by the solution x to Ux = c. Check if a solution really exists.
* Step 4: result = Qc;
*/
const int rows = m_lu.rows(), cols = m_lu.cols();
ei_assert(b.rows() == rows);
const int smalldim = std::min(rows, cols);
typename OtherDerived::PlainMatrixType c(b.rows(), b.cols());
// Step 1
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>()
.solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols)
{
c.corner(Eigen::BottomLeft, rows-cols, c.cols())
-= m_lu.corner(Eigen::BottomLeft, rows-cols, cols) * c.corner(Eigen::TopLeft, cols, c.cols());
}
// Step 3
if(!isSurjective())
{
// is c is in the image of U ?
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>()
.solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
// Step 4
for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = c.row(i);
for(int i = m_rank; i < m_lu.cols(); ++i) result->row(m_q.coeff(i)).setZero();
return true;
}
/** \lu_module
*
* \return the LU decomposition of \c *this.
*
* \sa class LU
*/
template<typename Derived>
inline const LU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::lu() const
{
return LU<PlainMatrixType>(eval());
}
#endif // EIGEN_LU_H

View File

@ -28,7 +28,7 @@
/** \ingroup LU_Module
*
* \class PartialLU
* \class PartialPivLU
*
* \brief LU decomposition of a matrix with partial pivoting, and related features
*
@ -38,27 +38,30 @@
* is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P
* is a permutation matrix.
*
* Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible matrices.
* So in this class, we plainly require that and take advantage of that to do some simplifications and optimizations.
* This class will assert that the matrix is square, but it won't (actually it can't) check that the matrix is invertible:
* it is your task to check that you only use this decomposition on invertible matrices.
* Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible
* matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class
* does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the
* matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
*
* The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided by class LU.
* The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided
* by class FullPivLU.
*
* This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class,
* such as rank computation. If you need these features, use class LU.
* such as rank computation. If you need these features, use class FullPivLU.
*
* This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses. On the other hand,
* it is \b not suitable to determine whether a given matrix is invertible.
* This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
* in the general case.
* On the other hand, it is \b not suitable to determine whether a given matrix is invertible.
*
* The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
*
* \sa MatrixBase::partialLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class LU
* \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
*/
template<typename MatrixType> class PartialLU
template<typename _MatrixType> class PartialPivLU
{
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
@ -75,63 +78,81 @@ template<typename MatrixType> class PartialLU
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via PartialLU::compute(const MatrixType&).
* perform decompositions via PartialPivLU::compute(const MatrixType&).
*/
PartialLU();
PartialPivLU();
/** Constructor.
*
* \param matrix the matrix of which to compute the LU decomposition.
*
* \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
* If you need to deal with non-full rank, use class LU instead.
* If you need to deal with non-full rank, use class FullPivLU instead.
*/
PartialLU(const MatrixType& matrix);
PartialPivLU(const MatrixType& matrix);
PartialLU& compute(const MatrixType& matrix);
PartialPivLU& compute(const MatrixType& matrix);
/** \returns the LU decomposition matrix: the upper-triangular part is U, the
* unit-lower-triangular part is L (at least for square matrices; in the non-square
* case, special care is needed, see the documentation of class LU).
* case, special care is needed, see the documentation of class FullPivLU).
*
* \sa matrixL(), matrixU()
*/
inline const MatrixType& matrixLU() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return m_lu;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class LU.
* see the examples given in the documentation of class FullPivLU.
*/
inline const IntColVectorType& permutationP() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return m_p;
}
/** This method finds the solution x to the equation Ax=b, where A is the matrix of which
* *this is the LU decomposition. Since if this partial pivoting decomposition the matrix is assumed
* to have full rank, such a solution is assumed to exist and to be unique.
*
* \warning Again, if your matrix may not have full rank, use class LU instead. See LU::solve().
/** This method returns the solution x to the equation Ax=b, where A is the matrix of which
* *this is the LU decomposition.
*
* \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
* the only requirement in order for the equation to make sense is that
* b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
* \param result a pointer to the vector or 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.
*
* Example: \include PartialLU_solve.cpp
* Output: \verbinclude PartialLU_solve.out
* \returns the solution.
*
* Example: \include PartialPivLU_solve.cpp
* Output: \verbinclude PartialPivLU_solve.out
*
* Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
* theoretically exists and is unique regardless of b.
*
* \sa TriangularView::solve(), inverse(), computeInverse()
*/
template<typename OtherDerived, typename ResultType>
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
template<typename Rhs>
inline const ei_solve_retval<PartialPivLU, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return ei_solve_retval<PartialPivLU, Rhs>(*this, b.derived());
}
/** \returns the inverse of the matrix of which *this is the LU decomposition.
*
* \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
* invertibility, use class FullPivLU instead.
*
* \sa MatrixBase::inverse(), LU::inverse()
*/
inline const ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
{
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> >
(*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
}
/** \returns the determinant of the matrix of which
* *this is the LU decomposition. It has only linear complexity
@ -148,33 +169,8 @@ template<typename MatrixType> class PartialLU
*/
typename ei_traits<MatrixType>::Scalar determinant() const;
/** Computes the inverse of the matrix of which *this is the LU decomposition.
*
* \param result a pointer to the matrix into which to store the inverse. Resized if needed.
*
* \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
* invertibility, use class LU instead.
*
* \sa MatrixBase::computeInverse(), inverse()
*/
inline void computeInverse(MatrixType *result) const
{
solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
}
/** \returns the inverse of the matrix of which *this is the LU decomposition.
*
* \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
* invertibility, use class LU instead.
*
* \sa computeInverse(), MatrixBase::inverse()
*/
inline MatrixType inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
}
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
protected:
MatrixType m_lu;
@ -184,7 +180,7 @@ template<typename MatrixType> class PartialLU
};
template<typename MatrixType>
PartialLU<MatrixType>::PartialLU()
PartialPivLU<MatrixType>::PartialPivLU()
: m_lu(),
m_p(),
m_det_p(0),
@ -193,7 +189,7 @@ PartialLU<MatrixType>::PartialLU()
}
template<typename MatrixType>
PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
: m_lu(),
m_p(),
m_det_p(0),
@ -202,9 +198,7 @@ PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
compute(matrix);
}
/** This is the blocked version of ei_lu_unblocked() */
/** This is the blocked version of ei_fullpivlu_unblocked() */
template<typename Scalar, int StorageOrder>
struct ei_partial_lu_impl
{
@ -216,6 +210,7 @@ struct ei_partial_lu_impl
typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
typedef typename MatrixType::RealScalar RealScalar;
/** \internal performs the LU decomposition in-place of the matrix \a lu
* using an unblocked algorithm.
@ -224,8 +219,12 @@ struct ei_partial_lu_impl
* vector \a row_transpositions which must have a size equal to the number
* of columns of the matrix \a lu, and an integer \a nb_transpositions
* which returns the actual number of transpositions.
*
* \returns false if some pivot is exactly zero, in which case the matrix is left with
* undefined coefficients (to avoid generating inf/nan values). Returns true
* otherwise.
*/
static void unblocked_lu(MatrixType& lu, int* row_transpositions, int& nb_transpositions)
static bool unblocked_lu(MatrixType& lu, int* row_transpositions, int& nb_transpositions)
{
const int rows = lu.rows();
const int size = std::min(lu.rows(),lu.cols());
@ -233,9 +232,22 @@ struct ei_partial_lu_impl
for(int k = 0; k < size; ++k)
{
int row_of_biggest_in_col;
lu.col(k).end(rows-k).cwise().abs().maxCoeff(&row_of_biggest_in_col);
RealScalar biggest_in_corner
= lu.col(k).end(rows-k).cwise().abs().maxCoeff(&row_of_biggest_in_col);
row_of_biggest_in_col += k;
if(biggest_in_corner == 0) // the pivot is exactly zero: the matrix is singular
{
// end quickly, avoid generating inf/nan values. Although in this unblocked_lu case
// the result is still valid, there's no need to boast about it because
// the blocked_lu code can't guarantee the same.
// before exiting, make sure to initialize the still uninitialized row_transpositions
// in a sane state without destroying what we already have.
for(int i = k; i < size; i++)
row_transpositions[i] = i;
return false;
}
row_transpositions[k] = row_of_biggest_in_col;
if(k != row_of_biggest_in_col)
@ -252,6 +264,7 @@ struct ei_partial_lu_impl
lu.corner(BottomRight,rrows,rsize).noalias() -= lu.col(k).end(rrows) * lu.row(k).end(rsize);
}
}
return true;
}
/** \internal performs the LU decomposition in-place of the matrix represented
@ -263,11 +276,15 @@ struct ei_partial_lu_impl
* of columns of the matrix \a lu, and an integer \a nb_transpositions
* which returns the actual number of transpositions.
*
* \returns false if some pivot is exactly zero, in which case the matrix is left with
* undefined coefficients (to avoid generating inf/nan values). Returns true
* otherwise.
*
* \note This very low level interface using pointers, etc. is to:
* 1 - reduce the number of instanciations to the strict minimum
* 2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > >
*/
static void blocked_lu(int rows, int cols, Scalar* lu_data, int luStride, int* row_transpositions, int& nb_transpositions, int maxBlockSize=256)
static bool blocked_lu(int rows, int cols, Scalar* lu_data, int luStride, int* row_transpositions, int& nb_transpositions, int maxBlockSize=256)
{
MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
MatrixType lu(lu1,0,0,rows,cols);
@ -277,8 +294,7 @@ struct ei_partial_lu_impl
// if the matrix is too small, no blocking:
if(size<=16)
{
unblocked_lu(lu, row_transpositions, nb_transpositions);
return;
return unblocked_lu(lu, row_transpositions, nb_transpositions);
}
// automatically adjust the number of subdivisions to the size
@ -311,12 +327,20 @@ struct ei_partial_lu_impl
int nb_transpositions_in_panel;
// recursively calls the blocked LU algorithm with a very small
// blocking size:
blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
row_transpositions+k, nb_transpositions_in_panel, 16);
if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
row_transpositions+k, nb_transpositions_in_panel, 16))
{
// end quickly with undefined coefficients, just avoid generating inf/nan values.
// before exiting, make sure to initialize the still uninitialized row_transpositions
// in a sane state without destroying what we already have.
for(int i=k; i<size; ++i)
row_transpositions[i] = i;
return false;
}
nb_transpositions += nb_transpositions_in_panel;
// update permutations and apply them to A10
for(int i=k;i<k+bs; ++i)
for(int i=k; i<k+bs; ++i)
{
int piv = (row_transpositions[i] += k);
A_0.row(i).swap(A_0.row(piv));
@ -334,6 +358,7 @@ struct ei_partial_lu_impl
A22 -= A21 * A12;
}
}
return true;
}
};
@ -351,12 +376,12 @@ void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& n
}
template<typename MatrixType>
PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
{
m_lu = matrix;
m_p.resize(matrix.rows());
ei_assert(matrix.rows() == matrix.cols() && "PartialLU is only for square (and moreover invertible) matrices");
ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
const int size = matrix.rows();
IntColVectorType rows_transpositions(size);
@ -374,54 +399,73 @@ PartialLU<MatrixType>& PartialLU<MatrixType>::compute(const MatrixType& matrix)
}
template<typename MatrixType>
typename ei_traits<MatrixType>::Scalar PartialLU<MatrixType>::determinant() const
typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
return Scalar(m_det_p) * m_lu.diagonal().prod();
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
void PartialLU<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
/***** Implementation of solve() *****************************************************/
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs>
: ei_solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "PartialLU is not initialized.");
EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
/* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
* So we proceed as follows:
* Step 1: compute c = Pb.
* Step 2: replace c by the solution x to Lx = c.
* Step 3: replace c by the solution x to Ux = c.
*/
const int size = dec().matrixLU().rows();
ei_assert(rhs().rows() == size);
/* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
* So we proceed as follows:
* Step 1: compute c = Pb.
* Step 2: replace c by the solution x to Lx = c.
* Step 3: replace c by the solution x to Ux = c.
*/
dst.resize(size, rhs().cols());
const int size = m_lu.rows();
ei_assert(b.rows() == size);
// Step 1
for(int i = 0; i < size; ++i) dst.row(dec().permutationP().coeff(i)) = rhs().row(i);
result->resize(size, b.cols());
// Step 2
dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst);
// Step 1
for(int i = 0; i < size; ++i) result->row(m_p.coeff(i)) = b.row(i);
// Step 3
dec().matrixLU().template triangularView<UpperTriangular>().solveInPlace(dst);
}
};
// Step 2
m_lu.template triangularView<UnitLowerTriangular>().solveInPlace(*result);
/******** MatrixBase methods *******/
// Step 3
m_lu.template triangularView<UpperTriangular>().solveInPlace(*result);
/** \lu_module
*
* \return the partial-pivoting LU decomposition of \c *this.
*
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::partialPivLu() const
{
return PartialPivLU<PlainMatrixType>(eval());
}
/** \lu_module
*
* \return the LU decomposition of \c *this.
* Synonym of partialPivLu().
*
* \sa class LU
* \return the partial-pivoting LU decomposition of \c *this.
*
* \sa class PartialPivLU
*/
template<typename Derived>
inline const PartialLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::partialLu() const
inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::lu() const
{
return PartialLU<PlainMatrixType>(eval());
return PartialPivLU<PlainMatrixType>(eval());
}
#endif // EIGEN_PARTIALLU_H

View File

@ -29,7 +29,7 @@
/** \ingroup QR_Module
* \nonstableyet
*
* \class ColPivotingHouseholderQR
* \class ColPivHouseholderQR
*
* \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
*
@ -38,21 +38,21 @@
* This class performs a rank-revealing QR decomposition using Householder transformations.
*
* This decomposition performs column pivoting in order to be rank-revealing and improve
* numerical stability. It is slower than HouseholderQR, and faster than FullPivotingHouseholderQR.
* numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
*
* \sa MatrixBase::colPivotingHouseholderQr()
* \sa MatrixBase::colPivHouseholderQr()
*/
template<typename MatrixType> class ColPivotingHouseholderQR
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;
@ -68,11 +68,11 @@ template<typename MatrixType> class ColPivotingHouseholderQR
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ColPivotingHouseholderQR::compute(const MatrixType&).
* perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
*/
ColPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
ColPivHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
ColPivotingHouseholderQR(const MatrixType& matrix)
ColPivHouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
m_isInitialized(false)
@ -83,22 +83,27 @@ template<typename MatrixType> class ColPivotingHouseholderQR
/** 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.
*
* Example: \include ColPivotingHouseholderQR_solve.cpp
* Output: \verbinclude ColPivotingHouseholderQR_solve.out
* \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_retval<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return ei_solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
HouseholderSequenceType matrixQ(void) const;
@ -106,15 +111,15 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
const MatrixType& matrixQR() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_qr;
}
ColPivotingHouseholderQR& compute(const MatrixType& matrix);
ColPivHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_cols_permutation;
}
@ -154,7 +159,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
inline int rank() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_rank;
}
@ -165,7 +170,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
inline int dimensionOfKernel() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_qr.cols() - m_rank;
}
@ -177,7 +182,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
inline bool isInjective() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_rank == m_qr.cols();
}
@ -189,7 +194,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
inline bool isSurjective() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_rank == m_qr.rows();
}
@ -200,40 +205,28 @@ template<typename MatrixType> class ColPivotingHouseholderQR
*/
inline bool isInvertible() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR 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 && "ColPivotingHouseholderQR 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_retval<ColPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return ei_solve_retval<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;
@ -247,23 +240,23 @@ template<typename MatrixType> class ColPivotingHouseholderQR
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
typename MatrixType::RealScalar ColPivotingHouseholderQR<MatrixType>::absDeterminant() const
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR 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>
typename MatrixType::RealScalar ColPivotingHouseholderQR<MatrixType>::logAbsDeterminant() const
typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
return m_qr.diagonal().cwise().abs().cwise().log().sum();
}
template<typename MatrixType>
ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
int rows = matrix.rows();
int cols = matrix.cols();
@ -331,56 +324,60 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool ColPivotingHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
: ei_solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
result->resize(m_qr.cols(), b.cols());
if(m_rank==0)
EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
if(b.squaredNorm() == RealScalar(0))
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>
typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
}
@ -388,13 +385,13 @@ typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivoti
/** \return the column-pivoting Householder QR decomposition of \c *this.
*
* \sa class ColPivotingHouseholderQR
* \sa class ColPivHouseholderQR
*/
template<typename Derived>
const ColPivotingHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::colPivotingHouseholderQr() const
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::colPivHouseholderQr() const
{
return ColPivotingHouseholderQR<PlainMatrixType>(eval());
return ColPivHouseholderQR<PlainMatrixType>(eval());
}

View File

@ -29,7 +29,7 @@
/** \ingroup QR_Module
* \nonstableyet
*
* \class FullPivotingHouseholderQR
* \class FullPivHouseholderQR
*
* \brief Householder rank-revealing QR decomposition of a matrix with full pivoting
*
@ -38,21 +38,21 @@
* This class performs a rank-revealing QR decomposition using Householder transformations.
*
* This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal
* numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivotingHouseholderQR.
* numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR.
*
* \sa MatrixBase::fullPivotingHouseholderQr()
* \sa MatrixBase::fullPivHouseholderQr()
*/
template<typename MatrixType> class FullPivotingHouseholderQR
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;
@ -65,11 +65,11 @@ template<typename MatrixType> class FullPivotingHouseholderQR
/** \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&).
* perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
*/
FullPivotingHouseholderQR() : m_isInitialized(false) {}
FullPivHouseholderQR() : m_isInitialized(false) {}
FullPivotingHouseholderQR(const MatrixType& matrix)
FullPivHouseholderQR(const MatrixType& matrix)
: m_isInitialized(false)
{
compute(matrix);
@ -78,22 +78,27 @@ template<typename MatrixType> class FullPivotingHouseholderQR
/** 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.
*
* Example: \include FullPivotingHouseholderQR_solve.cpp
* Output: \verbinclude FullPivotingHouseholderQR_solve.out
* \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_retval<FullPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return ei_solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
MatrixQType matrixQ(void) const;
@ -101,21 +106,21 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
const MatrixType& matrixQR() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_qr;
}
FullPivotingHouseholderQR& compute(const MatrixType& matrix);
FullPivHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_cols_permutation;
}
const IntColVectorType& rowsTranspositions() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_rows_transpositions;
}
@ -155,7 +160,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
inline int rank() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_rank;
}
@ -166,7 +171,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
inline int dimensionOfKernel() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_qr.cols() - m_rank;
}
@ -178,7 +183,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
inline bool isInjective() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_rank == m_qr.cols();
}
@ -190,7 +195,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
inline bool isSurjective() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return m_rank == m_qr.rows();
}
@ -201,40 +206,27 @@ template<typename MatrixType> class FullPivotingHouseholderQR
*/
inline bool isInvertible() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR 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
*/ inline const
ei_solve_retval<FullPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
inverse() const
{
MatrixType result;
computeInverse(&result);
return result;
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return ei_solve_retval<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;
@ -249,23 +241,23 @@ template<typename MatrixType> class FullPivotingHouseholderQR
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::absDeterminant() const
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR 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>
typename MatrixType::RealScalar FullPivotingHouseholderQR<MatrixType>::logAbsDeterminant() const
typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
return m_qr.diagonal().cwise().abs().cwise().log().sum();
}
template<typename MatrixType>
FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
int rows = matrix.rows();
int cols = matrix.cols();
@ -340,62 +332,63 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
return *this;
}
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool FullPivotingHouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
: ei_solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
result->resize(m_qr.cols(), b.cols());
if(m_rank==0)
EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
if(b.squaredNorm() == RealScalar(0))
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>
typename FullPivotingHouseholderQR<MatrixType>::MatrixQType FullPivotingHouseholderQR<MatrixType>::matrixQ() const
typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
ei_assert(m_isInitialized && "FullPivHouseholderQR 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), ...]
@ -417,13 +410,13 @@ typename FullPivotingHouseholderQR<MatrixType>::MatrixQType FullPivotingHousehol
/** \return the full-pivoting Householder QR decomposition of \c *this.
*
* \sa class FullPivotingHouseholderQR
* \sa class FullPivHouseholderQR
*/
template<typename Derived>
const FullPivotingHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::fullPivotingHouseholderQr() const
const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::fullPivHouseholderQr() const
{
return FullPivotingHouseholderQR<PlainMatrixType>(eval());
return FullPivHouseholderQR<PlainMatrixType>(eval());
}
#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H

View File

@ -39,24 +39,24 @@
* stored in a compact way compatible with LAPACK.
*
* Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
* If you want that feature, use FullPivotingHouseholderQR or ColPivotingHouseholderQR instead.
* If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead.
*
* This Householder QR decomposition is faster, but less numerically stable and less feature-full than
* FullPivotingHouseholderQR or ColPivotingHouseholderQR.
* FullPivHouseholderQR or ColPivHouseholderQR.
*
* \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_retval<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return ei_solve_retval<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>
struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
: ei_solve_retval_base<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);
EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
typename OtherDerived::PlainMatrixType c(b);
template<typename Dest> void evalTo(Dest& dst) const
{
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);
// 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());
typename Rhs::PlainMatrixType c(rhs());
m_qr.corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, rank, c.cols()));
// 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()
);
result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
}
dec().matrixQR()
.corner(TopLeft, rank, rank)
.template triangularView<UpperTriangular>()
.solveInPlace(c.corner(TopLeft, rank, c.cols()));
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>

View File

@ -233,7 +233,7 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true>
int diagSize = cols;
if(rows > cols)
{
FullPivotingHouseholderQR<MatrixType> qr(matrix);
FullPivHouseholderQR<MatrixType> qr(matrix);
work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>();
if(ComputeU) svd.m_matrixU = qr.matrixQ();
if(ComputeV)
@ -278,7 +278,7 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true>
typedef Matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime,
MatrixOptions,MaxColsAtCompileTime,MaxRowsAtCompileTime>
TransposeTypeWithSameStorageOrder;
FullPivotingHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint();
if(ComputeV) svd.m_matrixV = qr.matrixQ();
if(ComputeU)

View File

@ -25,6 +25,8 @@
#ifndef EIGEN_SVD_H
#define EIGEN_SVD_H
template<typename MatrixType, typename Rhs> struct ei_svd_solve_impl;
/** \ingroup SVD_Module
* \nonstableyet
*
@ -38,26 +40,27 @@
*
* \sa MatrixBase::SVD()
*/
template<typename MatrixType> class SVD
template<typename _MatrixType> class SVD
{
private:
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
PacketSize = ei_packet_traits<Scalar>::size,
AlignmentMask = int(PacketSize)-1,
MinSize = EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
MinSize = EIGEN_ENUM_MIN(RowsAtCompileTime, ColsAtCompileTime)
};
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVector;
typedef Matrix<Scalar, ColsAtCompileTime, 1> RowVector;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> SingularValuesType;
public:
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
typedef Matrix<Scalar, ColsAtCompileTime, 1> SingularValuesType;
/**
* \brief Default Constructor.
@ -76,8 +79,23 @@ template<typename MatrixType> class SVD
compute(matrix);
}
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
/** \returns a solution of \f$ A x = b \f$ using the current SVD decomposition of A.
*
* \param b the right-hand-side of the equation to solve.
*
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
*
* \sa MatrixBase::svd(),
*/
template<typename Rhs>
inline const ei_solve_retval<SVD, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return ei_solve_retval<SVD, Rhs>(*this, b.derived());
}
const MatrixUType& matrixU() const
{
@ -108,6 +126,18 @@ template<typename MatrixType> class SVD
template<typename ScalingType, typename RotationType>
void computeScalingRotation(ScalingType *positive, RotationType *unitary) const;
inline int rows() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_rows;
}
inline int cols() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_cols;
}
protected:
// Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow.
inline static Scalar pythag(Scalar a, Scalar b)
@ -133,6 +163,7 @@ template<typename MatrixType> class SVD
/** \internal */
SingularValuesType m_sigma;
bool m_isInitialized;
int m_rows, m_cols;
};
/** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
@ -144,8 +175,8 @@ template<typename MatrixType> class SVD
template<typename MatrixType>
SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
{
const int m = matrix.rows();
const int n = matrix.cols();
const int m = m_rows = matrix.rows();
const int n = m_cols = matrix.cols();
m_matU.resize(m, m);
m_matU.setZero();
@ -397,43 +428,35 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
return *this;
}
/** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.
* The parts of the solution corresponding to zero singular values are ignored.
*
* \sa MatrixBase::svd(), LU::solve(), LLT::solve()
*/
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* result) const
template<typename _MatrixType, typename Rhs>
struct ei_solve_retval<SVD<_MatrixType>, Rhs>
: ei_solve_retval_base<SVD<_MatrixType>, Rhs>
{
ei_assert(m_isInitialized && "SVD is not initialized.");
EIGEN_MAKE_SOLVE_HELPERS(SVD<_MatrixType>,Rhs)
const int rows = m_matU.rows();
ei_assert(b.rows() == rows);
result->resize(m_matV.rows(), b.cols());
Scalar maxVal = m_sigma.cwise().abs().maxCoeff();
for (int j=0; j<b.cols(); ++j)
template<typename Dest> void evalTo(Dest& dst) const
{
Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j);
ei_assert(rhs().rows() == dec().rows());
for (int i = 0; i <m_matU.cols(); ++i)
for (int j=0; j<cols(); ++j)
{
Scalar si = m_sigma.coeff(i);
if (ei_isMuchSmallerThan(ei_abs(si),maxVal))
aux.coeffRef(i) = 0;
else
aux.coeffRef(i) /= si;
Matrix<Scalar,MatrixType::RowsAtCompileTime,1> aux = dec().matrixU().adjoint() * rhs().col(j);
for (int i = 0; i < dec().rows(); ++i)
{
Scalar si = dec().singularValues().coeff(i);
if(si == RealScalar(0))
aux.coeffRef(i) = Scalar(0);
else
aux.coeffRef(i) /= si;
}
const int minsize = std::min(dec().rows(),dec().cols());
dst.col(j).start(minsize) = aux.start(minsize);
if(dec().cols()>dec().rows()) dst.col(j).end(cols()-minsize).setZero();
dst.col(j) = dec().matrixV() * dst.col(j);
}
const int cols = m_matV.rows();
const int minsize = std::min(rows,cols);
result->col(j).start(minsize) = aux.start(minsize);
if(cols>rows) result->col(j).end(cols-minsize).setZero();
result->col(j) = m_matV * result->col(j);
}
return true;
}
};
/** Computes the polar decomposition of the matrix, as a product unitary x positive.
*

View File

@ -39,7 +39,7 @@ enum {
*
* \param MatrixType the type of the matrix of which we are computing the LU factorization
*
* \sa class LU, class SparseLLT
* \sa class FullPivLU, class SparseLLT
*/
template<typename MatrixType, int Backend = DefaultBackend>
class SparseLU

View File

@ -0,0 +1,6 @@
FILE(GLOB Eigen_misc_SRCS "*.h")
INSTALL(FILES
${Eigen_misc_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/misc COMPONENT Devel
)

88
Eigen/src/misc/Image.h Normal file
View File

@ -0,0 +1,88 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_MISC_IMAGE_H
#define EIGEN_MISC_IMAGE_H
/** \class ei_image_retval_base
*
*/
template<typename DecompositionType>
struct ei_traits<ei_image_retval_base<DecompositionType> >
{
typedef typename DecompositionType::MatrixType MatrixType;
typedef Matrix<
typename MatrixType::Scalar,
MatrixType::RowsAtCompileTime, // the image is a subspace of the destination space, whose
// dimension is the number of rows of the original matrix
Dynamic, // we don't know at compile time the dimension of the image (the rank)
MatrixType::Options,
MatrixType::MaxRowsAtCompileTime, // the image matrix will consist of columns from the original matrix,
MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns.
> ReturnMatrixType;
};
template<typename _DecompositionType> struct ei_image_retval_base
: public ReturnByValue<ei_image_retval_base<_DecompositionType> >
{
typedef _DecompositionType DecompositionType;
typedef typename DecompositionType::MatrixType MatrixType;
ei_image_retval_base(const DecompositionType& dec, const MatrixType& originalMatrix)
: m_dec(dec), m_rank(dec.rank()),
m_cols(m_rank == 0 ? 1 : m_rank),
m_originalMatrix(originalMatrix)
{}
inline int rows() const { return m_dec.rows(); }
inline int cols() const { return m_cols; }
inline int rank() const { return m_rank; }
inline const DecompositionType& dec() const { return m_dec; }
inline const MatrixType& originalMatrix() const { return m_originalMatrix; }
template<typename Dest> inline void evalTo(Dest& dst) const
{
static_cast<const ei_image_retval<DecompositionType>*>(this)->evalTo(dst);
}
protected:
const DecompositionType& m_dec;
int m_rank, m_cols;
const MatrixType& m_originalMatrix;
};
#define EIGEN_MAKE_IMAGE_HELPERS(DecompositionType) \
typedef typename DecompositionType::MatrixType MatrixType; \
typedef typename MatrixType::Scalar Scalar; \
typedef typename MatrixType::RealScalar RealScalar; \
typedef ei_image_retval_base<DecompositionType> Base; \
using Base::dec; \
using Base::originalMatrix; \
using Base::rank; \
using Base::rows; \
using Base::cols; \
ei_image_retval(const DecompositionType& dec, const MatrixType& originalMatrix) \
: Base(dec, originalMatrix) {}
#endif // EIGEN_MISC_IMAGE_H

85
Eigen/src/misc/Kernel.h Normal file
View File

@ -0,0 +1,85 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_MISC_KERNEL_H
#define EIGEN_MISC_KERNEL_H
/** \class ei_kernel_retval_base
*
*/
template<typename DecompositionType>
struct ei_traits<ei_kernel_retval_base<DecompositionType> >
{
typedef typename DecompositionType::MatrixType MatrixType;
typedef Matrix<
typename MatrixType::Scalar,
MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix"
// is the number of cols of the original matrix
// so that the product "matrix * kernel = zero" makes sense
Dynamic, // we don't know at compile-time the dimension of the kernel
MatrixType::Options,
MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter
MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space,
// whose dimension is the number of columns of the original matrix
> ReturnMatrixType;
};
template<typename _DecompositionType> struct ei_kernel_retval_base
: public ReturnByValue<ei_kernel_retval_base<_DecompositionType> >
{
typedef _DecompositionType DecompositionType;
ei_kernel_retval_base(const DecompositionType& dec)
: m_dec(dec),
m_rank(dec.rank()),
m_cols(m_rank==dec.cols() ? 1 : dec.cols() - m_rank)
{}
inline int rows() const { return m_dec.cols(); }
inline int cols() const { return m_cols; }
inline int rank() const { return m_rank; }
inline const DecompositionType& dec() const { return m_dec; }
template<typename Dest> inline void evalTo(Dest& dst) const
{
static_cast<const ei_kernel_retval<DecompositionType>*>(this)->evalTo(dst);
}
protected:
const DecompositionType& m_dec;
int m_rank, m_cols;
};
#define EIGEN_MAKE_KERNEL_HELPERS(DecompositionType) \
typedef typename DecompositionType::MatrixType MatrixType; \
typedef typename MatrixType::Scalar Scalar; \
typedef typename MatrixType::RealScalar RealScalar; \
typedef ei_kernel_retval_base<DecompositionType> Base; \
using Base::dec; \
using Base::rank; \
using Base::rows; \
using Base::cols; \
ei_kernel_retval(const DecompositionType& dec) : Base(dec) {}
#endif // EIGEN_MISC_KERNEL_H

80
Eigen/src/misc/Solve.h Normal file
View File

@ -0,0 +1,80 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_MISC_SOLVE_H
#define EIGEN_MISC_SOLVE_H
/** \class ei_solve_retval_base
*
*/
template<typename DecompositionType, typename Rhs>
struct ei_traits<ei_solve_retval_base<DecompositionType, Rhs> >
{
typedef typename DecompositionType::MatrixType MatrixType;
typedef Matrix<typename Rhs::Scalar,
MatrixType::ColsAtCompileTime,
Rhs::ColsAtCompileTime,
Rhs::PlainMatrixType::Options,
MatrixType::MaxColsAtCompileTime,
Rhs::MaxColsAtCompileTime> ReturnMatrixType;
};
template<typename _DecompositionType, typename Rhs> struct ei_solve_retval_base
: public ReturnByValue<ei_solve_retval_base<_DecompositionType, Rhs> >
{
typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNestedCleaned;
typedef _DecompositionType DecompositionType;
ei_solve_retval_base(const DecompositionType& dec, const Rhs& rhs)
: m_dec(dec), m_rhs(rhs)
{}
inline int rows() const { return m_dec.cols(); }
inline int cols() const { return m_rhs.cols(); }
inline const DecompositionType& dec() const { return m_dec; }
inline const RhsNestedCleaned& rhs() const { return m_rhs; }
template<typename Dest> inline void evalTo(Dest& dst) const
{
static_cast<const ei_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
}
protected:
const DecompositionType& m_dec;
const typename Rhs::Nested m_rhs;
};
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType,Rhs) \
typedef typename DecompositionType::MatrixType MatrixType; \
typedef typename MatrixType::Scalar Scalar; \
typedef typename MatrixType::RealScalar RealScalar; \
typedef ei_solve_retval_base<DecompositionType,Rhs> Base; \
using Base::dec; \
using Base::rhs; \
using Base::rows; \
using Base::cols; \
ei_solve_retval(const DecompositionType& dec, const Rhs& rhs) \
: Base(dec, rhs) {}
#endif // EIGEN_MISC_SOLVE_H

View File

@ -218,7 +218,7 @@ public :
}
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
C = X.partialLu().matrixLU();
C = X.partialPivLu().matrixLU();
}
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){

View File

@ -98,7 +98,7 @@ int main(int argc, char *argv[])
BenchTimer timer;
timer.start();
LU<DenseMatrix> lu(m1);
FullPivLU<DenseMatrix> lu(m1);
timer.stop();
std::cout << "Eigen/dense:\t" << timer.value() << endl;

View File

@ -1,5 +1,3 @@
option(EIGEN_NO_ASSERTION_CHECKING "Disable checking of assertions" OFF)
# similar to set_target_properties but append the property instead of overwriting it
@ -18,30 +16,11 @@ macro(ei_add_property prop value)
set_property(GLOBAL PROPERTY ${prop} "${previous} ${value}")
endmacro(ei_add_property)
# Macro to add a test
#
# the unique parameter testname must correspond to a file
# <testname>.cpp which follows this pattern:
#
# #include "main.h"
# void test_<testname>() { ... }
#
# this macro add an executable test_<testname> as well as a ctest test
# named <testname>.
#
# it also adds another executable debug_<testname> that compiles in full debug mode
# and is not added to the test target. The idea is that when a test fails you want
# a quick way of rebuilding this specific test in full debug mode.
#
# On platforms with bash simply run:
# "ctest -V" or "ctest -V -R <testname>"
# On other platform use ctest as usual
#
macro(ei_add_test testname)
set(targetname test_${testname})
#internal. See documentation of ei_add_test for details.
macro(ei_add_test_internal testname testname_with_suffix)
set(targetname test_${testname_with_suffix})
if(NOT MSVC_IDE)
set(debug_targetname debug_${testname})
set(debug_targetname debug_${testname_with_suffix})
endif(NOT MSVC_IDE)
set(filename ${testname}.cpp)
@ -80,13 +59,21 @@ macro(ei_add_test testname)
endif(NOT EIGEN_NO_ASSERTION_CHECKING)
# let the user pass e.g. optimization flags, but don't apply them to the debug target
if(${ARGC} GREATER 1)
ei_add_target_property(${targetname} COMPILE_FLAGS "${ARGV1}")
endif(${ARGC} GREATER 1)
# let the user pass flags. Note that if the user passes an optimization flag here, it's important that
# we counter it by a no-optimization flag!
if(${ARGC} GREATER 2)
ei_add_target_property(${targetname} COMPILE_FLAGS "${ARGV2}")
if(NOT MSVC_IDE)
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "${ARGV2} ${EI_NO_OPTIMIZATION_FLAG}")
endif(NOT MSVC_IDE)
endif(${ARGC} GREATER 2)
# for the debug target, add full debug options
if(CMAKE_COMPILER_IS_GNUCXX)
# disable debug symbols: i get 2 GB of them!
# the user can still use the debug targets as needed.
# for MSVC, the equivalent (release mode) does the same.
ei_add_target_property(${targetname} COMPILE_FLAGS "-g0")
# O0 is in principle redundant here, but doesn't hurt
ei_add_target_property(${debug_targetname} COMPILE_FLAGS "-O0 -g3")
elseif(MSVC)
@ -101,27 +88,109 @@ macro(ei_add_test testname)
endif(NOT MSVC_IDE)
target_link_libraries(${targetname} ${EXTERNAL_LIBS})
if(${ARGC} GREATER 2)
string(STRIP "${ARGV2}" ARGV2_stripped)
string(LENGTH "${ARGV2_stripped}" ARGV2_stripped_length)
if(${ARGV2_stripped_length} GREATER 0)
target_link_libraries(${targetname} ${ARGV2})
if(${ARGC} GREATER 3)
string(STRIP "${ARGV3}" ARGV3_stripped)
string(LENGTH "${ARGV3_stripped}" ARGV3_stripped_length)
if(${ARGV3_stripped_length} GREATER 0)
target_link_libraries(${targetname} ${ARGV3})
if(NOT MSVC_IDE)
target_link_libraries(${debug_targetname} ${ARGV2})
target_link_libraries(${debug_targetname} ${ARGV3})
endif(NOT MSVC_IDE)
endif(${ARGV2_stripped_length} GREATER 0)
endif(${ARGC} GREATER 2)
endif(${ARGV3_stripped_length} GREATER 0)
endif(${ARGC} GREATER 3)
if(WIN32)
if(CYGWIN)
add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}")
add_test(${testname_with_suffix} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname_with_suffix}")
else(CYGWIN)
add_test(${testname} "${targetname}")
add_test(${testname_with_suffix} "${targetname}")
endif(CYGWIN)
else(WIN32)
add_test(${testname} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname}")
add_test(${testname_with_suffix} "${Eigen_SOURCE_DIR}/test/runtest.sh" "${testname_with_suffix}")
endif(WIN32)
endmacro(ei_add_test_internal)
# Macro to add a test
#
# the unique mandatory parameter testname must correspond to a file
# <testname>.cpp which follows this pattern:
#
# #include "main.h"
# void test_<testname>() { ... }
#
# Depending on the contents of that file, this macro can have 2 behaviors,
# see below.
#
# Optional parameters can be passed: the 2nd parameter is additional compile flags
# (optimization will be explicitly disabled for the debug test targets) and the 3rd
# parameter is libraries to link to.
#
# A. Default behavior
#
# this macro add an executable test_<testname> as well as a ctest test
# named <testname>.
#
# it also adds another executable debug_<testname> that compiles in full debug mode
# and is not added to the test target. The idea is that when a test fails you want
# a quick way of rebuilding this specific test in full debug mode.
#
# On platforms with bash simply run:
# "ctest -V" or "ctest -V -R <testname>"
# On other platform use ctest as usual
#
# B. Multi-part behavior
#
# If the source file matches the regexp
# CALL_SUBTEST(-9]+|EIGEN_TEST_PART_[0-9]+
# then it is interpreted as a multi-part test. The behavior then depends on the
# CMake option EIGEN_SPLIT_LARGE_TESTS, which is ON by default.
#
# If EIGEN_SPLIT_LARGE_TESTS is OFF, the behavior is the same as in A (the multi-part
# aspect is ignored).
#
# If EIGEN_SPLIT_LARGE_TESTS is ON, the test is split into multiple executables
# test_<testname>_<N>
# where N runs from 1 to the greatest occurence found in the source file. Each of these
# executables is built passing -DEIGEN_TEST_PART_N. This allows to split large tests
# into smaller executables.
#
# The same holds for the debug executables.
#
# Moreover, targets test_<testname> and debug_<testname> are still generated, they
# have the effect of building all the parts of the test.
#
# Again, ctest -R allows to run all matching tests.
#
macro(ei_add_test testname)
file(READ "${testname}.cpp" test_source)
set(parts 0)
string(REGEX MATCHALL "CALL_SUBTEST_[0-9]+|EIGEN_TEST_PART_[0-9]+"
occurences "${test_source}")
string(REGEX REPLACE "CALL_SUBTEST_|EIGEN_TEST_PART_" "" suffixes "${occurences}")
list(REMOVE_DUPLICATES suffixes)
if(EIGEN_SPLIT_LARGE_TESTS AND suffixes)
add_custom_target(test_${testname})
if(NOT MSVC_IDE)
add_custom_target(debug_${testname})
endif(NOT MSVC_IDE)
foreach(suffix ${suffixes})
ei_add_test_internal(${testname} ${testname}_${suffix}
"${ARGV1} -DEIGEN_TEST_PART_${suffix}=1" "${ARGV2}")
add_dependencies(test_${testname} test_${testname}_${suffix})
if(NOT MSVC_IDE)
add_dependencies(debug_${testname} debug_${testname}_${suffix})
endif(NOT MSVC_IDE)
endforeach(suffix)
else(EIGEN_SPLIT_LARGE_TESTS AND suffixes)
set(symbols_to_enable_all_parts "")
foreach(suffix ${suffixes})
set(symbols_to_enable_all_parts
"${symbols_to_enable_all_parts} -DEIGEN_TEST_PART_${suffix}=1")
endforeach(suffix)
ei_add_test_internal(${testname} ${testname} "${ARGV1} ${symbols_to_enable_all_parts}" "${ARGV2}")
endif(EIGEN_SPLIT_LARGE_TESTS AND suffixes)
endmacro(ei_add_test)
# print a summary of the different options
@ -204,12 +273,15 @@ if(CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${COVERAGE_FLAGS} -fno-inline-functions")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${COVERAGE_FLAGS} -O0 -g2")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${COVERAGE_FLAGS} -O0 -g3")
endif(CMAKE_SYSTEM_NAME MATCHES Linux)
set(EI_OFLAG "-O2")
set(EI_NO_OPTIMIZATION_FLAG "-O0")
elseif(MSVC)
set(CMAKE_CXX_FLAGS_DEBUG "/D_DEBUG /MDd /Zi /Ob0 /Od" CACHE STRING "Flags used by the compiler during debug builds." FORCE)
set(EI_OFLAG "/O2")
set(EI_NO_OPTIMIZATION_FLAG "/O0")
else(CMAKE_COMPILER_IS_GNUCXX)
set(EI_OFLAG "")
set(EI_NO_OPTIMIZATION_FLAG "")
endif(CMAKE_COMPILER_IS_GNUCXX)

View File

@ -55,8 +55,8 @@ matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1}
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.
MatrixBase::partialPivLu() function. This yields an object of the class PartialPivLU. Then, the
PartialPivLU::solve() method is called to compute a solution.
As an example, suppose we want to solve the following system of linear equations:
@ -69,9 +69,9 @@ As an example, suppose we want to solve the following system of linear equations
The following program solves this system:
<table class="tutorial_code"><tr><td>
\include Tutorial_PartialLU_solve.cpp
\include Tutorial_PartialPivLU_solve.cpp
</td><td>
output: \include Tutorial_PartialLU_solve.out
output: \include Tutorial_PartialPivLU_solve.out
</td></tr></table>
There are many situations in which we want to solve the same system of equations with different
@ -91,7 +91,7 @@ problem, and whether you want to solve it at all, after you solved the first pro
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.
PartialPivLU::PartialPivLU(const MatrixType&) to compute the %LU decomposition.
<table class="tutorial_code"><tr><td>
\include Tutorial_solve_reuse_decomposition.cpp
@ -102,7 +102,7 @@ output: \include Tutorial_solve_reuse_decomposition.out
\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
both cases, PartialPivLU::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 second equation. If we try
@ -114,10 +114,10 @@ to solve this inconsistent system with Eigen, we find:
output: \include Tutorial_solve_singular.out
</td></tr></table>
The %LU decomposition with \b full pivoting (class LU) and the singular value decomposition (class
The %LU decomposition with \b full pivoting (class FullPivLU) 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.
\sa LU_Module, MatrixBase::partialPivLu(), PartialPivLU::solve(), class PartialPivLU.
\subsection TutorialAdvSolvers_Cholesky Cholesky decomposition
@ -228,7 +228,7 @@ Note that the function inverse() is defined in the \ref LU_Module.
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).
decomposition with partial pivoting (PartialPivLU).
The solver based on the %SVD uses the class SVD. It can handle singular matrices. Here is an example
of its use:
@ -245,7 +245,7 @@ svdOfA.solve(b, &x);
\endcode
%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.
partial pivoting. It is defined in the class FullPivLU. The solver can also handle singular matrices.
\code
#include <Eigen/LU>
@ -254,13 +254,13 @@ MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x;
A.lu().solve(b, &x);
LU<MatrixXf> luOfA(A);
FullPivLU<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.
\sa class SVD, SVD::solve(), SVD_Module, class FullPivLU, LU::solve(), LU_Module.
@ -281,7 +281,7 @@ Alternatively, you can construct a named LU decomposition, which allows you to r
\code
#include <Eigen/LU>
MatrixXf A = MatrixXf::Random(20,20);
Eigen::LU<MatrixXf> lu(A);
Eigen::FullPivLU<MatrixXf> lu(A);
cout << "The rank of A is" << lu.rank() << endl;
if(lu.isInvertible()) {
cout << "A is invertible, its inverse is:" << endl << lu.inverse() << endl;
@ -292,7 +292,7 @@ else {
}
\endcode
\sa LU_Module, LU::solve(), class LU
\sa LU_Module, LU::solve(), class FullPivLU
<a href="#" class="top">top</a>\section TutorialAdvCholesky Cholesky
todo

View File

@ -215,7 +215,10 @@ ALIASES = "only_for_vectors=This is only for vectors (either row-
"eigenvalues_module=This is defined in the %Eigenvalues module. \code #include <Eigen/Eigenvalues> \endcode" \
"label=\bug" \
"redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \
"nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\""
"nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \
"note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \
"note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \
"note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values."
# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
# sources only. Doxygen will then generate output that is more tailored for C.

View File

@ -12,7 +12,6 @@ int main(int, char *[])
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);
Vector3f x = A.lu().solve(b);
cout << "The solution is:" << endl << x << endl;
}

View File

@ -6,4 +6,4 @@ cout << "Here is the matrix m:" << endl << m << endl;
cout << "Notice that the middle column is the sum of the two others, so the "
<< "columns are linearly dependent." << endl;
cout << "Here is a matrix whose columns have the same span but are linearly independent:"
<< endl << m.lu().image() << endl;
<< endl << m.fullPivLu().image(m) << endl;

View File

@ -1,6 +1,6 @@
MatrixXf m = MatrixXf::Random(3,5);
cout << "Here is the matrix m:" << endl << m << endl;
MatrixXf ker = m.lu().kernel();
MatrixXf ker = m.fullPivLu().kernel();
cout << "Here is a matrix whose columns form a basis of the kernel of m:"
<< endl << ker << endl;
cout << "By definition of the kernel, m*ker is zero:"

View File

@ -1,15 +1,11 @@
typedef Matrix<float,2,3> Matrix2x3;
typedef Matrix<float,3,2> Matrix3x2;
Matrix2x3 m = Matrix2x3::Random();
Matrix<float,2,3> m = Matrix<float,2,3>::Random();
Matrix2f y = Matrix2f::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is the matrix y:" << endl << y << endl;
Matrix3x2 x;
if(m.lu().solve(y, &x))
Matrix<float,3,2> x = m.fullPivLu().solve(y);
if((m*x).isApprox(y))
{
assert(y.isApprox(m*x));
cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
}
else
cout << "The equation mx=y does not have any solution." << endl;

View File

@ -4,6 +4,6 @@ Matrix3f y = Matrix3f::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is the matrix y:" << endl << y << endl;
Matrix3f x;
m.householderQr().solve(y, &x);
x = m.householderQr().solve(y);
assert(y.isApprox(m*x));
cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;

View File

@ -3,6 +3,6 @@ typedef Matrix<float,Dynamic,2> DataMatrix;
DataMatrix samples = DataMatrix::Random(12,2);
VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
// and let's solve samples * [x y]^T = elevations in least square sense:
Matrix<float,2,1> xy;
(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy);
Matrix<float,2,1> xy
= (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations));
cout << xy << endl;

View File

@ -1,13 +0,0 @@
MatrixXd m(3,3);
m << 1,1,0,
1,3,2,
0,1,1;
cout << "Here is the matrix m:" << endl << m << endl;
LU<Matrix3d> lu(m);
// allocate the matrix img with the correct size to avoid reallocation
MatrixXd img(m.rows(), lu.rank());
lu.computeImage(&img);
cout << "Notice that the middle column is the sum of the two others, so the "
<< "columns are linearly dependent." << endl;
cout << "Here is a matrix whose columns have the same span but are linearly independent:"
<< endl << img << endl;

View File

@ -1,10 +0,0 @@
MatrixXf m = MatrixXf::Random(3,5);
cout << "Here is the matrix m:" << endl << m << endl;
LU<MatrixXf> lu(m);
// allocate the matrix ker with the correct size to avoid reallocation
MatrixXf ker(m.rows(), lu.dimensionOfKernel());
lu.computeKernel(&ker);
cout << "Here is a matrix whose columns form a basis of the kernel of m:"
<< endl << ker << endl;
cout << "By definition of the kernel, m*ker is zero:"
<< endl << m*ker << endl;

View File

@ -1,5 +0,0 @@
Matrix3d m = Matrix3d::Random();
cout << "Here is the matrix m:" << endl << m << endl;
Matrix3d inv;
m.computeInverse(&inv);
cout << "Its inverse is:" << endl << inv << endl;

View File

@ -0,0 +1,13 @@
Matrix3d m = Matrix3d::Random();
cout << "Here is the matrix m:" << endl << m << endl;
Matrix3d inverse;
bool invertible;
double determinant;
m.computeInverseAndDetWithCheck(inverse,determinant,invertible);
cout << "Its determinant is " << determinant << endl;
if(invertible) {
cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
}
else {
cout << "It is not invertible." << endl;
}

View File

@ -1,8 +1,10 @@
Matrix3d m = Matrix3d::Random();
cout << "Here is the matrix m:" << endl << m << endl;
Matrix3d inv;
if(m.computeInverseWithCheck(&inv)) {
cout << "It is invertible, and its inverse is:" << endl << inv << endl;
Matrix3d inverse;
bool invertible;
m.computeInverseWithCheck(inverse,invertible);
if(invertible) {
cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
}
else {
cout << "It is not invertible." << endl;

View File

@ -2,7 +2,6 @@ MatrixXd A = MatrixXd::Random(3,3);
MatrixXd B = MatrixXd::Random(3,2);
cout << "Here is the invertible matrix A:" << endl << A << endl;
cout << "Here is the matrix B:" << endl << B << endl;
MatrixXd X;
if(A.lu().solve(B, &X))
MatrixXd X = A.lu().solve(B);
cout << "Here is the (unique) solution X to the equation AX=B:" << endl << X << endl;
cout << "Relative error: " << (A*X-B).norm() / B.norm() << endl;

View File

@ -3,7 +3,7 @@ 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);
X = A.lu().solve(B);
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;

View File

@ -1,13 +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
PartialPivLU<Matrix3f> luOfA(A); // compute LU decomposition of A
Vector3f b;
b << 3,3,4;
Vector3f x;
luOfA.solve(b, &x);
x = luOfA.solve(b);
cout << "The solution with right-hand side (3,3,4) is:" << endl;
cout << x << endl;
b << 1,1,1;
luOfA.solve(b, &x);
x = luOfA.solve(b);
cout << "The solution with right-hand side (1,1,1) is:" << endl;
cout << x << endl;

View File

@ -5,5 +5,5 @@ 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);
x = A.lu().solve(b);
cout << "The solution is:" << endl << x << endl;

View File

@ -2,7 +2,7 @@ typedef Matrix<double, 5, 3> Matrix5x3;
typedef Matrix<double, 5, 5> Matrix5x5;
Matrix5x3 m = Matrix5x3::Random();
cout << "Here is the matrix m:" << endl << m << endl;
Eigen::LU<Matrix5x3> lu(m);
Eigen::FullPivLU<Matrix5x3> lu(m);
cout << "Here is, up to permutations, its LU decomposition matrix:"
<< endl << lu.matrixLU() << endl;
cout << "Here is the L part:" << endl;

View File

@ -3,6 +3,8 @@ add_custom_target(btest)
include(EigenTesting)
ei_init_testing()
option(EIGEN_SPLIT_LARGE_TESTS "Split large tests into smaller executables" ON)
find_package(GSL)
if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
set(GSL_FOUND "")

View File

@ -108,16 +108,17 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
void test_adjoint()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( adjoint(Matrix<float, 1, 1>()) );
CALL_SUBTEST( adjoint(Matrix3d()) );
CALL_SUBTEST( adjoint(Matrix4f()) );
CALL_SUBTEST( adjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST( adjoint(MatrixXi(8, 12)) );
CALL_SUBTEST( adjoint(MatrixXf(21, 21)) );
CALL_SUBTEST_1( adjoint(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( adjoint(Matrix3d()) );
CALL_SUBTEST_3( adjoint(Matrix4f()) );
CALL_SUBTEST_4( adjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST_5( adjoint(MatrixXi(8, 12)) );
CALL_SUBTEST_6( adjoint(MatrixXf(21, 21)) );
}
// test a large matrix only once
CALL_SUBTEST( adjoint(Matrix<float, 100, 100>()) );
CALL_SUBTEST_7( adjoint(Matrix<float, 100, 100>()) );
#ifdef EIGEN_TEST_PART_4
{
MatrixXcf a(10,10), b(10,10);
VERIFY_RAISES_ASSERT(a = a.transpose());
@ -128,5 +129,6 @@ void test_adjoint()
VERIFY_RAISES_ASSERT(a = a.adjoint() + b);
VERIFY_RAISES_ASSERT(a = b + a.adjoint());
}
#endif
}

View File

@ -146,26 +146,26 @@ template<typename VectorType> void lpNorm(const VectorType& v)
void test_array()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( array(Matrix<float, 1, 1>()) );
CALL_SUBTEST( array(Matrix2f()) );
CALL_SUBTEST( array(Matrix4d()) );
CALL_SUBTEST( array(MatrixXcf(3, 3)) );
CALL_SUBTEST( array(MatrixXf(8, 12)) );
CALL_SUBTEST( array(MatrixXi(8, 12)) );
CALL_SUBTEST_1( array(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( array(Matrix2f()) );
CALL_SUBTEST_3( array(Matrix4d()) );
CALL_SUBTEST_4( array(MatrixXcf(3, 3)) );
CALL_SUBTEST_5( array(MatrixXf(8, 12)) );
CALL_SUBTEST_6( array(MatrixXi(8, 12)) );
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( comparisons(Matrix<float, 1, 1>()) );
CALL_SUBTEST( comparisons(Matrix2f()) );
CALL_SUBTEST( comparisons(Matrix4d()) );
CALL_SUBTEST( comparisons(MatrixXf(8, 12)) );
CALL_SUBTEST( comparisons(MatrixXi(8, 12)) );
CALL_SUBTEST_1( comparisons(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( comparisons(Matrix2f()) );
CALL_SUBTEST_3( comparisons(Matrix4d()) );
CALL_SUBTEST_5( comparisons(MatrixXf(8, 12)) );
CALL_SUBTEST_6( comparisons(MatrixXi(8, 12)) );
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( lpNorm(Matrix<float, 1, 1>()) );
CALL_SUBTEST( lpNorm(Vector2f()) );
CALL_SUBTEST( lpNorm(Vector3d()) );
CALL_SUBTEST( lpNorm(Vector4f()) );
CALL_SUBTEST( lpNorm(VectorXf(16)) );
CALL_SUBTEST( lpNorm(VectorXcd(10)) );
CALL_SUBTEST_1( lpNorm(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( lpNorm(Vector2f()) );
CALL_SUBTEST_7( lpNorm(Vector3d()) );
CALL_SUBTEST_8( lpNorm(Vector4f()) );
CALL_SUBTEST_5( lpNorm(VectorXf(16)) );
CALL_SUBTEST_4( lpNorm(VectorXcf(10)) );
}
}

View File

@ -76,11 +76,11 @@ template<typename MatrixType> void replicate(const MatrixType& m)
void test_array_replicate()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( replicate(Matrix<float, 1, 1>()) );
CALL_SUBTEST( replicate(Vector2f()) );
CALL_SUBTEST( replicate(Vector3d()) );
CALL_SUBTEST( replicate(Vector4f()) );
CALL_SUBTEST( replicate(VectorXf(16)) );
CALL_SUBTEST( replicate(VectorXcd(10)) );
CALL_SUBTEST_1( replicate(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( replicate(Vector2f()) );
CALL_SUBTEST_3( replicate(Vector3d()) );
CALL_SUBTEST_4( replicate(Vector4f()) );
CALL_SUBTEST_5( replicate(VectorXf(16)) );
CALL_SUBTEST_6( replicate(VectorXcd(10)) );
}
}

View File

@ -163,19 +163,20 @@ template<typename MatrixType> void reverse(const MatrixType& m)
void test_array_reverse()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( reverse(Matrix<float, 1, 1>()) );
CALL_SUBTEST( reverse(Matrix2f()) );
CALL_SUBTEST( reverse(Matrix4f()) );
CALL_SUBTEST( reverse(Matrix4d()) );
CALL_SUBTEST( reverse(MatrixXcf(3, 3)) );
CALL_SUBTEST( reverse(MatrixXi(6, 3)) );
CALL_SUBTEST( reverse(MatrixXcd(20, 20)) );
CALL_SUBTEST( reverse(Matrix<float, 100, 100>()) );
CALL_SUBTEST( reverse(Matrix<float,Dynamic,Dynamic,RowMajor>(6,3)) );
CALL_SUBTEST_1( reverse(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( reverse(Matrix2f()) );
CALL_SUBTEST_3( reverse(Matrix4f()) );
CALL_SUBTEST_4( reverse(Matrix4d()) );
CALL_SUBTEST_5( reverse(MatrixXcf(3, 3)) );
CALL_SUBTEST_6( reverse(MatrixXi(6, 3)) );
CALL_SUBTEST_7( reverse(MatrixXcd(20, 20)) );
CALL_SUBTEST_8( reverse(Matrix<float, 100, 100>()) );
CALL_SUBTEST_9( reverse(Matrix<float,Dynamic,Dynamic,RowMajor>(6,3)) );
}
#ifdef EIGEN_TEST_PART_3
Vector4f x; x << 1, 2, 3, 4;
Vector4f y; y << 4, 3, 2, 1;
VERIFY(x.reverse()[1] == 3);
VERIFY(x.reverse() == y);
#endif
}

View File

@ -79,6 +79,6 @@ void test_bandmatrix()
int cols = ei_random<int>(1,10);
int sups = ei_random<int>(0,cols-1);
int subs = ei_random<int>(0,rows-1);
CALL_SUBTEST( bandmatrix(BandMatrix<float>(rows,cols,sups,subs)) );
CALL_SUBTEST(bandmatrix(BandMatrix<float>(rows,cols,sups,subs)) );
}
}

View File

@ -146,17 +146,17 @@ void casting()
void test_basicstuff()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( basicStuff(Matrix<float, 1, 1>()) );
CALL_SUBTEST( basicStuff(Matrix4d()) );
CALL_SUBTEST( basicStuff(MatrixXcf(3, 3)) );
CALL_SUBTEST( basicStuff(MatrixXi(8, 12)) );
CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) );
CALL_SUBTEST( basicStuff(Matrix<float, 100, 100>()) );
CALL_SUBTEST( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
CALL_SUBTEST_1( basicStuff(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( basicStuff(Matrix4d()) );
CALL_SUBTEST_3( basicStuff(MatrixXcf(3, 3)) );
CALL_SUBTEST_4( basicStuff(MatrixXi(8, 12)) );
CALL_SUBTEST_5( basicStuff(MatrixXcd(20, 20)) );
CALL_SUBTEST_6( basicStuff(Matrix<float, 100, 100>()) );
CALL_SUBTEST_7( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) );
CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) );
CALL_SUBTEST_3( basicStuffComplex(MatrixXcf(21, 17)) );
CALL_SUBTEST_5( basicStuffComplex(MatrixXcd(2, 3)) );
}
CALL_SUBTEST(casting());
CALL_SUBTEST_2(casting());
}

View File

@ -22,7 +22,10 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_NO_ASSERTION_CHECKING
#define EIGEN_NO_ASSERTION_CHECKING
#endif
#include "main.h"
#include <Eigen/Cholesky>
#include <Eigen/QR>
@ -93,17 +96,17 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
{
LLT<SquareMatrixType,LowerTriangular> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.matrixL().toDense() * chollo.matrixL().adjoint().toDense());
chollo.solve(vecB, &vecX);
vecX = chollo.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
chollo.solve(matB, &matX);
matX = chollo.solve(matB);
VERIFY_IS_APPROX(symm * matX, matB);
// test the upper mode
LLT<SquareMatrixType,UpperTriangular> cholup(symmUp);
VERIFY_IS_APPROX(symm, cholup.matrixL().toDense() * cholup.matrixL().adjoint().toDense());
cholup.solve(vecB, &vecX);
vecX = cholup.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
cholup.solve(matB, &matX);
matX = cholup.solve(matB);
VERIFY_IS_APPROX(symm * matX, matB);
}
@ -118,9 +121,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
LDLT<SquareMatrixType> ldlt(symm);
// TODO(keir): This doesn't make sense now that LDLT pivots.
//VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
ldlt.solve(vecB, &vecX);
vecX = ldlt.solve(vecB);
VERIFY_IS_APPROX(symm * vecX, vecB);
ldlt.solve(matB, &matX);
matX = ldlt.solve(matB);
VERIFY_IS_APPROX(symm * matX, matB);
}
@ -132,7 +135,7 @@ template<typename MatrixType> void cholesky_verify_assert()
LLT<MatrixType> llt;
VERIFY_RAISES_ASSERT(llt.matrixL())
VERIFY_RAISES_ASSERT(llt.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(llt.solve(tmp))
VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))
LDLT<MatrixType> ldlt;
@ -141,24 +144,24 @@ template<typename MatrixType> void cholesky_verify_assert()
VERIFY_RAISES_ASSERT(ldlt.vectorD())
VERIFY_RAISES_ASSERT(ldlt.isPositive())
VERIFY_RAISES_ASSERT(ldlt.isNegative())
VERIFY_RAISES_ASSERT(ldlt.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(ldlt.solve(tmp))
VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp))
}
void test_cholesky()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
CALL_SUBTEST( cholesky(Matrix2d()) );
CALL_SUBTEST( cholesky(Matrix3f()) );
CALL_SUBTEST( cholesky(Matrix4d()) );
CALL_SUBTEST( cholesky(MatrixXd(200,200)) );
CALL_SUBTEST( cholesky(MatrixXcd(100,100)) );
CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
CALL_SUBTEST_2( cholesky(MatrixXd(1,1)) );
CALL_SUBTEST_3( cholesky(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) );
CALL_SUBTEST_2( cholesky(MatrixXd(200,200)) );
CALL_SUBTEST_6( cholesky(MatrixXcd(100,100)) );
}
CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
}

View File

@ -110,20 +110,20 @@ void run_vector_tests()
void test_conservative_resize()
{
run_matrix_tests<int, Eigen::RowMajor>();
run_matrix_tests<int, Eigen::ColMajor>();
run_matrix_tests<float, Eigen::RowMajor>();
run_matrix_tests<float, Eigen::ColMajor>();
run_matrix_tests<double, Eigen::RowMajor>();
run_matrix_tests<double, Eigen::ColMajor>();
run_matrix_tests<std::complex<float>, Eigen::RowMajor>();
run_matrix_tests<std::complex<float>, Eigen::ColMajor>();
run_matrix_tests<std::complex<double>, Eigen::RowMajor>();
run_matrix_tests<std::complex<double>, Eigen::ColMajor>();
CALL_SUBTEST_1((run_matrix_tests<int, Eigen::RowMajor>()));
CALL_SUBTEST_1((run_matrix_tests<int, Eigen::ColMajor>()));
CALL_SUBTEST_2((run_matrix_tests<float, Eigen::RowMajor>()));
CALL_SUBTEST_2((run_matrix_tests<float, Eigen::ColMajor>()));
CALL_SUBTEST_3((run_matrix_tests<double, Eigen::RowMajor>()));
CALL_SUBTEST_3((run_matrix_tests<double, Eigen::ColMajor>()));
CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::RowMajor>()));
CALL_SUBTEST_4((run_matrix_tests<std::complex<float>, Eigen::ColMajor>()));
CALL_SUBTEST_5((run_matrix_tests<std::complex<double>, Eigen::RowMajor>()));
CALL_SUBTEST_6((run_matrix_tests<std::complex<double>, Eigen::ColMajor>()));
run_vector_tests<int>();
run_vector_tests<float>();
run_vector_tests<double>();
run_vector_tests<std::complex<float> >();
run_vector_tests<std::complex<double> >();
CALL_SUBTEST_1((run_vector_tests<int>()));
CALL_SUBTEST_2((run_vector_tests<float>()));
CALL_SUBTEST_3((run_vector_tests<double>()));
CALL_SUBTEST_4((run_vector_tests<std::complex<float> >()));
CALL_SUBTEST_5((run_vector_tests<std::complex<double> >()));
}

View File

@ -163,11 +163,11 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
void test_cwiseop()
{
for(int i = 0; i < g_repeat ; i++) {
CALL_SUBTEST( cwiseops(Matrix<float, 1, 1>()) );
CALL_SUBTEST( cwiseops(Matrix4d()) );
CALL_SUBTEST( cwiseops(MatrixXf(3, 3)) );
CALL_SUBTEST( cwiseops(MatrixXf(22, 22)) );
CALL_SUBTEST( cwiseops(MatrixXi(8, 12)) );
CALL_SUBTEST( cwiseops(MatrixXd(20, 20)) );
CALL_SUBTEST_1( cwiseops(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( cwiseops(Matrix4d()) );
CALL_SUBTEST_3( cwiseops(MatrixXf(3, 3)) );
CALL_SUBTEST_4( cwiseops(MatrixXf(22, 22)) );
CALL_SUBTEST_5( cwiseops(MatrixXi(8, 12)) );
CALL_SUBTEST_6( cwiseops(MatrixXd(20, 20)) );
}
}

View File

@ -65,12 +65,12 @@ template<typename MatrixType> void determinant(const MatrixType& m)
void test_determinant()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( determinant(Matrix<float, 1, 1>()) );
CALL_SUBTEST( determinant(Matrix<double, 2, 2>()) );
CALL_SUBTEST( determinant(Matrix<double, 3, 3>()) );
CALL_SUBTEST( determinant(Matrix<double, 4, 4>()) );
CALL_SUBTEST( determinant(Matrix<std::complex<double>, 10, 10>()) );
CALL_SUBTEST( determinant(MatrixXd(20, 20)) );
CALL_SUBTEST_1( determinant(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( determinant(Matrix<double, 2, 2>()) );
CALL_SUBTEST_3( determinant(Matrix<double, 3, 3>()) );
CALL_SUBTEST_4( determinant(Matrix<double, 4, 4>()) );
CALL_SUBTEST_5( determinant(Matrix<std::complex<double>, 10, 10>()) );
CALL_SUBTEST_6( determinant(MatrixXd(20, 20)) );
}
CALL_SUBTEST( determinant(MatrixXd(200, 200)) );
CALL_SUBTEST_6( determinant(MatrixXd(200, 200)) );
}

View File

@ -95,14 +95,14 @@ template<typename MatrixType> void diagonalmatrices(const MatrixType& m)
void test_diagonalmatrices()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( diagonalmatrices(Matrix<float, 1, 1>()) );
CALL_SUBTEST( diagonalmatrices(Matrix3f()) );
CALL_SUBTEST( diagonalmatrices(Matrix<double,3,3,RowMajor>()) );
CALL_SUBTEST( diagonalmatrices(Matrix4d()) );
CALL_SUBTEST( diagonalmatrices(Matrix<float,4,4,RowMajor>()) );
CALL_SUBTEST( diagonalmatrices(MatrixXcf(3, 5)) );
CALL_SUBTEST( diagonalmatrices(MatrixXi(10, 8)) );
CALL_SUBTEST( diagonalmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
CALL_SUBTEST( diagonalmatrices(MatrixXf(21, 24)) );
CALL_SUBTEST_1( diagonalmatrices(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( diagonalmatrices(Matrix3f()) );
CALL_SUBTEST_3( diagonalmatrices(Matrix<double,3,3,RowMajor>()) );
CALL_SUBTEST_4( diagonalmatrices(Matrix4d()) );
CALL_SUBTEST_5( diagonalmatrices(Matrix<float,4,4,RowMajor>()) );
CALL_SUBTEST_6( diagonalmatrices(MatrixXcf(3, 5)) );
CALL_SUBTEST_7( diagonalmatrices(MatrixXi(10, 8)) );
CALL_SUBTEST_8( diagonalmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
CALL_SUBTEST_9( diagonalmatrices(MatrixXf(21, 24)) );
}
}

View File

@ -112,11 +112,11 @@ void test_dynalloc()
for (int i=0; i<g_repeat*100; ++i)
{
CALL_SUBTEST( check_dynaligned<Vector4f>() );
CALL_SUBTEST( check_dynaligned<Vector2d>() );
CALL_SUBTEST( check_dynaligned<Matrix4f>() );
CALL_SUBTEST( check_dynaligned<Vector4d>() );
CALL_SUBTEST( check_dynaligned<Vector4i>() );
CALL_SUBTEST(check_dynaligned<Vector4f>() );
CALL_SUBTEST(check_dynaligned<Vector2d>() );
CALL_SUBTEST(check_dynaligned<Matrix4f>() );
CALL_SUBTEST(check_dynaligned<Vector4d>() );
CALL_SUBTEST(check_dynaligned<Vector4i>() );
}
// check static allocation, who knows ?

View File

@ -58,7 +58,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
void test_eigensolver_complex()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( eigensolver(Matrix4cf()) );
CALL_SUBTEST( eigensolver(MatrixXcd(14,14)) );
CALL_SUBTEST_1( eigensolver(Matrix4cf()) );
CALL_SUBTEST_2( eigensolver(MatrixXcd(14,14)) );
}
}

View File

@ -75,18 +75,18 @@ template<typename MatrixType> void eigensolver_verify_assert()
void test_eigensolver_generic()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( eigensolver(Matrix4f()) );
CALL_SUBTEST( eigensolver(MatrixXd(17,17)) );
CALL_SUBTEST_1( eigensolver(Matrix4f()) );
CALL_SUBTEST_2( eigensolver(MatrixXd(17,17)) );
// some trivial but implementation-wise tricky cases
CALL_SUBTEST( eigensolver(MatrixXd(1,1)) );
CALL_SUBTEST( eigensolver(MatrixXd(2,2)) );
CALL_SUBTEST( eigensolver(Matrix<double,1,1>()) );
CALL_SUBTEST( eigensolver(Matrix<double,2,2>()) );
CALL_SUBTEST_2( eigensolver(MatrixXd(1,1)) );
CALL_SUBTEST_2( eigensolver(MatrixXd(2,2)) );
CALL_SUBTEST_3( eigensolver(Matrix<double,1,1>()) );
CALL_SUBTEST_4( eigensolver(Matrix2d()) );
}
CALL_SUBTEST( eigensolver_verify_assert<Matrix3f>() );
CALL_SUBTEST( eigensolver_verify_assert<Matrix3d>() );
CALL_SUBTEST( eigensolver_verify_assert<MatrixXf>() );
CALL_SUBTEST( eigensolver_verify_assert<MatrixXd>() );
CALL_SUBTEST_1( eigensolver_verify_assert<Matrix4f>() );
CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXd>() );
CALL_SUBTEST_4( eigensolver_verify_assert<Matrix2d>() );
CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
}

View File

@ -117,17 +117,17 @@ void test_eigensolver_selfadjoint()
{
for(int i = 0; i < g_repeat; i++) {
// very important to test a 3x3 matrix since we provide a special path for it
CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(10,10)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(17,17)) );
CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(10,10)) );
CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(19,19)) );
CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(17,17)) );
// some trivial but implementation-wise tricky cases
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(1,1)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(2,2)) );
CALL_SUBTEST( selfadjointeigensolver(Matrix<double,1,1>()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix<double,2,2>()) );
CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
}
}

View File

@ -75,8 +75,8 @@ template<typename BoxType> void alignedbox(const BoxType& _box)
void test_geo_alignedbox()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( alignedbox(AlignedBox<float,2>()) );
CALL_SUBTEST( alignedbox(AlignedBox<float,3>()) );
CALL_SUBTEST( alignedbox(AlignedBox<double,4>()) );
CALL_SUBTEST_1( alignedbox(AlignedBox<float,2>()) );
CALL_SUBTEST_2( alignedbox(AlignedBox<float,3>()) );
CALL_SUBTEST_3( alignedbox(AlignedBox<double,4>()) );
}
}

View File

@ -64,7 +64,7 @@ template<typename Scalar> void eulerangles(void)
void test_geo_eulerangles()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( eulerangles<float>() );
CALL_SUBTEST( eulerangles<double>() );
CALL_SUBTEST_1( eulerangles<float>() );
CALL_SUBTEST_2( eulerangles<double>() );
}
}

View File

@ -104,8 +104,8 @@ template<typename Scalar,int Size> void homogeneous(void)
void test_geo_homogeneous()
{
for(int i = 0; i < g_repeat; i++) {
// CALL_SUBTEST(( homogeneous<float,1>() ));
CALL_SUBTEST(( homogeneous<double,3>() ));
// CALL_SUBTEST(( homogeneous<double,8>() ));
CALL_SUBTEST_1(( homogeneous<float,1>() ));
CALL_SUBTEST_2(( homogeneous<double,3>() ));
CALL_SUBTEST_3(( homogeneous<double,8>() ));
}
}

View File

@ -131,11 +131,11 @@ template<typename Scalar> void lines()
void test_geo_hyperplane()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( hyperplane(Hyperplane<float,2>()) );
CALL_SUBTEST( hyperplane(Hyperplane<float,3>()) );
CALL_SUBTEST( hyperplane(Hyperplane<double,4>()) );
CALL_SUBTEST( hyperplane(Hyperplane<std::complex<double>,5>()) );
CALL_SUBTEST( lines<float>() );
CALL_SUBTEST( lines<double>() );
CALL_SUBTEST_1( hyperplane(Hyperplane<float,2>()) );
CALL_SUBTEST_2( hyperplane(Hyperplane<float,3>()) );
CALL_SUBTEST_3( hyperplane(Hyperplane<double,4>()) );
CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) );
CALL_SUBTEST_1( lines<float>() );
CALL_SUBTEST_2( lines<double>() );
}
}

View File

@ -113,15 +113,15 @@ template<typename Scalar, int Size> void orthomethods(int size=Size)
void test_geo_orthomethods()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( orthomethods_3<float>() );
CALL_SUBTEST( orthomethods_3<double>() );
CALL_SUBTEST( (orthomethods<float,2>()) );
CALL_SUBTEST( (orthomethods<double,2>()) );
CALL_SUBTEST( (orthomethods<float,3>()) );
CALL_SUBTEST( (orthomethods<double,3>()) );
CALL_SUBTEST( (orthomethods<float,7>()) );
CALL_SUBTEST( (orthomethods<std::complex<double>,8>()) );
CALL_SUBTEST( (orthomethods<float,Dynamic>(36)) );
CALL_SUBTEST( (orthomethods<double,Dynamic>(35)) );
CALL_SUBTEST_1( orthomethods_3<float>() );
CALL_SUBTEST_2( orthomethods_3<double>() );
CALL_SUBTEST_1( (orthomethods<float,2>()) );
CALL_SUBTEST_2( (orthomethods<double,2>()) );
CALL_SUBTEST_1( (orthomethods<float,3>()) );
CALL_SUBTEST_2( (orthomethods<double,3>()) );
CALL_SUBTEST_3( (orthomethods<float,7>()) );
CALL_SUBTEST_4( (orthomethods<std::complex<double>,8>()) );
CALL_SUBTEST_5( (orthomethods<float,Dynamic>(36)) );
CALL_SUBTEST_6( (orthomethods<double,Dynamic>(35)) );
}
}

View File

@ -69,9 +69,9 @@ template<typename LineType> void parametrizedline(const LineType& _line)
void test_geo_parametrizedline()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( parametrizedline(ParametrizedLine<float,2>()) );
CALL_SUBTEST( parametrizedline(ParametrizedLine<float,3>()) );
CALL_SUBTEST( parametrizedline(ParametrizedLine<double,4>()) );
CALL_SUBTEST( parametrizedline(ParametrizedLine<std::complex<double>,5>()) );
CALL_SUBTEST_1( parametrizedline(ParametrizedLine<float,2>()) );
CALL_SUBTEST_2( parametrizedline(ParametrizedLine<float,3>()) );
CALL_SUBTEST_3( parametrizedline(ParametrizedLine<double,4>()) );
CALL_SUBTEST_4( parametrizedline(ParametrizedLine<std::complex<double>,5>()) );
}
}

View File

@ -92,9 +92,12 @@ template<typename Scalar> void quaternion(void)
VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
VERIFY_IS_APPROX( v1.normalized(),(q2.setFromTwoVectors(v1, v1)*v1).normalized());
VERIFY_IS_APPROX(-v1.normalized(),(q2.setFromTwoVectors(v1,-v1)*v1).normalized());
v3 = v1.cwise()+eps;
VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
if (ei_is_same_type<Scalar,double>::ret)
{
v3 = v1.cwise()+eps;
VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
}
// inverse and conjugate
VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
@ -110,7 +113,7 @@ template<typename Scalar> void quaternion(void)
void test_geo_quaternion()
{
for(int i = 0; i < g_repeat; i++) {
// CALL_SUBTEST( quaternion<float>() );
CALL_SUBTEST( quaternion<double>() );
CALL_SUBTEST_1( quaternion<float>() );
CALL_SUBTEST_2( quaternion<double>() );
}
}

View File

@ -298,16 +298,19 @@ template<typename Scalar, int Mode> void transformations(void)
VERIFY_IS_APPROX((t0 * v1).template start<3>(), AlignedScaling3(v0) * v1);
// test transform inversion
if(Mode!=AffineCompact)
{
t0.setIdentity();
t0.translate(v0);
t0.linear().setRandom();
VERIFY_IS_APPROX(t0.inverse(Affine).matrix(), t0.matrix().inverse());
t0.setIdentity();
t0.translate(v0).rotate(q1);
VERIFY_IS_APPROX(t0.inverse(Isometry).matrix(), t0.matrix().inverse());
}
t0.setIdentity();
t0.translate(v0);
t0.linear().setRandom();
Matrix4 t044 = Matrix4::Zero();
t044(3,3) = 1;
t044.block(0,0,t0.matrix().rows(),4) = t0.matrix();
VERIFY_IS_APPROX(t0.inverse(Affine).matrix(), t044.inverse().block(0,0,t0.matrix().rows(),4));
t0.setIdentity();
t0.translate(v0).rotate(q1);
t044 = Matrix4::Zero();
t044(3,3) = 1;
t044.block(0,0,t0.matrix().rows(),4) = t0.matrix();
VERIFY_IS_APPROX(t0.inverse(Isometry).matrix(), t044.inverse().block(0,0,t0.matrix().rows(),4));
// test extract rotation and aligned scaling
// t0.setIdentity();
@ -354,9 +357,8 @@ template<typename Scalar, int Mode> void transformations(void)
void test_geo_transformations()
{
for(int i = 0; i < g_repeat; i++) {
// CALL_SUBTEST( transformations<float>() );
CALL_SUBTEST(( transformations<double,Affine>() ));
CALL_SUBTEST(( transformations<double,AffineCompact>() ));
CALL_SUBTEST(( transformations<double,Projective>() ));
CALL_SUBTEST_1(( transformations<double,Affine>() ));
CALL_SUBTEST_2(( transformations<float,AffineCompact>() ));
CALL_SUBTEST_3(( transformations<double,Projective>() ));
}
}

View File

@ -92,12 +92,12 @@ template<typename MatrixType> void householder(const MatrixType& m)
void test_householder()
{
for(int i = 0; i < 2*g_repeat; i++) {
CALL_SUBTEST( householder(Matrix<double,2,2>()) );
CALL_SUBTEST( householder(Matrix<float,2,3>()) );
CALL_SUBTEST( householder(Matrix<double,3,5>()) );
CALL_SUBTEST( householder(Matrix<float,4,4>()) );
CALL_SUBTEST( householder(MatrixXd(10,12)) );
CALL_SUBTEST( householder(MatrixXcf(16,17)) );
CALL_SUBTEST_1( householder(Matrix<double,2,2>()) );
CALL_SUBTEST_2( householder(Matrix<float,2,3>()) );
CALL_SUBTEST_3( householder(Matrix<double,3,5>()) );
CALL_SUBTEST_4( householder(Matrix<float,4,4>()) );
CALL_SUBTEST_5( householder(MatrixXd(10,12)) );
CALL_SUBTEST_6( householder(MatrixXcf(16,17)) );
}
}

View File

@ -53,9 +53,6 @@ template<typename MatrixType> void inverse(const MatrixType& m)
m2 = m1.inverse();
VERIFY_IS_APPROX(m1, m2.inverse() );
m1.computeInverse(&m2);
VERIFY_IS_APPROX(m1, m2.inverse() );
VERIFY_IS_APPROX((Scalar(2)*m2).inverse(), m2.inverse()*Scalar(0.5));
VERIFY_IS_APPROX(identity, m1.inverse() * m1 );
@ -66,36 +63,53 @@ template<typename MatrixType> void inverse(const MatrixType& m)
// since for the general case we implement separately row-major and col-major, test that
VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose());
//computeInverseWithCheck tests
#if !defined(EIGEN_TEST_PART_5) && !defined(EIGEN_TEST_PART_6)
//computeInverseAndDetWithCheck tests
//First: an invertible matrix
bool invertible = m1.computeInverseWithCheck(&m2);
bool invertible;
RealScalar det;
m2.setZero();
m1.computeInverseAndDetWithCheck(m2, det, invertible);
VERIFY(invertible);
VERIFY_IS_APPROX(identity, m1*m2);
VERIFY_IS_APPROX(det, m1.determinant());
m2.setZero();
m1.computeInverseWithCheck(m2, invertible);
VERIFY(invertible);
VERIFY_IS_APPROX(identity, m1*m2);
//Second: a rank one matrix (not invertible, except for 1x1 matrices)
VectorType v3 = VectorType::Random(rows);
MatrixType m3 = v3*v3.transpose(), m4(rows,cols);
invertible = m3.computeInverseWithCheck( &m4 );
m3.computeInverseAndDetWithCheck(m4, det, invertible);
VERIFY( rows==1 ? invertible : !invertible );
VERIFY_IS_APPROX(det, m3.determinant());
m3.computeInverseWithCheck(m4, invertible);
VERIFY( rows==1 ? invertible : !invertible );
#endif
}
void test_inverse()
{
int s;
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( inverse(Matrix<double,1,1>()) );
CALL_SUBTEST( inverse(Matrix2d()) );
CALL_SUBTEST( inverse(Matrix3f()) );
CALL_SUBTEST( inverse(Matrix4f()) );
CALL_SUBTEST_1( inverse(Matrix<double,1,1>()) );
CALL_SUBTEST_2( inverse(Matrix2d()) );
CALL_SUBTEST_3( inverse(Matrix3f()) );
CALL_SUBTEST_4( inverse(Matrix4f()) );
s = ei_random<int>(50,320);
CALL_SUBTEST( inverse(MatrixXf(s,s)) );
CALL_SUBTEST_5( inverse(MatrixXf(s,s)) );
s = ei_random<int>(25,100);
CALL_SUBTEST( inverse(MatrixXcd(s,s)) );
CALL_SUBTEST_6( inverse(MatrixXcd(s,s)) );
}
#ifdef EIGEN_TEST_PART_4
// test some tricky cases for 4x4 matrices
VERIFY_IS_APPROX((Matrix4f() << 0,0,1,0, 1,0,0,0, 0,1,0,0, 0,0,0,1).finished().inverse(),
(Matrix4f() << 0,1,0,0, 0,0,1,0, 1,0,0,0, 0,0,0,1).finished());
VERIFY_IS_APPROX((Matrix4f() << 1,0,0,0, 0,0,1,0, 0,0,0,1, 0,1,0,0).finished().inverse(),
(Matrix4f() << 1,0,0,0, 0,0,0,1, 0,1,0,0, 0,0,1,0).finished());
#endif
}

View File

@ -80,27 +80,27 @@ void test_jacobisvd()
Matrix2cd m;
m << 0, 1,
0, 1;
CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) ));
CALL_SUBTEST_1(( svd<Matrix2cd,0>(m, false) ));
m << 1, 0,
1, 0;
CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) ));
CALL_SUBTEST_1(( svd<Matrix2cd,0>(m, false) ));
Matrix2d n;
n << 1, 1,
1, -1;
CALL_SUBTEST(( svd<Matrix2d,0>(n, false) ));
CALL_SUBTEST(( svd<Matrix3f,0>() ));
CALL_SUBTEST(( svd<Matrix4d,Square>() ));
CALL_SUBTEST(( svd<Matrix<float,3,5> , AtLeastAsManyColsAsRows>() ));
CALL_SUBTEST(( svd<Matrix<double,Dynamic,2> , AtLeastAsManyRowsAsCols>(Matrix<double,Dynamic,2>(10,2)) ));
CALL_SUBTEST_2(( svd<Matrix2d,0>(n, false) ));
CALL_SUBTEST_3(( svd<Matrix3f,0>() ));
CALL_SUBTEST_4(( svd<Matrix4d,Square>() ));
CALL_SUBTEST_5(( svd<Matrix<float,3,5> , AtLeastAsManyColsAsRows>() ));
CALL_SUBTEST_6(( svd<Matrix<double,Dynamic,2> , AtLeastAsManyRowsAsCols>(Matrix<double,Dynamic,2>(10,2)) ));
CALL_SUBTEST(( svd<MatrixXf,Square>(MatrixXf(50,50)) ));
CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyRowsAsCols>(MatrixXcd(14,7)) ));
CALL_SUBTEST_7(( svd<MatrixXf,Square>(MatrixXf(50,50)) ));
CALL_SUBTEST_8(( svd<MatrixXcd,AtLeastAsManyRowsAsCols>(MatrixXcd(14,7)) ));
}
CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) ));
CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) ));
CALL_SUBTEST_9(( svd<MatrixXf,0>(MatrixXf(300,200)) ));
CALL_SUBTEST_10(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) ));
CALL_SUBTEST(( svd_verify_assert<Matrix3f>() ));
CALL_SUBTEST(( svd_verify_assert<Matrix3d>() ));
CALL_SUBTEST(( svd_verify_assert<MatrixXf>() ));
CALL_SUBTEST(( svd_verify_assert<MatrixXd>() ));
CALL_SUBTEST_3(( svd_verify_assert<Matrix3f>() ));
CALL_SUBTEST_3(( svd_verify_assert<Matrix3d>() ));
CALL_SUBTEST_9(( svd_verify_assert<MatrixXf>() ));
CALL_SUBTEST_11(( svd_verify_assert<MatrixXd>() ));
}

View File

@ -87,13 +87,13 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
void test_linearstructure()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( linearStructure(Matrix<float, 1, 1>()) );
CALL_SUBTEST( linearStructure(Matrix2f()) );
CALL_SUBTEST( linearStructure(Vector3d()) );
CALL_SUBTEST( linearStructure(Matrix4d()) );
CALL_SUBTEST( linearStructure(MatrixXcf(3, 3)) );
CALL_SUBTEST( linearStructure(MatrixXf(8, 12)) );
CALL_SUBTEST( linearStructure(MatrixXi(8, 12)) );
CALL_SUBTEST( linearStructure(MatrixXcd(20, 20)) );
CALL_SUBTEST_1( linearStructure(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( linearStructure(Matrix2f()) );
CALL_SUBTEST_3( linearStructure(Vector3d()) );
CALL_SUBTEST_4( linearStructure(Matrix4d()) );
CALL_SUBTEST_5( linearStructure(MatrixXcf(3, 3)) );
CALL_SUBTEST_6( linearStructure(MatrixXf(8, 12)) );
CALL_SUBTEST_7( linearStructure(MatrixXi(8, 12)) );
CALL_SUBTEST_8( linearStructure(MatrixXcd(20, 20)) );
}
}

View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-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
@ -30,15 +30,44 @@ template<typename MatrixType> void lu_non_invertible()
/* this test covers the following files:
LU.h
*/
int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200);
int rows, cols, cols2;
if(MatrixType::RowsAtCompileTime==Dynamic)
{
rows = ei_random<int>(20,200);
}
else
{
rows = MatrixType::RowsAtCompileTime;
}
if(MatrixType::ColsAtCompileTime==Dynamic)
{
cols = ei_random<int>(20,200);
cols2 = ei_random<int>(20,200);
}
else
{
cols2 = cols = MatrixType::ColsAtCompileTime;
}
typedef typename ei_kernel_retval_base<FullPivLU<MatrixType> >::ReturnMatrixType KernelMatrixType;
typedef typename ei_image_retval_base<FullPivLU<MatrixType> >::ReturnMatrixType ImageMatrixType;
typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType;
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime>
CMatrixType;
int rank = ei_random<int>(1, std::min(rows, cols)-1);
MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
// The image of the zero matrix should consist of a single (zero) column vector
VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1));
MatrixType m1(rows, cols), m3(rows, cols2);
CMatrixType m2(cols, cols2);
createRandomMatrixOfRank(rank, rows, cols, m1);
LU<MatrixType> lu(m1);
typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel();
typename LU<MatrixType>::ImageResultType m1image = lu.image();
FullPivLU<MatrixType> lu(m1);
std::cout << lu.kernel().rows() << " " << lu.kernel().cols() << std::endl;
KernelMatrixType m1kernel = lu.kernel();
ImageMatrixType m1image = lu.image(m1);
VERIFY(rank == lu.rank());
VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
@ -46,22 +75,16 @@ template<typename MatrixType> void lu_non_invertible()
VERIFY(!lu.isInvertible());
VERIFY(!lu.isSurjective());
VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
VERIFY(m1image.lu().rank() == rank);
MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
VERIFY(m1image.fullPivLu().rank() == rank);
DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
sidebyside << m1, m1image;
VERIFY(sidebyside.lu().rank() == rank);
m2 = MatrixType::Random(cols,cols2);
VERIFY(sidebyside.fullPivLu().rank() == rank);
m2 = CMatrixType::Random(cols,cols2);
m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
VERIFY(lu.solve(m3, &m2));
m2 = CMatrixType::Random(cols,cols2);
// test that the code, which does resize(), may be applied to an xpr
m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
m3 = MatrixType::Random(rows,cols2);
VERIFY(!lu.solve(m3, &m2));
typedef Matrix<typename MatrixType::Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
SquareMatrixType m4(rows, rows), m5(rows, rows);
createRandomMatrixOfRank(rows/2, rows, rows, m4);
VERIFY(!m4.computeInverseWithCheck(&m5));
}
template<typename MatrixType> void lu_invertible()
@ -82,69 +105,70 @@ template<typename MatrixType> void lu_invertible()
m1 += a * a.adjoint();
}
LU<MatrixType> lu(m1);
FullPivLU<MatrixType> lu(m1);
VERIFY(0 == lu.dimensionOfKernel());
VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
VERIFY(size == lu.rank());
VERIFY(lu.isInjective());
VERIFY(lu.isSurjective());
VERIFY(lu.isInvertible());
VERIFY(lu.image().lu().isInvertible());
VERIFY(lu.image(m1).fullPivLu().isInvertible());
m3 = MatrixType::Random(size,size);
lu.solve(m3, &m2);
m2 = lu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
VERIFY_IS_APPROX(m2, lu.inverse()*m3);
m3 = MatrixType::Random(size,size);
VERIFY(lu.solve(m3, &m2));
}
template<typename MatrixType> void lu_verify_assert()
{
MatrixType tmp;
LU<MatrixType> lu;
FullPivLU<MatrixType> lu;
VERIFY_RAISES_ASSERT(lu.matrixLU())
VERIFY_RAISES_ASSERT(lu.permutationP())
VERIFY_RAISES_ASSERT(lu.permutationQ())
VERIFY_RAISES_ASSERT(lu.computeKernel(&tmp))
VERIFY_RAISES_ASSERT(lu.computeImage(&tmp))
VERIFY_RAISES_ASSERT(lu.kernel())
VERIFY_RAISES_ASSERT(lu.image())
VERIFY_RAISES_ASSERT(lu.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(lu.image(tmp))
VERIFY_RAISES_ASSERT(lu.solve(tmp))
VERIFY_RAISES_ASSERT(lu.determinant())
VERIFY_RAISES_ASSERT(lu.rank())
VERIFY_RAISES_ASSERT(lu.dimensionOfKernel())
VERIFY_RAISES_ASSERT(lu.isInjective())
VERIFY_RAISES_ASSERT(lu.isSurjective())
VERIFY_RAISES_ASSERT(lu.isInvertible())
VERIFY_RAISES_ASSERT(lu.computeInverse(&tmp))
VERIFY_RAISES_ASSERT(lu.inverse())
PartialLU<MatrixType> plu;
PartialPivLU<MatrixType> plu;
VERIFY_RAISES_ASSERT(plu.matrixLU())
VERIFY_RAISES_ASSERT(plu.permutationP())
VERIFY_RAISES_ASSERT(plu.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(plu.solve(tmp))
VERIFY_RAISES_ASSERT(plu.determinant())
VERIFY_RAISES_ASSERT(plu.computeInverse(&tmp))
VERIFY_RAISES_ASSERT(plu.inverse())
}
void test_lu()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( lu_non_invertible<MatrixXf>() );
CALL_SUBTEST( lu_non_invertible<MatrixXd>() );
CALL_SUBTEST( lu_non_invertible<MatrixXcf>() );
CALL_SUBTEST( lu_non_invertible<MatrixXcd>() );
CALL_SUBTEST( lu_invertible<MatrixXf>() );
CALL_SUBTEST( lu_invertible<MatrixXd>() );
CALL_SUBTEST( lu_invertible<MatrixXcf>() );
CALL_SUBTEST( lu_invertible<MatrixXcd>() );
}
CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() );
CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() );
CALL_SUBTEST( lu_verify_assert<Matrix3f>() );
CALL_SUBTEST( lu_verify_assert<Matrix3d>() );
CALL_SUBTEST( lu_verify_assert<MatrixXf>() );
CALL_SUBTEST( lu_verify_assert<MatrixXd>() );
CALL_SUBTEST( lu_verify_assert<MatrixXcf>() );
CALL_SUBTEST( lu_verify_assert<MatrixXcd>() );
CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) );
CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) );
CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() );
CALL_SUBTEST_3( lu_invertible<MatrixXf>() );
CALL_SUBTEST_3( lu_verify_assert<MatrixXf>() );
CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() );
CALL_SUBTEST_4( lu_invertible<MatrixXd>() );
CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() );
CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() );
CALL_SUBTEST_5( lu_invertible<MatrixXcf>() );
CALL_SUBTEST_5( lu_verify_assert<MatrixXcf>() );
CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() );
CALL_SUBTEST_6( lu_invertible<MatrixXcd>() );
CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
}
}

View File

@ -170,6 +170,102 @@ namespace Eigen
g_test_stack.pop_back(); \
} while (0)
#ifdef EIGEN_TEST_PART_1
#define CALL_SUBTEST_1(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_1(FUNC)
#endif
#ifdef EIGEN_TEST_PART_2
#define CALL_SUBTEST_2(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_2(FUNC)
#endif
#ifdef EIGEN_TEST_PART_3
#define CALL_SUBTEST_3(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_3(FUNC)
#endif
#ifdef EIGEN_TEST_PART_4
#define CALL_SUBTEST_4(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_4(FUNC)
#endif
#ifdef EIGEN_TEST_PART_5
#define CALL_SUBTEST_5(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_5(FUNC)
#endif
#ifdef EIGEN_TEST_PART_6
#define CALL_SUBTEST_6(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_6(FUNC)
#endif
#ifdef EIGEN_TEST_PART_7
#define CALL_SUBTEST_7(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_7(FUNC)
#endif
#ifdef EIGEN_TEST_PART_8
#define CALL_SUBTEST_8(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_8(FUNC)
#endif
#ifdef EIGEN_TEST_PART_9
#define CALL_SUBTEST_9(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_9(FUNC)
#endif
#ifdef EIGEN_TEST_PART_10
#define CALL_SUBTEST_10(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_10(FUNC)
#endif
#ifdef EIGEN_TEST_PART_11
#define CALL_SUBTEST_11(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_11(FUNC)
#endif
#ifdef EIGEN_TEST_PART_12
#define CALL_SUBTEST_12(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_12(FUNC)
#endif
#ifdef EIGEN_TEST_PART_13
#define CALL_SUBTEST_13(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_13(FUNC)
#endif
#ifdef EIGEN_TEST_PART_14
#define CALL_SUBTEST_14(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_14(FUNC)
#endif
#ifdef EIGEN_TEST_PART_15
#define CALL_SUBTEST_15(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_15(FUNC)
#endif
#ifdef EIGEN_TEST_PART_16
#define CALL_SUBTEST_16(FUNC) CALL_SUBTEST(FUNC)
#else
#define CALL_SUBTEST_16(FUNC)
#endif
namespace Eigen {
template<typename T> inline typename NumTraits<T>::Real test_precision();

View File

@ -81,16 +81,16 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
void test_map()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( map_class(Matrix<float, 1, 1>()) );
CALL_SUBTEST( map_class(Vector4d()) );
CALL_SUBTEST( map_class(RowVector4f()) );
CALL_SUBTEST( map_class(VectorXcf(8)) );
CALL_SUBTEST( map_class(VectorXi(12)) );
CALL_SUBTEST_1( map_class(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( map_class(Vector4d()) );
CALL_SUBTEST_3( map_class(RowVector4f()) );
CALL_SUBTEST_4( map_class(VectorXcf(8)) );
CALL_SUBTEST_5( map_class(VectorXi(12)) );
CALL_SUBTEST( map_static_methods(Matrix<double, 1, 1>()) );
CALL_SUBTEST( map_static_methods(Vector3f()) );
CALL_SUBTEST( map_static_methods(RowVector3d()) );
CALL_SUBTEST( map_static_methods(VectorXcd(8)) );
CALL_SUBTEST( map_static_methods(VectorXf(12)) );
CALL_SUBTEST_6( map_static_methods(Matrix<double, 1, 1>()) );
CALL_SUBTEST_7( map_static_methods(Vector3f()) );
CALL_SUBTEST_8( map_static_methods(RowVector3d()) );
CALL_SUBTEST_9( map_static_methods(VectorXcd(8)) );
CALL_SUBTEST_10( map_static_methods(VectorXf(12)) );
}
}

View File

@ -54,10 +54,10 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m)
void test_miscmatrices()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( miscMatrices(Matrix<float, 1, 1>()) );
CALL_SUBTEST( miscMatrices(Matrix4d()) );
CALL_SUBTEST( miscMatrices(MatrixXcf(3, 3)) );
CALL_SUBTEST( miscMatrices(MatrixXi(8, 12)) );
CALL_SUBTEST( miscMatrices(MatrixXcd(20, 20)) );
CALL_SUBTEST_1( miscMatrices(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( miscMatrices(Matrix4d()) );
CALL_SUBTEST_3( miscMatrices(MatrixXcf(3, 3)) );
CALL_SUBTEST_4( miscMatrices(MatrixXi(8, 12)) );
CALL_SUBTEST_5( miscMatrices(MatrixXcd(20, 20)) );
}
}

View File

@ -174,10 +174,10 @@ template<int SizeAtCompileType> void mixingtypes_small()
void test_mixingtypes()
{
// check that our operator new is indeed called:
CALL_SUBTEST(mixingtypes<3>());
CALL_SUBTEST(mixingtypes<4>());
CALL_SUBTEST(mixingtypes<Dynamic>(20));
CALL_SUBTEST_1(mixingtypes<3>());
CALL_SUBTEST_2(mixingtypes<4>());
CALL_SUBTEST_3(mixingtypes<Dynamic>(20));
CALL_SUBTEST(mixingtypes_small<4>());
CALL_SUBTEST(mixingtypes_large(20));
CALL_SUBTEST_4(mixingtypes_small<4>());
CALL_SUBTEST_5(mixingtypes_large(20));
}

View File

@ -77,7 +77,7 @@ void test_nomalloc()
{
// check that our operator new is indeed called:
VERIFY_RAISES_ASSERT(MatrixXd dummy = MatrixXd::Random(3,3));
CALL_SUBTEST( nomalloc(Matrix<float, 1, 1>()) );
CALL_SUBTEST( nomalloc(Matrix4d()) );
CALL_SUBTEST( nomalloc(Matrix<float,32,32>()) );
CALL_SUBTEST(nomalloc(Matrix<float, 1, 1>()) );
CALL_SUBTEST(nomalloc(Matrix4d()) );
CALL_SUBTEST(nomalloc(Matrix<float,32,32>()) );
}

View File

@ -233,12 +233,12 @@ template<typename Scalar> void packetmath_real()
void test_packetmath()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( packetmath<float>() );
CALL_SUBTEST( packetmath<double>() );
CALL_SUBTEST( packetmath<int>() );
CALL_SUBTEST( packetmath<std::complex<float> >() );
CALL_SUBTEST_1( packetmath<float>() );
CALL_SUBTEST_2( packetmath<double>() );
CALL_SUBTEST_3( packetmath<int>() );
CALL_SUBTEST_1( packetmath<std::complex<float> >() );
CALL_SUBTEST( packetmath_real<float>() );
CALL_SUBTEST( packetmath_real<double>() );
CALL_SUBTEST_1( packetmath_real<float>() );
CALL_SUBTEST_2( packetmath_real<double>() );
}
}

View File

@ -119,8 +119,8 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
void test_product_extra()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) );
CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) );
CALL_SUBTEST_1( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) );
CALL_SUBTEST_2( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
CALL_SUBTEST_3( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) );
}
}

View File

@ -27,13 +27,14 @@
void test_product_large()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( product(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST( product(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST( product(MatrixXi(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST( product(MatrixXcf(ei_random<int>(1,50), ei_random<int>(1,50))) );
CALL_SUBTEST( product(Matrix<float,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST_1( product(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST_2( product(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST_3( product(MatrixXi(ei_random<int>(1,320), ei_random<int>(1,320))) );
CALL_SUBTEST_4( product(MatrixXcf(ei_random<int>(1,50), ei_random<int>(1,50))) );
CALL_SUBTEST_5( product(Matrix<float,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,320), ei_random<int>(1,320))) );
}
#if defined EIGEN_TEST_PART_6
{
// test a specific issue in DiagonalProduct
int N = 1000000;
@ -48,4 +49,5 @@ void test_product_large()
MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
VERIFY_IS_APPROX((a = a * b), (c * b).eval());
}
#endif
}

View File

@ -117,8 +117,8 @@ void test_product_notemporary()
int s;
for(int i = 0; i < g_repeat; i++) {
s = ei_random<int>(16,320);
CALL_SUBTEST( product_notemporary(MatrixXf(s, s)) );
CALL_SUBTEST_1( product_notemporary(MatrixXf(s, s)) );
s = ei_random<int>(16,120);
CALL_SUBTEST( product_notemporary(MatrixXcd(s,s)) );
CALL_SUBTEST_2( product_notemporary(MatrixXcd(s,s)) );
}
}

View File

@ -78,13 +78,13 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
void test_product_selfadjoint()
{
for(int i = 0; i < g_repeat ; i++) {
CALL_SUBTEST( product_selfadjoint(Matrix<float, 1, 1>()) );
CALL_SUBTEST( product_selfadjoint(Matrix<float, 2, 2>()) );
CALL_SUBTEST( product_selfadjoint(Matrix3d()) );
CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) );
CALL_SUBTEST( product_selfadjoint(MatrixXd(14,14)) );
CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
CALL_SUBTEST_1( product_selfadjoint(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( product_selfadjoint(Matrix<float, 2, 2>()) );
CALL_SUBTEST_3( product_selfadjoint(Matrix3d()) );
CALL_SUBTEST_4( product_selfadjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST_5( product_selfadjoint(MatrixXcd(21,21)) );
CALL_SUBTEST_6( product_selfadjoint(MatrixXd(14,14)) );
CALL_SUBTEST_7( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
CALL_SUBTEST_8( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
}
}

View File

@ -28,16 +28,18 @@
void test_product_small()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( product(Matrix<float, 3, 2>()) );
CALL_SUBTEST( product(Matrix<int, 3, 5>()) );
CALL_SUBTEST( product(Matrix3d()) );
CALL_SUBTEST( product(Matrix4d()) );
CALL_SUBTEST( product(Matrix4f()) );
CALL_SUBTEST_1( product(Matrix<float, 3, 2>()) );
CALL_SUBTEST_2( product(Matrix<int, 3, 5>()) );
CALL_SUBTEST_3( product(Matrix3d()) );
CALL_SUBTEST_4( product(Matrix4d()) );
CALL_SUBTEST_5( product(Matrix4f()) );
}
#ifdef EIGEN_TEST_PART_6
{
// test compilation of (outer_product) * vector
Vector3f v = Vector3f::Random();
VERIFY_IS_APPROX( (v * v.transpose()) * v, (v * v.transpose()).eval() * v);
}
#endif
}

View File

@ -108,10 +108,10 @@ void test_product_symm()
{
for(int i = 0; i < g_repeat ; i++)
{
CALL_SUBTEST(( symm<float,Dynamic,Dynamic>(ei_random<int>(10,320),ei_random<int>(10,320)) ));
CALL_SUBTEST(( symm<std::complex<double>,Dynamic,Dynamic>(ei_random<int>(10,320),ei_random<int>(10,320)) ));
CALL_SUBTEST_1(( symm<float,Dynamic,Dynamic>(ei_random<int>(10,320),ei_random<int>(10,320)) ));
CALL_SUBTEST_2(( symm<std::complex<double>,Dynamic,Dynamic>(ei_random<int>(10,320),ei_random<int>(10,320)) ));
CALL_SUBTEST(( symm<float,Dynamic,1>(ei_random<int>(10,320)) ));
CALL_SUBTEST(( symm<std::complex<double>,Dynamic,1>(ei_random<int>(10,320)) ));
CALL_SUBTEST_3(( symm<float,Dynamic,1>(ei_random<int>(10,320)) ));
CALL_SUBTEST_4(( symm<std::complex<double>,Dynamic,1>(ei_random<int>(10,320)) ));
}
}

View File

@ -75,8 +75,8 @@ void test_product_syrk()
{
int s;
s = ei_random<int>(10,320);
CALL_SUBTEST( syrk(MatrixXf(s, s)) );
CALL_SUBTEST_1( syrk(MatrixXf(s, s)) );
s = ei_random<int>(10,320);
CALL_SUBTEST( syrk(MatrixXcd(s, s)) );
CALL_SUBTEST_2( syrk(MatrixXcd(s, s)) );
}
}

View File

@ -63,7 +63,7 @@ void test_product_trmm()
{
for(int i = 0; i < g_repeat ; i++)
{
trmm<float>(ei_random<int>(1,320),ei_random<int>(1,320));
trmm<std::complex<double> >(ei_random<int>(1,320),ei_random<int>(1,320));
CALL_SUBTEST_1((trmm<float>(ei_random<int>(1,320),ei_random<int>(1,320))));
CALL_SUBTEST_2((trmm<std::complex<double> >(ei_random<int>(1,320),ei_random<int>(1,320))));
}
}

View File

@ -82,11 +82,11 @@ template<typename MatrixType> void trmv(const MatrixType& m)
void test_product_trmv()
{
for(int i = 0; i < g_repeat ; i++) {
CALL_SUBTEST( trmv(Matrix<float, 1, 1>()) );
CALL_SUBTEST( trmv(Matrix<float, 2, 2>()) );
CALL_SUBTEST( trmv(Matrix3d()) );
CALL_SUBTEST( trmv(Matrix<std::complex<float>,23, 23>()) );
CALL_SUBTEST( trmv(MatrixXcd(17,17)) );
CALL_SUBTEST( trmv(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) );
CALL_SUBTEST_1( trmv(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( trmv(Matrix<float, 2, 2>()) );
CALL_SUBTEST_3( trmv(Matrix3d()) );
CALL_SUBTEST_4( trmv(Matrix<std::complex<float>,23, 23>()) );
CALL_SUBTEST_5( trmv(MatrixXcd(17,17)) );
CALL_SUBTEST_6( trmv(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) );
}
}

View File

@ -59,7 +59,7 @@ void test_product_trsm()
{
for(int i = 0; i < g_repeat ; i++)
{
trsm<float>(ei_random<int>(1,320),ei_random<int>(1,320));
trsm<std::complex<double> >(ei_random<int>(1,320),ei_random<int>(1,320));
CALL_SUBTEST_1((trsm<float>(ei_random<int>(1,320),ei_random<int>(1,320))));
CALL_SUBTEST_2((trsm<std::complex<double> >(ei_random<int>(1,320),ei_random<int>(1,320))));
}
}

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())
@ -115,24 +115,24 @@ template<typename MatrixType> void qr_verify_assert()
void test_qr()
{
for(int i = 0; i < 1; i++) {
CALL_SUBTEST( qr(MatrixXf(47,40)) );
CALL_SUBTEST( qr(MatrixXcd(17,7)) );
CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));
CALL_SUBTEST_1( qr(MatrixXf(47,40)) );
CALL_SUBTEST_2( qr(MatrixXcd(17,7)) );
CALL_SUBTEST_3(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
CALL_SUBTEST_4(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));
}
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_1( qr_invertible<MatrixXf>() );
CALL_SUBTEST_6( qr_invertible<MatrixXd>() );
CALL_SUBTEST_7( qr_invertible<MatrixXcf>() );
CALL_SUBTEST_8( 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>());
CALL_SUBTEST_9(qr_verify_assert<Matrix3f>());
CALL_SUBTEST_10(qr_verify_assert<Matrix3d>());
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
CALL_SUBTEST_6(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_7(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_8(qr_verify_assert<MatrixXcd>());
}

View File

@ -36,7 +36,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
createRandomMatrixOfRank(rank,rows,cols,m1);
ColPivotingHouseholderQR<MatrixType> qr(m1);
ColPivHouseholderQR<MatrixType> qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
@ -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()
@ -74,7 +72,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1);
Matrix<Scalar,Rows,Cols> m1;
createRandomMatrixOfRank(rank,Rows,Cols,m1);
ColPivotingHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
@ -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()
@ -118,9 +114,9 @@ template<typename MatrixType> void qr_invertible()
m1 += a * a.adjoint();
}
ColPivotingHouseholderQR<MatrixType> qr(m1);
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
@ -138,15 +134,14 @@ template<typename MatrixType> void qr_verify_assert()
{
MatrixType tmp;
ColPivotingHouseholderQR<MatrixType> qr;
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())
@ -155,24 +150,24 @@ template<typename MatrixType> void qr_verify_assert()
void test_qr_colpivoting()
{
for(int i = 0; i < 1; i++) {
CALL_SUBTEST( qr<MatrixXf>() );
CALL_SUBTEST( qr<MatrixXd>() );
CALL_SUBTEST( qr<MatrixXcd>() );
CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
CALL_SUBTEST_1( qr<MatrixXf>() );
CALL_SUBTEST_2( qr<MatrixXd>() );
CALL_SUBTEST_3( qr<MatrixXcd>() );
CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
}
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_1( qr_invertible<MatrixXf>() );
CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
CALL_SUBTEST_6( qr_invertible<MatrixXcf>() );
CALL_SUBTEST_3( 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>());
CALL_SUBTEST_7(qr_verify_assert<Matrix3f>());
CALL_SUBTEST_8(qr_verify_assert<Matrix3d>());
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
}

Some files were not shown because too many files have changed in this diff Show More