mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-10 02:39:03 +08:00
introduce ei_xxx_return_value and ei_xxx_impl for xxx in solve,kernel,impl
put them in a new internal 'misc' directory
This commit is contained in:
Benoit Jacob
parent
f975b9bd3e
commit
da363d997f
@ -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
|
||||
|
||||
3
Eigen/LU
3
Eigen/LU
@ -18,6 +18,9 @@ namespace Eigen {
|
||||
* \endcode
|
||||
*/
|
||||
|
||||
#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"
|
||||
|
||||
@ -10,3 +10,4 @@ ADD_SUBDIRECTORY(Sparse)
|
||||
ADD_SUBDIRECTORY(Jacobi)
|
||||
ADD_SUBDIRECTORY(Householder)
|
||||
ADD_SUBDIRECTORY(Eigenvalues)
|
||||
ADD_SUBDIRECTORY(misc)
|
||||
@ -38,8 +38,10 @@ struct ei_traits<ReturnByValue<Derived> >
|
||||
// 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 better fix is probably rather to make Block work directly
|
||||
// on the nested type, right?
|
||||
// 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
|
||||
};
|
||||
|
||||
@ -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_return_value;
|
||||
template<typename DecompositionType, typename Rhs, typename Dest> struct ei_solve_impl;
|
||||
template<typename DecompositionType> struct ei_kernel_return_value;
|
||||
template<typename DecompositionType, typename Dest> struct ei_kernel_impl;
|
||||
template<typename DecompositionType> struct ei_image_return_value;
|
||||
template<typename DecompositionType, typename Dest> struct ei_image_impl;
|
||||
|
||||
template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int Supers=Dynamic, int Subs=Dynamic, int Options=0> class BandMatrix;
|
||||
|
||||
|
||||
|
||||
@ -25,10 +25,6 @@
|
||||
#ifndef EIGEN_LU_H
|
||||
#define EIGEN_LU_H
|
||||
|
||||
template<typename MatrixType, typename Rhs> struct ei_fullpivlu_solve_impl;
|
||||
template<typename MatrixType> struct ei_fullpivlu_kernel_impl;
|
||||
template<typename MatrixType> struct ei_fullpivlu_image_impl;
|
||||
|
||||
/** \ingroup LU_Module
|
||||
*
|
||||
* \class FullPivLU
|
||||
@ -59,10 +55,10 @@ template<typename MatrixType> struct ei_fullpivlu_image_impl;
|
||||
*
|
||||
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
|
||||
*/
|
||||
template<typename MatrixType> class FullPivLU
|
||||
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;
|
||||
@ -70,17 +66,12 @@ template<typename MatrixType> class FullPivLU
|
||||
typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
|
||||
|
||||
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
|
||||
MatrixType::MaxColsAtCompileTime,
|
||||
MatrixType::MaxRowsAtCompileTime)
|
||||
};
|
||||
|
||||
/**
|
||||
* \brief Default Constructor.
|
||||
*
|
||||
* The default constructor is useful in cases in which the user intends to
|
||||
* perform decompositions via LU::compute(const MatrixType&).
|
||||
*/
|
||||
* \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.
|
||||
@ -167,10 +158,10 @@ template<typename MatrixType> class FullPivLU
|
||||
*
|
||||
* \sa image()
|
||||
*/
|
||||
inline const ei_fullpivlu_kernel_impl<MatrixType> kernel() const
|
||||
inline const ei_kernel_return_value<FullPivLU> kernel() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "LU is not initialized.");
|
||||
return ei_fullpivlu_kernel_impl<MatrixType>(*this);
|
||||
return ei_kernel_return_value<FullPivLU>(*this);
|
||||
}
|
||||
|
||||
/** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
|
||||
@ -192,12 +183,11 @@ template<typename MatrixType> class FullPivLU
|
||||
*
|
||||
* \sa kernel()
|
||||
*/
|
||||
template<typename OriginalMatrixType>
|
||||
inline const ei_fullpivlu_image_impl<MatrixType>
|
||||
image(const MatrixBase<OriginalMatrixType>& originalMatrix) const
|
||||
inline const ei_image_return_value<FullPivLU>
|
||||
image(const MatrixType& originalMatrix) const
|
||||
{
|
||||
ei_assert(m_isInitialized && "LU is not initialized.");
|
||||
return ei_fullpivlu_image_impl<MatrixType>(*this, originalMatrix.derived());
|
||||
return ei_image_return_value<FullPivLU>(*this, originalMatrix);
|
||||
}
|
||||
|
||||
/** \return a solution x to the equation Ax=b, where A is the matrix of which
|
||||
@ -220,11 +210,11 @@ template<typename MatrixType> class FullPivLU
|
||||
* \sa TriangularView::solve(), kernel(), inverse()
|
||||
*/
|
||||
template<typename Rhs>
|
||||
inline const ei_fullpivlu_solve_impl<MatrixType, Rhs>
|
||||
inline const ei_solve_return_value<FullPivLU, Rhs>
|
||||
solve(const MatrixBase<Rhs>& b) const
|
||||
{
|
||||
ei_assert(m_isInitialized && "LU is not initialized.");
|
||||
return ei_fullpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived());
|
||||
return ei_solve_return_value<FullPivLU, Rhs>(*this, b.derived());
|
||||
}
|
||||
|
||||
/** \returns the determinant of the matrix of which
|
||||
@ -365,14 +355,17 @@ template<typename MatrixType> class FullPivLU
|
||||
*
|
||||
* \sa MatrixBase::inverse()
|
||||
*/
|
||||
inline const ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
|
||||
inline const ei_solve_return_value<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_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
|
||||
return ei_solve_return_value<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;
|
||||
@ -492,42 +485,21 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
|
||||
|
||||
/********* Implementation of kernel() **************************************************/
|
||||
|
||||
template<typename MatrixType>
|
||||
struct ei_traits<ei_fullpivlu_kernel_impl<MatrixType> >
|
||||
template<typename MatrixType, typename Dest>
|
||||
struct ei_kernel_impl<FullPivLU<MatrixType>, Dest>
|
||||
: ei_kernel_return_value<FullPivLU<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 MatrixType>
|
||||
struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl<MatrixType> >
|
||||
{
|
||||
typedef FullPivLU<MatrixType> LUType;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
const LUType& m_lu;
|
||||
int m_rank, m_cols;
|
||||
|
||||
ei_fullpivlu_kernel_impl(const LUType& lu)
|
||||
: m_lu(lu),
|
||||
m_rank(lu.rank()),
|
||||
m_cols(m_rank==lu.matrixLU().cols() ? 1 : lu.matrixLU().cols() - m_rank){}
|
||||
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
|
||||
MatrixType::MaxColsAtCompileTime,
|
||||
MatrixType::MaxRowsAtCompileTime)
|
||||
};
|
||||
|
||||
inline int rows() const { return m_lu.matrixLU().cols(); }
|
||||
inline int cols() const { return m_cols; }
|
||||
|
||||
template<typename Dest> void evalTo(Dest& dst) const
|
||||
void evalTo(Dest& dst) const
|
||||
{
|
||||
const int cols = m_lu.matrixLU().cols(), dimker = cols - m_rank;
|
||||
const FullPivLU<MatrixType>& dec = this->m_dec;
|
||||
const int cols = dec.matrixLU().cols(), rank = this->m_rank, dimker = cols - rank;
|
||||
if(dimker == 0)
|
||||
{
|
||||
// The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
|
||||
@ -549,89 +521,73 @@ struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl<
|
||||
*
|
||||
* 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 dimKer linearly
|
||||
* dimKer zero's. Let us use that to construct dimKer linearly
|
||||
* independent vectors in Ker U.
|
||||
*/
|
||||
|
||||
Matrix<int, Dynamic, 1, 0, LUType::MaxSmallDimAtCompileTime, 1> pivots(m_rank);
|
||||
RealScalar premultiplied_threshold = m_lu.maxPivot() * m_lu.threshold();
|
||||
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank);
|
||||
RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold();
|
||||
int p = 0;
|
||||
for(int i = 0; i < m_lu.nonzeroPivots(); ++i)
|
||||
if(ei_abs(m_lu.matrixLU().coeff(i,i)) > premultiplied_threshold)
|
||||
for(int i = 0; i < dec.nonzeroPivots(); ++i)
|
||||
if(ei_abs(dec.matrixLU().coeff(i,i)) > premultiplied_threshold)
|
||||
pivots.coeffRef(p++) = i;
|
||||
ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!");
|
||||
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,
|
||||
LUType::MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
|
||||
m(m_lu.matrixLU().block(0, 0, m_rank, cols));
|
||||
for(int i = 0; i < m_rank; ++i)
|
||||
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) = m_lu.matrixLU().row(pivots.coeff(i)).end(cols-i);
|
||||
m.row(i).end(cols-i) = dec.matrixLU().row(pivots.coeff(i)).end(cols-i);
|
||||
}
|
||||
m.block(0, 0, m_rank, m_rank).template triangularView<StrictlyLowerTriangular>().setZero();
|
||||
for(int i = 0; i < m_rank; ++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, m_rank, m_rank)
|
||||
m.corner(TopLeft, rank, rank)
|
||||
.template triangularView<UpperTriangular>().solveInPlace(
|
||||
m.corner(TopRight, m_rank, dimker)
|
||||
m.corner(TopRight, rank, dimker)
|
||||
);
|
||||
|
||||
// now we must undo the column permutation that we had applied!
|
||||
for(int i = m_rank-1; i >= 0; --i)
|
||||
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 < m_rank; ++i) dst.row(m_lu.permutationQ().coeff(i)) = -m.row(i).end(dimker);
|
||||
for(int i = m_rank; i < cols; ++i) dst.row(m_lu.permutationQ().coeff(i)).setZero();
|
||||
for(int k = 0; k < dimker; ++k) dst.coeffRef(m_lu.permutationQ().coeff(m_rank+k), k) = Scalar(1);
|
||||
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_traits<ei_fullpivlu_image_impl<MatrixType> >
|
||||
template<typename MatrixType, typename Dest>
|
||||
struct ei_image_impl<FullPivLU<MatrixType>, Dest>
|
||||
: ei_image_return_value<FullPivLU<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 MatrixType>
|
||||
struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<MatrixType> >
|
||||
{
|
||||
typedef FullPivLU<MatrixType> LUType;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
const LUType& m_lu;
|
||||
int m_rank, m_cols;
|
||||
const MatrixType& m_originalMatrix;
|
||||
|
||||
ei_fullpivlu_image_impl(const LUType& lu, const MatrixType& originalMatrix)
|
||||
: m_lu(lu), m_rank(lu.rank()),
|
||||
m_cols(m_rank == 0 ? 1 : m_rank),
|
||||
m_originalMatrix(originalMatrix) {}
|
||||
enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
|
||||
MatrixType::MaxColsAtCompileTime,
|
||||
MatrixType::MaxRowsAtCompileTime)
|
||||
};
|
||||
|
||||
inline int rows() const { return m_lu.matrixLU().rows(); }
|
||||
inline int cols() const { return m_cols; }
|
||||
|
||||
template<typename Dest> void evalTo(Dest& dst) const
|
||||
void evalTo(Dest& dst) const
|
||||
{
|
||||
if(m_rank == 0)
|
||||
const int rank = this->m_rank;
|
||||
const FullPivLU<MatrixType>& dec = this->m_dec;
|
||||
const MatrixType& originalMatrix = this->m_originalMatrix;
|
||||
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
|
||||
@ -640,61 +596,40 @@ struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<Ma
|
||||
return;
|
||||
}
|
||||
|
||||
Matrix<int, Dynamic, 1, 0, LUType::MaxSmallDimAtCompileTime, 1> pivots(m_rank);
|
||||
RealScalar premultiplied_threshold = m_lu.maxPivot() * m_lu.threshold();
|
||||
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank);
|
||||
RealScalar premultiplied_threshold = dec.maxPivot() * dec.threshold();
|
||||
int p = 0;
|
||||
for(int i = 0; i < m_lu.nonzeroPivots(); ++i)
|
||||
if(ei_abs(m_lu.matrixLU().coeff(i,i)) > premultiplied_threshold)
|
||||
for(int i = 0; i < dec.nonzeroPivots(); ++i)
|
||||
if(ei_abs(dec.matrixLU().coeff(i,i)) > premultiplied_threshold)
|
||||
pivots.coeffRef(p++) = i;
|
||||
ei_assert(p == m_rank && "You hit a bug in Eigen! Please report (backtrace and matrix)!");
|
||||
ei_internal_assert(p == rank);
|
||||
|
||||
for(int i = 0; i < m_rank; ++i)
|
||||
dst.col(i) = m_originalMatrix.col(m_lu.permutationQ().coeff(pivots.coeff(i)));
|
||||
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_traits<ei_fullpivlu_solve_impl<MatrixType,Rhs> >
|
||||
template<typename MatrixType, typename Rhs, typename Dest>
|
||||
struct ei_solve_impl<FullPivLU<MatrixType>, Rhs, Dest>
|
||||
: ei_solve_return_value<FullPivLU<MatrixType>, Rhs>
|
||||
{
|
||||
typedef Matrix<typename Rhs::Scalar,
|
||||
MatrixType::ColsAtCompileTime,
|
||||
Rhs::ColsAtCompileTime,
|
||||
Rhs::PlainMatrixType::Options,
|
||||
MatrixType::MaxColsAtCompileTime,
|
||||
Rhs::MaxColsAtCompileTime> ReturnMatrixType;
|
||||
};
|
||||
|
||||
template<typename MatrixType, typename Rhs>
|
||||
struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<MatrixType, Rhs> >
|
||||
{
|
||||
typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
|
||||
typedef FullPivLU<MatrixType> LUType;
|
||||
const LUType& m_lu;
|
||||
const typename Rhs::Nested m_rhs;
|
||||
|
||||
ei_fullpivlu_solve_impl(const LUType& lu, const Rhs& rhs)
|
||||
: m_lu(lu), m_rhs(rhs)
|
||||
{}
|
||||
|
||||
inline int rows() const { return m_lu.matrixLU().cols(); }
|
||||
inline int cols() const { return m_rhs.cols(); }
|
||||
|
||||
template<typename Dest> void evalTo(Dest& dst) const
|
||||
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;
|
||||
*/
|
||||
* 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 = m_lu.matrixLU().rows(),
|
||||
cols = m_lu.matrixLU().cols(),
|
||||
nonzero_pivots = m_lu.nonzeroPivots();
|
||||
ei_assert(m_rhs.rows() == rows);
|
||||
const FullPivLU<MatrixType>& dec = this->m_dec;
|
||||
const Rhs& rhs = this->m_rhs;
|
||||
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)
|
||||
@ -702,36 +637,36 @@ struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<Ma
|
||||
dst.setZero();
|
||||
return;
|
||||
}
|
||||
|
||||
typename Rhs::PlainMatrixType c(m_rhs.rows(), m_rhs.cols());
|
||||
|
||||
typename Rhs::PlainMatrixType c(rhs.rows(), rhs.cols());
|
||||
|
||||
// Step 1
|
||||
for(int i = 0; i < rows; ++i)
|
||||
c.row(m_lu.permutationP().coeff(i)) = m_rhs.row(i);
|
||||
c.row(dec.permutationP().coeff(i)) = rhs.row(i);
|
||||
|
||||
// Step 2
|
||||
m_lu.matrixLU()
|
||||
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())
|
||||
-= m_lu.matrixLU().corner(Eigen::BottomLeft, rows-cols, cols)
|
||||
* c.corner(Eigen::TopLeft, cols, c.cols());
|
||||
-= dec.matrixLU().corner(Eigen::BottomLeft, rows-cols, cols)
|
||||
* c.corner(Eigen::TopLeft, cols, c.cols());
|
||||
}
|
||||
|
||||
// Step 3
|
||||
m_lu.matrixLU()
|
||||
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(m_lu.permutationQ().coeff(i)) = c.row(i);
|
||||
for(int i = nonzero_pivots; i < m_lu.matrixLU().cols(); ++i)
|
||||
dst.row(m_lu.permutationQ().coeff(i)).setZero();
|
||||
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();
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
6
Eigen/src/misc/CMakeLists.txt
Normal file
6
Eigen/src/misc/CMakeLists.txt
Normal 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
|
||||
)
|
||||
72
Eigen/src/misc/Image.h
Normal file
72
Eigen/src/misc/Image.h
Normal file
@ -0,0 +1,72 @@
|
||||
// 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_return_value
|
||||
*
|
||||
*/
|
||||
template<typename DecompositionType>
|
||||
struct ei_traits<ei_image_return_value<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_return_value
|
||||
: public ReturnByValue<ei_image_return_value<_DecompositionType> >
|
||||
{
|
||||
typedef _DecompositionType DecompositionType;
|
||||
typedef typename DecompositionType::MatrixType MatrixType;
|
||||
|
||||
const DecompositionType& m_dec;
|
||||
int m_rank, m_cols;
|
||||
const MatrixType& m_originalMatrix;
|
||||
|
||||
ei_image_return_value(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; }
|
||||
|
||||
template<typename Dest> inline void evalTo(Dest& dst) const
|
||||
{
|
||||
static_cast<const ei_image_impl<DecompositionType, Dest> *>
|
||||
(this)->evalTo(dst);
|
||||
}
|
||||
};
|
||||
|
||||
#endif // EIGEN_MISC_IMAGE_H
|
||||
71
Eigen/src/misc/Kernel.h
Normal file
71
Eigen/src/misc/Kernel.h
Normal file
@ -0,0 +1,71 @@
|
||||
// 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_return_value
|
||||
*
|
||||
*/
|
||||
template<typename DecompositionType>
|
||||
struct ei_traits<ei_kernel_return_value<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_return_value
|
||||
: public ReturnByValue<ei_kernel_return_value<_DecompositionType> >
|
||||
{
|
||||
typedef _DecompositionType DecompositionType;
|
||||
const DecompositionType& m_dec;
|
||||
int m_rank, m_cols;
|
||||
|
||||
ei_kernel_return_value(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; }
|
||||
|
||||
template<typename Dest> inline void evalTo(Dest& dst) const
|
||||
{
|
||||
static_cast<const ei_kernel_impl<DecompositionType, Dest> *>
|
||||
(this)->evalTo(dst);
|
||||
}
|
||||
};
|
||||
|
||||
#endif // EIGEN_MISC_KERNEL_H
|
||||
65
Eigen/src/misc/Solve.h
Normal file
65
Eigen/src/misc/Solve.h
Normal file
@ -0,0 +1,65 @@
|
||||
// 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_return_value
|
||||
*
|
||||
*/
|
||||
template<typename DecompositionType, typename Rhs>
|
||||
struct ei_traits<ei_solve_return_value<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_return_value
|
||||
: public ReturnByValue<ei_solve_return_value<_DecompositionType, Rhs> >
|
||||
{
|
||||
typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNestedCleaned;
|
||||
typedef _DecompositionType DecompositionType;
|
||||
const DecompositionType& m_dec;
|
||||
const typename Rhs::Nested m_rhs;
|
||||
|
||||
ei_solve_return_value(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(); }
|
||||
|
||||
template<typename Dest> inline void evalTo(Dest& dst) const
|
||||
{
|
||||
static_cast<const ei_solve_impl<DecompositionType, RhsNestedCleaned, Dest> *>
|
||||
(this)->evalTo(dst);
|
||||
}
|
||||
};
|
||||
|
||||
#endif // EIGEN_MISC_SOLVE_H
|
||||
@ -49,8 +49,8 @@ template<typename MatrixType> void lu_non_invertible()
|
||||
cols2 = cols = MatrixType::ColsAtCompileTime;
|
||||
}
|
||||
|
||||
typedef typename ei_fullpivlu_kernel_impl<MatrixType>::ReturnMatrixType KernelMatrixType;
|
||||
typedef typename ei_fullpivlu_image_impl <MatrixType>::ReturnMatrixType ImageMatrixType;
|
||||
typedef typename ei_kernel_return_value<FullPivLU<MatrixType> >::ReturnMatrixType KernelMatrixType;
|
||||
typedef typename ei_image_return_value<FullPivLU<MatrixType> >::ReturnMatrixType ImageMatrixType;
|
||||
typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType;
|
||||
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime>
|
||||
CMatrixType;
|
||||
@ -65,10 +65,10 @@ template<typename MatrixType> void lu_non_invertible()
|
||||
createRandomMatrixOfRank(rank, rows, cols, m1);
|
||||
|
||||
FullPivLU<MatrixType> lu(m1);
|
||||
std::cout << lu.kernel().rows() << " " << lu.kernel().cols() << std::endl;
|
||||
KernelMatrixType m1kernel = lu.kernel();
|
||||
ImageMatrixType m1image = lu.image(m1);
|
||||
ImageMatrixType m1image = lu.image(m1);
|
||||
|
||||
// std::cerr << rank << " " << lu.rank() << std::endl;
|
||||
VERIFY(rank == lu.rank());
|
||||
VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
|
||||
VERIFY(!lu.isInjective());
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user