Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 26bc714..03b6af7 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -10,7 +10,18 @@
#ifndef EIGEN_LU_H
#define EIGEN_LU_H
-namespace Eigen {
+namespace Eigen {
+
+namespace internal {
+template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
+ : traits<_MatrixType>
+{
+ typedef MatrixXpr XprKind;
+ typedef SolverStorage StorageKind;
+ enum { Flags = 0 };
+};
+
+} // end namespace internal
/** \ingroup LU_Module
*
@@ -18,7 +29,7 @@
*
* \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
+ * \tparam _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 \f$ A = P^{-1} L U Q^{-1} \f$ where L is unit-lower-triangular, U is
@@ -41,27 +52,28 @@
* \include class_FullPivLU.cpp
* Output: \verbinclude class_FullPivLU.out
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
*/
template<typename _MatrixType> class FullPivLU
+ : public SolverBase<FullPivLU<_MatrixType> >
{
public:
typedef _MatrixType MatrixType;
+ typedef SolverBase<FullPivLU> Base;
+
+ EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU)
+ // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
enum {
- RowsAtCompileTime = MatrixType::RowsAtCompileTime,
- ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
- typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
- typedef typename MatrixType::Index Index;
- typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
- typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
+ typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
+ typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
+ typedef typename MatrixType::PlainObject PlainObject;
/**
* \brief Default Constructor.
@@ -84,7 +96,17 @@
* \param matrix the matrix of which to compute the LU decomposition.
* It is required to be nonzero.
*/
- FullPivLU(const MatrixType& matrix);
+ template<typename InputType>
+ explicit FullPivLU(const EigenBase<InputType>& matrix);
+
+ /** \brief Constructs a LU factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
+ *
+ * \sa FullPivLU(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit FullPivLU(EigenBase<InputType>& matrix);
/** Computes the LU decomposition of the given matrix.
*
@@ -93,7 +115,12 @@
*
* \returns a reference to *this
*/
- FullPivLU& compute(const MatrixType& matrix);
+ template<typename InputType>
+ FullPivLU& compute(const EigenBase<InputType>& matrix) {
+ m_lu = matrix.derived();
+ computeInPlace();
+ return *this;
+ }
/** \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
@@ -129,7 +156,7 @@
*
* \sa permutationQ()
*/
- inline const PermutationPType& permutationP() const
+ EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
return m_p;
@@ -166,7 +193,7 @@
}
/** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
- * will form a basis of the kernel.
+ * will form a basis of the image (column-space).
*
* \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
@@ -210,12 +237,22 @@
*
* \sa TriangularView::solve(), kernel(), inverse()
*/
+ // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
template<typename Rhs>
- inline const internal::solve_retval<FullPivLU, Rhs>
+ inline const Solve<FullPivLU, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
- return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
+ return Solve<FullPivLU, Rhs>(*this, b.derived());
+ }
+
+ /** \returns an estimate of the reciprocal condition number of the matrix of which \c *this is
+ the LU decomposition.
+ */
+ inline RealScalar rcond() const
+ {
+ eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
+ return internal::rcond_estimate_helper(m_l1_norm, *this);
}
/** \returns the determinant of the matrix of which
@@ -360,33 +397,46 @@
*
* \sa MatrixBase::inverse()
*/
- inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
+ inline const Inverse<FullPivLU> inverse() const
{
eigen_assert(m_isInitialized && "LU is not initialized.");
eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
- return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
- (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
+ return Inverse<FullPivLU>(*this);
}
MatrixType reconstructedMatrix() const;
- inline Index rows() const { return m_lu.rows(); }
- inline Index cols() const { return m_lu.cols(); }
+ EIGEN_DEVICE_FUNC inline Index rows() const { return m_lu.rows(); }
+ EIGEN_DEVICE_FUNC inline Index cols() const { return m_lu.cols(); }
+
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ template<typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl(const RhsType &rhs, DstType &dst) const;
+
+ template<bool Conjugate, typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
+ #endif
protected:
-
+
static void check_template_parameters()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
-
+
+ void computeInPlace();
+
MatrixType m_lu;
PermutationPType m_p;
PermutationQType m_q;
IntColVectorType m_rowsTranspositions;
IntRowVectorType m_colsTranspositions;
- Index m_det_pq, m_nonzero_pivots;
+ Index m_nonzero_pivots;
+ RealScalar m_l1_norm;
RealScalar m_maxpivot, m_prescribedThreshold;
+ signed char m_det_pq;
bool m_isInitialized, m_usePrescribedThreshold;
};
@@ -409,7 +459,8 @@
}
template<typename MatrixType>
-FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
+template<typename InputType>
+FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix)
: m_lu(matrix.rows(), matrix.cols()),
m_p(matrix.rows()),
m_q(matrix.cols()),
@@ -418,28 +469,41 @@
m_isInitialized(false),
m_usePrescribedThreshold(false)
{
- compute(matrix);
+ compute(matrix.derived());
}
template<typename MatrixType>
-FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
+template<typename InputType>
+FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix)
+ : m_lu(matrix.derived()),
+ m_p(matrix.rows()),
+ m_q(matrix.cols()),
+ m_rowsTranspositions(matrix.rows()),
+ m_colsTranspositions(matrix.cols()),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false)
+{
+ computeInPlace();
+}
+
+template<typename MatrixType>
+void FullPivLU<MatrixType>::computeInPlace()
{
check_template_parameters();
-
- // the permutations are stored as int indices, so just to be sure:
- eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
-
- m_isInitialized = true;
- m_lu = matrix;
- const Index size = matrix.diagonalSize();
- const Index rows = matrix.rows();
- const Index cols = matrix.cols();
+ // the permutations are stored as int indices, so just to be sure:
+ eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest());
+
+ m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
+
+ const Index size = m_lu.diagonalSize();
+ const Index rows = m_lu.rows();
+ const Index cols = m_lu.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.
- m_rowsTranspositions.resize(matrix.rows());
- m_colsTranspositions.resize(matrix.cols());
+ m_rowsTranspositions.resize(m_lu.rows());
+ m_colsTranspositions.resize(m_lu.cols());
Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
@@ -451,14 +515,16 @@
// biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
Index row_of_biggest_in_corner, col_of_biggest_in_corner;
- RealScalar biggest_in_corner;
+ typedef internal::scalar_score_coeff_op<Scalar> Scoring;
+ typedef typename Scoring::result_type Score;
+ Score biggest_in_corner;
biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
- .cwiseAbs()
+ .unaryExpr(Scoring())
.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(biggest_in_corner==RealScalar(0))
+ if(biggest_in_corner==Score(0))
{
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.
@@ -471,7 +537,8 @@
break;
}
- if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
+ RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
+ if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
// Now that we've found the pivot, we need to apply the row/col swaps to
// bring it to the location (k,k).
@@ -508,7 +575,8 @@
m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
- return *this;
+
+ m_isInitialized = true;
}
template<typename MatrixType>
@@ -671,64 +739,136 @@
/***** Implementation of solve() *****************************************************/
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<FullPivLU<_MatrixType>, Rhs>
- : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
+} // end namespace internal
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename _MatrixType>
+template<typename RhsType, typename DstType>
+void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
- EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
+ /* 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;
+ */
- template<typename Dest> void evalTo(Dest& dst) const
+ const Index rows = this->rows(),
+ cols = this->cols(),
+ nonzero_pivots = this->rank();
+ eigen_assert(rhs.rows() == rows);
+ const Index smalldim = (std::min)(rows, cols);
+
+ if(nonzero_pivots == 0)
{
- /* 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;
- */
+ dst.setZero();
+ return;
+ }
- const Index rows = dec().rows(), cols = dec().cols(),
- nonzero_pivots = dec().nonzeroPivots();
- eigen_assert(rhs().rows() == rows);
- const Index smalldim = (std::min)(rows, cols);
+ typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
- if(nonzero_pivots == 0)
- {
- dst.setZero();
- return;
- }
+ // Step 1
+ c = permutationP() * rhs;
- typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
+ // Step 2
+ m_lu.topLeftCorner(smalldim,smalldim)
+ .template triangularView<UnitLower>()
+ .solveInPlace(c.topRows(smalldim));
+ if(rows>cols)
+ c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
- // Step 1
- c = dec().permutationP() * rhs();
+ // Step 3
+ m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
+ .template triangularView<Upper>()
+ .solveInPlace(c.topRows(nonzero_pivots));
+ // Step 4
+ for(Index i = 0; i < nonzero_pivots; ++i)
+ dst.row(permutationQ().indices().coeff(i)) = c.row(i);
+ for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
+ dst.row(permutationQ().indices().coeff(i)).setZero();
+}
+
+template<typename _MatrixType>
+template<bool Conjugate, typename RhsType, typename DstType>
+void FullPivLU<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
+{
+ /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
+ * and since permutations are real and unitary, we can write this
+ * as A^T = Q U^T L^T P,
+ * So we proceed as follows:
+ * Step 1: compute c = Q^T rhs.
+ * Step 2: replace c by the solution x to U^T x = c. May or may not exist.
+ * Step 3: replace c by the solution x to L^T x = c.
+ * Step 4: result = P^T c.
+ * If Conjugate is true, replace "^T" by "^*" above.
+ */
+
+ const Index rows = this->rows(), cols = this->cols(),
+ nonzero_pivots = this->rank();
+ eigen_assert(rhs.rows() == cols);
+ const Index smalldim = (std::min)(rows, cols);
+
+ if(nonzero_pivots == 0)
+ {
+ dst.setZero();
+ return;
+ }
+
+ typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
+
+ // Step 1
+ c = permutationQ().inverse() * rhs;
+
+ if (Conjugate) {
// Step 2
- dec().matrixLU()
- .topLeftCorner(smalldim,smalldim)
- .template triangularView<UnitLower>()
- .solveInPlace(c.topRows(smalldim));
- if(rows>cols)
- {
- c.bottomRows(rows-cols)
- -= dec().matrixLU().bottomRows(rows-cols)
- * c.topRows(cols);
- }
-
- // Step 3
- dec().matrixLU()
- .topLeftCorner(nonzero_pivots, nonzero_pivots)
+ m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
+ .adjoint()
.solveInPlace(c.topRows(nonzero_pivots));
+ // Step 3
+ m_lu.topLeftCorner(smalldim, smalldim)
+ .template triangularView<UnitLower>()
+ .adjoint()
+ .solveInPlace(c.topRows(smalldim));
+ } else {
+ // Step 2
+ m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
+ .template triangularView<Upper>()
+ .transpose()
+ .solveInPlace(c.topRows(nonzero_pivots));
+ // Step 3
+ m_lu.topLeftCorner(smalldim, smalldim)
+ .template triangularView<UnitLower>()
+ .transpose()
+ .solveInPlace(c.topRows(smalldim));
+ }
- // Step 4
- for(Index i = 0; i < nonzero_pivots; ++i)
- dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
- for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
- dst.row(dec().permutationQ().indices().coeff(i)).setZero();
+ // Step 4
+ PermutationPType invp = permutationP().inverse().eval();
+ for(Index i = 0; i < smalldim; ++i)
+ dst.row(invp.indices().coeff(i)) = c.row(i);
+ for(Index i = smalldim; i < rows; ++i)
+ dst.row(invp.indices().coeff(i)).setZero();
+}
+
+#endif
+
+namespace internal {
+
+
+/***** Implementation of inverse() *****************************************************/
+template<typename DstXprType, typename MatrixType>
+struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense>
+{
+ typedef FullPivLU<MatrixType> LuType;
+ typedef Inverse<LuType> SrcXprType;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &)
+ {
+ dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
}
};
-
} // end namespace internal
/******* MatrixBase methods *****************************************************************/