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/QR/CMakeLists.txt b/Eigen/src/QR/CMakeLists.txt
deleted file mode 100644
index 96f43d7..0000000
--- a/Eigen/src/QR/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_QR_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_QR_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/QR COMPONENT Devel
- )
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 567eab7..a7b47d5 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -11,7 +11,16 @@
#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
-namespace Eigen {
+namespace Eigen {
+
+namespace internal {
+template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
+ : traits<_MatrixType>
+{
+ enum { Flags = 0 };
+};
+
+} // end namespace internal
/** \ingroup QR_Module
*
@@ -19,19 +28,21 @@
*
* \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
*
- * \param MatrixType the type of the matrix of which we are computing the QR decomposition
+ * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
- * such that
+ * such that
* \f[
* \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
* \f]
- * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
+ * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
* upper triangular matrix.
*
* This decomposition performs column pivoting in order to be rank-revealing and improve
* numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::colPivHouseholderQr()
*/
template<typename _MatrixType> class ColPivHouseholderQR
@@ -42,25 +53,25 @@
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
+ // FIXME should be int
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
-
+ typedef typename MatrixType::PlainObject PlainObject;
+
private:
-
- typedef typename PermutationType::Index PermIndexType;
-
+
+ typedef typename PermutationType::StorageIndex PermIndexType;
+
public:
/**
@@ -75,7 +86,8 @@
m_colsPermutation(),
m_colsTranspositions(),
m_temp(),
- m_colSqNorms(),
+ m_colNormsUpdated(),
+ m_colNormsDirect(),
m_isInitialized(false),
m_usePrescribedThreshold(false) {}
@@ -91,7 +103,8 @@
m_colsPermutation(PermIndexType(cols)),
m_colsTranspositions(cols),
m_temp(cols),
- m_colSqNorms(cols),
+ m_colNormsUpdated(cols),
+ m_colNormsDirect(cols),
m_isInitialized(false),
m_usePrescribedThreshold(false) {}
@@ -99,25 +112,48 @@
*
* This constructor computes the QR factorization of the matrix \a matrix by calling
* the method compute(). It is a short cut for:
- *
+ *
* \code
* ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
* qr.compute(matrix);
* \endcode
- *
+ *
* \sa compute()
*/
- ColPivHouseholderQR(const MatrixType& matrix)
+ template<typename InputType>
+ explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
m_colsPermutation(PermIndexType(matrix.cols())),
m_colsTranspositions(matrix.cols()),
m_temp(matrix.cols()),
- m_colSqNorms(matrix.cols()),
+ m_colNormsUpdated(matrix.cols()),
+ m_colNormsDirect(matrix.cols()),
m_isInitialized(false),
m_usePrescribedThreshold(false)
{
- compute(matrix);
+ compute(matrix.derived());
+ }
+
+ /** \brief Constructs a QR factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
+ *
+ * \sa ColPivHouseholderQR(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
+ : m_qr(matrix.derived()),
+ m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
+ m_colsPermutation(PermIndexType(matrix.cols())),
+ m_colsTranspositions(matrix.cols()),
+ m_temp(matrix.cols()),
+ m_colNormsUpdated(matrix.cols()),
+ m_colNormsDirect(matrix.cols()),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false)
+ {
+ computeInPlace();
}
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
@@ -127,9 +163,6 @@
*
* \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
@@ -138,17 +171,17 @@
* Output: \verbinclude ColPivHouseholderQR_solve.out
*/
template<typename Rhs>
- inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
+ inline const Solve<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
+ return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
- HouseholderSequenceType householderQ(void) const;
- HouseholderSequenceType matrixQ(void) const
+ HouseholderSequenceType householderQ() const;
+ HouseholderSequenceType matrixQ() const
{
- return householderQ();
+ return householderQ();
}
/** \returns a reference to the matrix where the Householder QR decomposition is stored
@@ -158,14 +191,14 @@
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_qr;
}
-
- /** \returns a reference to the matrix where the result Householder QR is stored
- * \warning The strict lower part of this matrix contains internal values.
+
+ /** \returns a reference to the matrix where the result Householder QR is stored
+ * \warning The strict lower part of this matrix contains internal values.
* Only the upper triangular part should be referenced. To get it, use
* \code matrixR().template triangularView<Upper>() \endcode
- * For rank-deficient matrices, use
- * \code
- * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
+ * For rank-deficient matrices, use
+ * \code
+ * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
* \endcode
*/
const MatrixType& matrixR() const
@@ -173,8 +206,9 @@
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return m_qr;
}
-
- ColPivHouseholderQR& compute(const MatrixType& matrix);
+
+ template<typename InputType>
+ ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
/** \returns a const reference to the column permutation matrix */
const PermutationType& colsPermutation() const
@@ -284,20 +318,17 @@
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*/
- inline const
- internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
- inverse() const
+ inline const Inverse<ColPivHouseholderQR> inverse() const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
- (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
+ return Inverse<ColPivHouseholderQR>(*this);
}
inline Index rows() const { return m_qr.rows(); }
inline Index cols() const { return m_qr.cols(); }
-
+
/** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
- *
+ *
* For advanced uses only.
*/
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
@@ -370,12 +401,12 @@
* diagonal coefficient of R.
*/
RealScalar maxPivot() const { return m_maxpivot; }
-
+
/** \brief Reports whether the QR factorization was succesful.
*
- * \note This function always returns \c Success. It is provided for compatibility
+ * \note This function always returns \c Success. It is provided for compatibility
* with other factorization routines.
- * \returns \c Success
+ * \returns \c Success
*/
ComputationInfo info() const
{
@@ -383,19 +414,30 @@
return Success;
}
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ template<typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl(const RhsType &rhs, DstType &dst) const;
+ #endif
+
protected:
-
+
+ friend class CompleteOrthogonalDecomposition<MatrixType>;
+
static void check_template_parameters()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
-
+
+ void computeInPlace();
+
MatrixType m_qr;
HCoeffsType m_hCoeffs;
PermutationType m_colsPermutation;
IntRowVectorType m_colsTranspositions;
RowVectorType m_temp;
- RealRowVectorType m_colSqNorms;
+ RealRowVectorType m_colNormsUpdated;
+ RealRowVectorType m_colNormsDirect;
bool m_isInitialized, m_usePrescribedThreshold;
RealScalar m_prescribedThreshold, m_maxpivot;
Index m_nonzero_pivots;
@@ -426,51 +468,57 @@
* \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
*/
template<typename MatrixType>
-ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
+template<typename InputType>
+ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
+{
+ m_qr = matrix.derived();
+ computeInPlace();
+ return *this;
+}
+
+template<typename MatrixType>
+void ColPivHouseholderQR<MatrixType>::computeInPlace()
{
check_template_parameters();
-
- using std::abs;
- Index rows = matrix.rows();
- Index cols = matrix.cols();
- Index size = matrix.diagonalSize();
-
- // the column permutation is stored as int indices, so just to be sure:
- eigen_assert(cols<=NumTraits<int>::highest());
- m_qr = matrix;
+ // the column permutation is stored as int indices, so just to be sure:
+ eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
+
+ using std::abs;
+
+ Index rows = m_qr.rows();
+ Index cols = m_qr.cols();
+ Index size = m_qr.diagonalSize();
+
m_hCoeffs.resize(size);
m_temp.resize(cols);
- m_colsTranspositions.resize(matrix.cols());
+ m_colsTranspositions.resize(m_qr.cols());
Index number_of_transpositions = 0;
- m_colSqNorms.resize(cols);
- for(Index k = 0; k < cols; ++k)
- m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
+ m_colNormsUpdated.resize(cols);
+ m_colNormsDirect.resize(cols);
+ for (Index k = 0; k < cols; ++k) {
+ // colNormsDirect(k) caches the most recent directly computed norm of
+ // column k.
+ m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
+ m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
+ }
- RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
+ RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
+ RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
for(Index k = 0; k < size; ++k)
{
- // first, we look up in our table m_colSqNorms which column has the biggest squared norm
+ // first, we look up in our table m_colNormsUpdated which column has the biggest norm
Index biggest_col_index;
- RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
+ RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
biggest_col_index += k;
- // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
- // the actual squared norm of the selected column.
- // Note that not doing so does result in solve() sometimes returning inf/nan values
- // when running the unit test with 1000 repetitions.
- biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
-
- // we store that back into our table: it can't hurt to correct our table.
- m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
-
// Track the number of meaningful pivots but do not stop the decomposition to make
// sure that the initial matrix is properly reproduced. See bug 941.
if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
@@ -480,7 +528,8 @@
m_colsTranspositions.coeffRef(k) = biggest_col_index;
if(k != biggest_col_index) {
m_qr.col(k).swap(m_qr.col(biggest_col_index));
- std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
+ std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
+ std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
++number_of_transpositions;
}
@@ -498,8 +547,28 @@
m_qr.bottomRightCorner(rows-k, cols-k-1)
.applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
- // update our table of squared norms of the columns
- m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
+ // update our table of norms of the columns
+ for (Index j = k + 1; j < cols; ++j) {
+ // The following implements the stable norm downgrade step discussed in
+ // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
+ // and used in LAPACK routines xGEQPF and xGEQP3.
+ // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
+ if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
+ RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
+ temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
+ temp = temp < RealScalar(0) ? RealScalar(0) : temp;
+ RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
+ m_colNormsDirect.coeffRef(j));
+ if (temp2 <= norm_downdate_threshold) {
+ // The updated norm has become too inaccurate so re-compute the column
+ // norm directly.
+ m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
+ m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
+ } else {
+ m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
+ }
+ }
+ }
}
m_colsPermutation.setIdentity(PermIndexType(cols));
@@ -508,46 +577,50 @@
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
-
- return *this;
}
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename _MatrixType>
+template<typename RhsType, typename DstType>
+void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
+{
+ eigen_assert(rhs.rows() == rows());
+
+ const Index nonzero_pivots = nonzeroPivots();
+
+ if(nonzero_pivots == 0)
+ {
+ dst.setZero();
+ return;
+ }
+
+ typename RhsType::PlainObject c(rhs);
+
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+ c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
+ .setLength(nonzero_pivots)
+ .transpose()
+ );
+
+ m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
+ .template triangularView<Upper>()
+ .solveInPlace(c.topRows(nonzero_pivots));
+
+ for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
+ for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
+}
+#endif
+
namespace internal {
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
- : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
+template<typename DstXprType, typename MatrixType>
+struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
{
- EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
+ typedef ColPivHouseholderQR<MatrixType> QrType;
+ typedef Inverse<QrType> SrcXprType;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
{
- eigen_assert(rhs().rows() == dec().rows());
-
- const Index cols = dec().cols(),
- nonzero_pivots = dec().nonzeroPivots();
-
- if(nonzero_pivots == 0)
- {
- dst.setZero();
- return;
- }
-
- typename Rhs::PlainObject c(rhs());
-
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
- c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
- .setLength(dec().nonzeroPivots())
- .transpose()
- );
-
- dec().matrixR()
- .topLeftCorner(nonzero_pivots, nonzero_pivots)
- .template triangularView<Upper>()
- .solveInPlace(c.topRows(nonzero_pivots));
-
- for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
- for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
+ dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
}
};
diff --git a/Eigen/src/QR/ColPivHouseholderQR_MKL.h b/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
similarity index 64%
rename from Eigen/src/QR/ColPivHouseholderQR_MKL.h
rename to Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
index b5b1983..4e9651f 100644
--- a/Eigen/src/QR/ColPivHouseholderQR_MKL.h
+++ b/Eigen/src/QR/ColPivHouseholderQR_LAPACKE.h
@@ -25,37 +25,34 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Householder QR decomposition of a matrix with column pivoting based on
* LAPACKE_?geqp3 function.
********************************************************************************
*/
-#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_MKL_H
-#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_MKL_H
-
-#include "Eigen/src/Core/util/MKL_support.h"
+#ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
+#define EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
namespace Eigen {
-/** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_QR_COLPIV(EIGTYPE, MKLTYPE, MKLPREFIX, EIGCOLROW, MKLCOLROW) \
-template<> inline \
+#define EIGEN_LAPACKE_QR_COLPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \
+template<> template<typename InputType> inline \
ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >& \
ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> >::compute( \
- const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix) \
+ const EigenBase<InputType>& matrix) \
\
{ \
using std::abs; \
typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \
- typedef MatrixType::Scalar Scalar; \
typedef MatrixType::RealScalar RealScalar; \
Index rows = matrix.rows();\
Index cols = matrix.cols();\
- Index size = matrix.diagonalSize();\
\
m_qr = matrix;\
+ Index size = m_qr.diagonalSize();\
m_hCoeffs.resize(size);\
\
m_colsTranspositions.resize(cols);\
@@ -66,34 +63,35 @@
m_colsPermutation.resize(cols); \
m_colsPermutation.indices().setZero(); \
\
- lapack_int lda = m_qr.outerStride(), i; \
- lapack_int matrix_order = MKLCOLROW; \
- LAPACKE_##MKLPREFIX##geqp3( matrix_order, rows, cols, (MKLTYPE*)m_qr.data(), lda, (lapack_int*)m_colsPermutation.indices().data(), (MKLTYPE*)m_hCoeffs.data()); \
+ lapack_int lda = internal::convert_index<lapack_int,Index>(m_qr.outerStride()); \
+ lapack_int matrix_order = LAPACKE_COLROW; \
+ LAPACKE_##LAPACKE_PREFIX##geqp3( matrix_order, internal::convert_index<lapack_int,Index>(rows), internal::convert_index<lapack_int,Index>(cols), \
+ (LAPACKE_TYPE*)m_qr.data(), lda, (lapack_int*)m_colsPermutation.indices().data(), (LAPACKE_TYPE*)m_hCoeffs.data()); \
m_isInitialized = true; \
m_maxpivot=m_qr.diagonal().cwiseAbs().maxCoeff(); \
m_hCoeffs.adjointInPlace(); \
RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); \
lapack_int *perm = m_colsPermutation.indices().data(); \
- for(i=0;i<size;i++) { \
+ for(Index i=0;i<size;i++) { \
m_nonzero_pivots += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);\
} \
- for(i=0;i<cols;i++) perm[i]--;\
+ for(Index i=0;i<cols;i++) perm[i]--;\
\
/*m_det_pq = (number_of_transpositions%2) ? -1 : 1; // TODO: It's not needed now; fix upon availability in Eigen */ \
\
return *this; \
}
-EIGEN_MKL_QR_COLPIV(double, double, d, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_QR_COLPIV(float, float, s, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_QR_COLPIV(dcomplex, MKL_Complex16, z, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_QR_COLPIV(scomplex, MKL_Complex8, c, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(double, double, d, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(float, float, s, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, ColMajor, LAPACK_COL_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float, c, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_MKL_QR_COLPIV(double, double, d, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_QR_COLPIV(float, float, s, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_QR_COLPIV(dcomplex, MKL_Complex16, z, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_MKL_QR_COLPIV(scomplex, MKL_Complex8, c, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(double, double, d, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(float, float, s, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(dcomplex, lapack_complex_double, z, RowMajor, LAPACK_ROW_MAJOR)
+EIGEN_LAPACKE_QR_COLPIV(scomplex, lapack_complex_float, c, RowMajor, LAPACK_ROW_MAJOR)
} // end namespace Eigen
-#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_MKL_H
+#endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_LAPACKE_H
diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
new file mode 100644
index 0000000..34c637b
--- /dev/null
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -0,0 +1,562 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
+#define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
+
+namespace Eigen {
+
+namespace internal {
+template <typename _MatrixType>
+struct traits<CompleteOrthogonalDecomposition<_MatrixType> >
+ : traits<_MatrixType> {
+ enum { Flags = 0 };
+};
+
+} // end namespace internal
+
+/** \ingroup QR_Module
+ *
+ * \class CompleteOrthogonalDecomposition
+ *
+ * \brief Complete orthogonal decomposition (COD) of a matrix.
+ *
+ * \param MatrixType the type of the matrix of which we are computing the COD.
+ *
+ * This class performs a rank-revealing complete orthogonal decomposition of a
+ * matrix \b A into matrices \b P, \b Q, \b T, and \b Z such that
+ * \f[
+ * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \,
+ * \begin{bmatrix} \mathbf{T} & \mathbf{0} \\
+ * \mathbf{0} & \mathbf{0} \end{bmatrix} \, \mathbf{Z}
+ * \f]
+ * by using Householder transformations. Here, \b P is a permutation matrix,
+ * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of
+ * size rank-by-rank. \b A may be rank deficient.
+ *
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
+ * \sa MatrixBase::completeOrthogonalDecomposition()
+ */
+template <typename _MatrixType>
+class CompleteOrthogonalDecomposition {
+ public:
+ typedef _MatrixType MatrixType;
+ enum {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+ };
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
+ typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
+ PermutationType;
+ typedef typename internal::plain_row_type<MatrixType, Index>::type
+ IntRowVectorType;
+ typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
+ typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
+ RealRowVectorType;
+ typedef HouseholderSequence<
+ MatrixType, typename internal::remove_all<
+ typename HCoeffsType::ConjugateReturnType>::type>
+ HouseholderSequenceType;
+ typedef typename MatrixType::PlainObject PlainObject;
+
+ private:
+ typedef typename PermutationType::Index PermIndexType;
+
+ public:
+ /**
+ * \brief Default Constructor.
+ *
+ * The default constructor is useful in cases in which the user intends to
+ * perform decompositions via
+ * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&).
+ */
+ CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
+
+ /** \brief Default Constructor with memory preallocation
+ *
+ * Like the default constructor but with preallocation of the internal data
+ * according to the specified problem \a size.
+ * \sa CompleteOrthogonalDecomposition()
+ */
+ CompleteOrthogonalDecomposition(Index rows, Index cols)
+ : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
+
+ /** \brief Constructs a complete orthogonal decomposition from a given
+ * matrix.
+ *
+ * This constructor computes the complete orthogonal decomposition of the
+ * matrix \a matrix by calling the method compute(). The default
+ * threshold for rank determination will be used. It is a short cut for:
+ *
+ * \code
+ * CompleteOrthogonalDecomposition<MatrixType> cod(matrix.rows(),
+ * matrix.cols());
+ * cod.setThreshold(Default);
+ * cod.compute(matrix);
+ * \endcode
+ *
+ * \sa compute()
+ */
+ template <typename InputType>
+ explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix)
+ : m_cpqr(matrix.rows(), matrix.cols()),
+ m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
+ m_temp(matrix.cols())
+ {
+ compute(matrix.derived());
+ }
+
+ /** \brief Constructs a complete orthogonal decomposition from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
+ *
+ * \sa CompleteOrthogonalDecomposition(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix)
+ : m_cpqr(matrix.derived()),
+ m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
+ m_temp(matrix.cols())
+ {
+ computeInPlace();
+ }
+
+
+ /** This method computes the minimum-norm solution X to a least squares
+ * problem \f[\mathrm{minimize} \|A X - B\|, \f] where \b A is the matrix of
+ * which \c *this is the complete orthogonal decomposition.
+ *
+ * \param b the right-hand sides of the problem to solve.
+ *
+ * \returns a solution.
+ *
+ */
+ template <typename Rhs>
+ inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve(
+ const MatrixBase<Rhs>& b) const {
+ eigen_assert(m_cpqr.m_isInitialized &&
+ "CompleteOrthogonalDecomposition is not initialized.");
+ return Solve<CompleteOrthogonalDecomposition, Rhs>(*this, b.derived());
+ }
+
+ HouseholderSequenceType householderQ(void) const;
+ HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
+
+ /** \returns the matrix \b Z.
+ */
+ MatrixType matrixZ() const {
+ MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
+ applyZAdjointOnTheLeftInPlace(Z);
+ return Z.adjoint();
+ }
+
+ /** \returns a reference to the matrix where the complete orthogonal
+ * decomposition is stored
+ */
+ const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
+
+ /** \returns a reference to the matrix where the complete orthogonal
+ * decomposition is stored.
+ * \warning The strict lower part and \code cols() - rank() \endcode right
+ * columns of this matrix contains internal values.
+ * Only the upper triangular part should be referenced. To get it, use
+ * \code matrixT().template triangularView<Upper>() \endcode
+ * For rank-deficient matrices, use
+ * \code
+ * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
+ * \endcode
+ */
+ const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
+
+ template <typename InputType>
+ CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
+ // Compute the column pivoted QR factorization A P = Q R.
+ m_cpqr.compute(matrix);
+ computeInPlace();
+ return *this;
+ }
+
+ /** \returns a const reference to the column permutation matrix */
+ const PermutationType& colsPermutation() const {
+ return m_cpqr.colsPermutation();
+ }
+
+ /** \returns the absolute value of the determinant of the matrix of which
+ * *this is the complete orthogonal decomposition. It has only linear
+ * complexity (that is, O(n) where n is the dimension of the square matrix)
+ * as the complete orthogonal decomposition has already been computed.
+ *
+ * \note This is only for square matrices.
+ *
+ * \warning a determinant can be very big or small, so for matrices
+ * of large enough dimension, there is a risk of overflow/underflow.
+ * One way to work around that is to use logAbsDeterminant() instead.
+ *
+ * \sa logAbsDeterminant(), MatrixBase::determinant()
+ */
+ typename MatrixType::RealScalar absDeterminant() const;
+
+ /** \returns the natural log of the absolute value of the determinant of the
+ * matrix of which *this is the complete orthogonal decomposition. It has
+ * only linear complexity (that is, O(n) where n is the dimension of the
+ * square matrix) as the complete orthogonal decomposition has already been
+ * computed.
+ *
+ * \note This is only for square matrices.
+ *
+ * \note This method is useful to work around the risk of overflow/underflow
+ * that's inherent to determinant computation.
+ *
+ * \sa absDeterminant(), MatrixBase::determinant()
+ */
+ typename MatrixType::RealScalar logAbsDeterminant() const;
+
+ /** \returns the rank of the matrix of which *this is the complete orthogonal
+ * 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 Index rank() const { return m_cpqr.rank(); }
+
+ /** \returns the dimension of the kernel of the matrix of which *this is the
+ * complete orthogonal 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 Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
+
+ /** \returns true if the matrix of which *this is the 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 { return m_cpqr.isInjective(); }
+
+ /** \returns true if the matrix of which *this is the 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 { return m_cpqr.isSurjective(); }
+
+ /** \returns true if the matrix of which *this is the complete orthogonal
+ * 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 { return m_cpqr.isInvertible(); }
+
+ /** \returns the pseudo-inverse of the matrix of which *this is the complete
+ * orthogonal decomposition.
+ * \warning: Do not compute \c this->pseudoInverse()*rhs to solve a linear systems.
+ * It is more efficient and numerically stable to call \c this->solve(rhs).
+ */
+ inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const
+ {
+ return Inverse<CompleteOrthogonalDecomposition>(*this);
+ }
+
+ inline Index rows() const { return m_cpqr.rows(); }
+ inline Index cols() const { return m_cpqr.cols(); }
+
+ /** \returns a const reference to the vector of Householder coefficients used
+ * to represent the factor \c Q.
+ *
+ * For advanced uses only.
+ */
+ inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
+
+ /** \returns a const reference to the vector of Householder coefficients
+ * used to represent the factor \c Z.
+ *
+ * For advanced uses only.
+ */
+ const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
+
+ /** 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.
+ * Most be called before calling compute().
+ *
+ * When it needs to get the threshold value, Eigen calls threshold(). By
+ * default, this uses a formula to automatically determine a reasonable
+ * threshold. 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)
+ */
+ CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) {
+ m_cpqr.setThreshold(threshold);
+ return *this;
+ }
+
+ /** 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 qr.setThreshold(Eigen::Default); \endcode
+ *
+ * See the documentation of setThreshold(const RealScalar&).
+ */
+ CompleteOrthogonalDecomposition& setThreshold(Default_t) {
+ m_cpqr.setThreshold(Default);
+ return *this;
+ }
+
+ /** Returns the threshold that will be used by certain methods such as rank().
+ *
+ * See the documentation of setThreshold(const RealScalar&).
+ */
+ RealScalar threshold() const { return m_cpqr.threshold(); }
+
+ /** \returns the number of nonzero pivots in the complete orthogonal
+ * 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 Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
+
+ /** \returns the absolute value of the biggest pivot, i.e. the biggest
+ * diagonal coefficient of R.
+ */
+ inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
+
+ /** \brief Reports whether the complete orthogonal decomposition was
+ * succesful.
+ *
+ * \note This function always returns \c Success. It is provided for
+ * compatibility
+ * with other factorization routines.
+ * \returns \c Success
+ */
+ ComputationInfo info() const {
+ eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
+ return Success;
+ }
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+ template <typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC void _solve_impl(const RhsType& rhs, DstType& dst) const;
+#endif
+
+ protected:
+ static void check_template_parameters() {
+ EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
+ }
+
+ void computeInPlace();
+
+ /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$.
+ */
+ template <typename Rhs>
+ void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
+
+ ColPivHouseholderQR<MatrixType> m_cpqr;
+ HCoeffsType m_zCoeffs;
+ RowVectorType m_temp;
+};
+
+template <typename MatrixType>
+typename MatrixType::RealScalar
+CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {
+ return m_cpqr.absDeterminant();
+}
+
+template <typename MatrixType>
+typename MatrixType::RealScalar
+CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
+ return m_cpqr.logAbsDeterminant();
+}
+
+/** Performs the complete orthogonal decomposition of the given matrix \a
+ * matrix. The result of the factorization is stored into \c *this, and a
+ * reference to \c *this is returned.
+ *
+ * \sa class CompleteOrthogonalDecomposition,
+ * CompleteOrthogonalDecomposition(const MatrixType&)
+ */
+template <typename MatrixType>
+void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
+{
+ check_template_parameters();
+
+ // the column permutation is stored as int indices, so just to be sure:
+ eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
+
+ const Index rank = m_cpqr.rank();
+ const Index cols = m_cpqr.cols();
+ const Index rows = m_cpqr.rows();
+ m_zCoeffs.resize((std::min)(rows, cols));
+ m_temp.resize(cols);
+
+ if (rank < cols) {
+ // We have reduced the (permuted) matrix to the form
+ // [R11 R12]
+ // [ 0 R22]
+ // where R11 is r-by-r (r = rank) upper triangular, R12 is
+ // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
+ // We now compute the complete orthogonal decomposition by applying
+ // Householder transformations from the right to the upper trapezoidal
+ // matrix X = [R11 R12] to zero out R12 and obtain the factorization
+ // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
+ // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
+ // We store the data representing Z in R12 and m_zCoeffs.
+ for (Index k = rank - 1; k >= 0; --k) {
+ if (k != rank - 1) {
+ // Given the API for Householder reflectors, it is more convenient if
+ // we swap the leading parts of columns k and r-1 (zero-based) to form
+ // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
+ m_cpqr.m_qr.col(k).head(k + 1).swap(
+ m_cpqr.m_qr.col(rank - 1).head(k + 1));
+ }
+ // Construct Householder reflector Z(k) to zero out the last row of X_k,
+ // i.e. choose Z(k) such that
+ // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
+ RealScalar beta;
+ m_cpqr.m_qr.row(k)
+ .tail(cols - rank + 1)
+ .makeHouseholderInPlace(m_zCoeffs(k), beta);
+ m_cpqr.m_qr(k, rank - 1) = beta;
+ if (k > 0) {
+ // Apply Z(k) to the first k rows of X_k
+ m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
+ .applyHouseholderOnTheRight(
+ m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
+ &m_temp(0));
+ }
+ if (k != rank - 1) {
+ // Swap X(0:k,k) back to its proper location.
+ m_cpqr.m_qr.col(k).head(k + 1).swap(
+ m_cpqr.m_qr.col(rank - 1).head(k + 1));
+ }
+ }
+ }
+}
+
+template <typename MatrixType>
+template <typename Rhs>
+void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
+ Rhs& rhs) const {
+ const Index cols = this->cols();
+ const Index nrhs = rhs.cols();
+ const Index rank = this->rank();
+ Matrix<typename MatrixType::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
+ for (Index k = 0; k < rank; ++k) {
+ if (k != rank - 1) {
+ rhs.row(k).swap(rhs.row(rank - 1));
+ }
+ rhs.middleRows(rank - 1, cols - rank + 1)
+ .applyHouseholderOnTheLeft(
+ matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
+ &temp(0));
+ if (k != rank - 1) {
+ rhs.row(k).swap(rhs.row(rank - 1));
+ }
+ }
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template <typename _MatrixType>
+template <typename RhsType, typename DstType>
+void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
+ const RhsType& rhs, DstType& dst) const {
+ eigen_assert(rhs.rows() == this->rows());
+
+ const Index rank = this->rank();
+ if (rank == 0) {
+ dst.setZero();
+ return;
+ }
+
+ // Compute c = Q^* * rhs
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is
+ // Q^* = (H_0 H_1 ...)^T
+ typename RhsType::PlainObject c(rhs);
+ c.applyOnTheLeft(
+ householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
+
+ // Solve T z = c(1:rank, :)
+ dst.topRows(rank) = matrixT()
+ .topLeftCorner(rank, rank)
+ .template triangularView<Upper>()
+ .solve(c.topRows(rank));
+
+ const Index cols = this->cols();
+ if (rank < cols) {
+ // Compute y = Z^* * [ z ]
+ // [ 0 ]
+ dst.bottomRows(cols - rank).setZero();
+ applyZAdjointOnTheLeftInPlace(dst);
+ }
+
+ // Undo permutation to get x = P^{-1} * y.
+ dst = colsPermutation() * dst;
+}
+#endif
+
+namespace internal {
+
+template<typename DstXprType, typename MatrixType>
+struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
+{
+ typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
+ typedef Inverse<CodType> SrcXprType;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
+ {
+ dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
+ }
+};
+
+} // end namespace internal
+
+/** \returns the matrix Q as a sequence of householder transformations */
+template <typename MatrixType>
+typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
+CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
+ return m_cpqr.householderQ();
+}
+
+/** \return the complete orthogonal decomposition of \c *this.
+ *
+ * \sa class CompleteOrthogonalDecomposition
+ */
+template <typename Derived>
+const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject>
+MatrixBase<Derived>::completeOrthogonalDecomposition() const {
+ return CompleteOrthogonalDecomposition<PlainObject>(eval());
+}
+
+} // end namespace Eigen
+
+#endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 0b39966..e489bdd 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -15,6 +15,12 @@
namespace internal {
+template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
+ : traits<_MatrixType>
+{
+ enum { Flags = 0 };
+};
+
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
template<typename MatrixType>
@@ -23,7 +29,7 @@
typedef typename MatrixType::PlainObject ReturnType;
};
-}
+} // end namespace internal
/** \ingroup QR_Module
*
@@ -31,19 +37,21 @@
*
* \brief Householder rank-revealing QR decomposition of a matrix with full pivoting
*
- * \param MatrixType the type of the matrix of which we are computing the QR decomposition
+ * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
*
- * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
+ * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b P', \b Q and \b R
* such that
* \f[
- * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
+ * \mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R}
* \f]
- * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
- * upper triangular matrix.
+ * by using Householder transformations. Here, \b P and \b P' are permutation matrices, \b Q a unitary matrix
+ * and \b R an upper triangular matrix.
*
* 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 ColPivHouseholderQR.
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::fullPivHouseholderQr()
*/
template<typename _MatrixType> class FullPivHouseholderQR
@@ -54,21 +62,22 @@
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
+ // FIXME should be int
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
- typedef Matrix<Index, 1,
+ typedef Matrix<StorageIndex, 1,
EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
+ typedef typename MatrixType::PlainObject PlainObject;
/** \brief Default Constructor.
*
@@ -113,7 +122,8 @@
*
* \sa compute()
*/
- FullPivHouseholderQR(const MatrixType& matrix)
+ template<typename InputType>
+ explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
@@ -123,7 +133,27 @@
m_isInitialized(false),
m_usePrescribedThreshold(false)
{
- compute(matrix);
+ compute(matrix.derived());
+ }
+
+ /** \brief Constructs a QR factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
+ *
+ * \sa FullPivHouseholderQR(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit FullPivHouseholderQR(EigenBase<InputType>& matrix)
+ : m_qr(matrix.derived()),
+ m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
+ m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
+ m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
+ m_cols_permutation(matrix.cols()),
+ m_temp(matrix.cols()),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false)
+ {
+ computeInPlace();
}
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
@@ -134,9 +164,6 @@
* \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A,
* and an arbitrary solution otherwise.
*
- * \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
@@ -145,11 +172,11 @@
* Output: \verbinclude FullPivHouseholderQR_solve.out
*/
template<typename Rhs>
- inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
+ inline const Solve<FullPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
- return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
+ return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
/** \returns Expression object representing the matrix Q
@@ -164,7 +191,8 @@
return m_qr;
}
- FullPivHouseholderQR& compute(const MatrixType& matrix);
+ template<typename InputType>
+ FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
/** \returns a const reference to the column permutation matrix */
const PermutationType& colsPermutation() const
@@ -280,13 +308,11 @@
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
- */ inline const
- internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
- inverse() const
+ */
+ inline const Inverse<FullPivHouseholderQR> inverse() const
{
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
- return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
- (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
+ return Inverse<FullPivHouseholderQR>(*this);
}
inline Index rows() const { return m_qr.rows(); }
@@ -366,6 +392,12 @@
* diagonal coefficient of U.
*/
RealScalar maxPivot() const { return m_maxpivot; }
+
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ template<typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl(const RhsType &rhs, DstType &dst) const;
+ #endif
protected:
@@ -374,6 +406,8 @@
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
+ void computeInPlace();
+
MatrixType m_qr;
HCoeffsType m_hCoeffs;
IntDiagSizeVectorType m_rows_transpositions;
@@ -411,16 +445,25 @@
* \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
*/
template<typename MatrixType>
-FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
+template<typename InputType>
+FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
+{
+ m_qr = matrix.derived();
+ computeInPlace();
+ return *this;
+}
+
+template<typename MatrixType>
+void FullPivHouseholderQR<MatrixType>::computeInPlace()
{
check_template_parameters();
-
+
using std::abs;
- Index rows = matrix.rows();
- Index cols = matrix.cols();
+ Index rows = m_qr.rows();
+ Index cols = m_qr.cols();
Index size = (std::min)(rows,cols);
- m_qr = matrix;
+
m_hCoeffs.resize(size);
m_temp.resize(cols);
@@ -439,13 +482,15 @@
for (Index k = 0; k < size; ++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;
- biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
- .cwiseAbs()
- .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+ Score score = m_qr.bottomRightCorner(rows-k, cols-k)
+ .unaryExpr(Scoring())
+ .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
row_of_biggest_in_corner += k;
col_of_biggest_in_corner += k;
+ RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
if(k==0) biggest = biggest_in_corner;
// if the corner is negligible, then we have less than full rank, and we can finish early
@@ -489,50 +534,55 @@
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
-
- return *this;
}
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
- : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename _MatrixType>
+template<typename RhsType, typename DstType>
+void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
- EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
+ eigen_assert(rhs.rows() == rows());
+ const Index l_rank = rank();
- template<typename Dest> void evalTo(Dest& dst) const
+ // FIXME introduce nonzeroPivots() and use it here. and more generally,
+ // make the same improvements in this dec as in FullPivLU.
+ if(l_rank==0)
{
- const Index rows = dec().rows(), cols = dec().cols();
- eigen_assert(rhs().rows() == rows);
+ dst.setZero();
+ return;
+ }
- // FIXME introduce nonzeroPivots() and use it here. and more generally,
- // make the same improvements in this dec as in FullPivLU.
- if(dec().rank()==0)
- {
- dst.setZero();
- return;
- }
+ typename RhsType::PlainObject c(rhs);
- typename Rhs::PlainObject c(rhs());
+ Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
+ for (Index k = 0; k < l_rank; ++k)
+ {
+ Index remainingSize = rows()-k;
+ c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
+ c.bottomRightCorner(remainingSize, rhs.cols())
+ .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
+ m_hCoeffs.coeff(k), &temp.coeffRef(0));
+ }
- Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
- for (Index k = 0; k < dec().rank(); ++k)
- {
- Index remainingSize = rows-k;
- c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
- c.bottomRightCorner(remainingSize, rhs().cols())
- .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
- dec().hCoeffs().coeff(k), &temp.coeffRef(0));
- }
+ m_qr.topLeftCorner(l_rank, l_rank)
+ .template triangularView<Upper>()
+ .solveInPlace(c.topRows(l_rank));
- dec().matrixQR()
- .topLeftCorner(dec().rank(), dec().rank())
- .template triangularView<Upper>()
- .solveInPlace(c.topRows(dec().rank()));
+ for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
+ for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
+}
+#endif
- for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
- for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
+namespace internal {
+
+template<typename DstXprType, typename MatrixType>
+struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
+{
+ typedef FullPivHouseholderQR<MatrixType> QrType;
+ typedef Inverse<QrType> SrcXprType;
+ static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
+ {
+ dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
}
};
@@ -546,7 +596,6 @@
: public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
{
public:
- typedef typename MatrixType::Index Index;
typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
@@ -558,7 +607,7 @@
: m_qr(qr),
m_hCoeffs(hCoeffs),
m_rowsTranspositions(rowsTranspositions)
- {}
+ {}
template <typename ResultType>
void evalTo(ResultType& result) const
@@ -588,8 +637,8 @@
}
}
- Index rows() const { return m_qr.rows(); }
- Index cols() const { return m_qr.rows(); }
+ Index rows() const { return m_qr.rows(); }
+ Index cols() const { return m_qr.rows(); }
protected:
typename MatrixType::Nested m_qr;
@@ -597,6 +646,11 @@
typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
};
+// template<typename MatrixType>
+// struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
+// : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
+// {};
+
} // end namespace internal
template<typename MatrixType>
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 343a664..3513d99 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -21,7 +21,7 @@
*
* \brief Householder QR decomposition of a matrix
*
- * \param MatrixType the type of the matrix of which we are computing the QR decomposition
+ * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R
* such that
@@ -37,6 +37,8 @@
* This Householder QR decomposition is faster, but less numerically stable and less feature-full than
* FullPivHouseholderQR or ColPivHouseholderQR.
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::householderQr()
*/
template<typename _MatrixType> class HouseholderQR
@@ -47,13 +49,13 @@
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
+ // FIXME should be int
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
@@ -91,13 +93,32 @@
*
* \sa compute()
*/
- HouseholderQR(const MatrixType& matrix)
+ template<typename InputType>
+ explicit HouseholderQR(const EigenBase<InputType>& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
m_temp(matrix.cols()),
m_isInitialized(false)
{
- compute(matrix);
+ compute(matrix.derived());
+ }
+
+
+ /** \brief Constructs a QR factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
+ * \c MatrixType is a Eigen::Ref.
+ *
+ * \sa HouseholderQR(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit HouseholderQR(EigenBase<InputType>& matrix)
+ : m_qr(matrix.derived()),
+ m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
+ m_temp(matrix.cols()),
+ m_isInitialized(false)
+ {
+ computeInPlace();
}
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
@@ -107,9 +128,6 @@
*
* \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
@@ -118,11 +136,11 @@
* Output: \verbinclude HouseholderQR_solve.out
*/
template<typename Rhs>
- inline const internal::solve_retval<HouseholderQR, Rhs>
+ inline const Solve<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
- return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
+ return Solve<HouseholderQR, Rhs>(*this, b.derived());
}
/** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations.
@@ -148,7 +166,12 @@
return m_qr;
}
- HouseholderQR& compute(const MatrixType& matrix);
+ template<typename InputType>
+ HouseholderQR& compute(const EigenBase<InputType>& matrix) {
+ m_qr = matrix.derived();
+ computeInPlace();
+ return *this;
+ }
/** \returns the absolute value of the determinant of the matrix of which
* *this is the QR decomposition. It has only linear complexity
@@ -187,6 +210,12 @@
* For advanced uses only.
*/
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ template<typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl(const RhsType &rhs, DstType &dst) const;
+ #endif
protected:
@@ -194,6 +223,8 @@
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
+
+ void computeInPlace();
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@@ -224,7 +255,6 @@
template<typename MatrixQR, typename HCoeffs>
void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
{
- typedef typename MatrixQR::Index Index;
typedef typename MatrixQR::Scalar Scalar;
typedef typename MatrixQR::RealScalar RealScalar;
Index rows = mat.rows();
@@ -263,11 +293,9 @@
struct householder_qr_inplace_blocked
{
// This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
- static void run(MatrixQR& mat, HCoeffs& hCoeffs,
- typename MatrixQR::Index maxBlockSize=32,
+ static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
typename MatrixQR::Scalar* tempData = 0)
{
- typedef typename MatrixQR::Index Index;
typedef typename MatrixQR::Scalar Scalar;
typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
@@ -289,8 +317,8 @@
for (k = 0; k < size; k += blockSize)
{
Index bs = (std::min)(size-k,blockSize); // actual size of the block
- Index tcols = cols - k - bs; // trailing columns
- Index brows = rows-k; // rows of the block
+ Index tcols = cols - k - bs; // trailing columns
+ Index brows = rows-k; // rows of the block
// partition the matrix:
// A00 | A01 | A02
@@ -308,44 +336,39 @@
if(tcols)
{
BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
- apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
+ apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
}
}
}
};
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
- : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
-{
- EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- const Index rows = dec().rows(), cols = dec().cols();
- const Index rank = (std::min)(rows, cols);
- eigen_assert(rhs().rows() == rows);
-
- typename Rhs::PlainObject c(rhs());
-
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
- c.applyOnTheLeft(householderSequence(
- dec().matrixQR().leftCols(rank),
- dec().hCoeffs().head(rank)).transpose()
- );
-
- dec().matrixQR()
- .topLeftCorner(rank, rank)
- .template triangularView<Upper>()
- .solveInPlace(c.topRows(rank));
-
- dst.topRows(rank) = c.topRows(rank);
- dst.bottomRows(cols-rank).setZero();
- }
-};
-
} // end namespace internal
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename _MatrixType>
+template<typename RhsType, typename DstType>
+void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
+{
+ const Index rank = (std::min)(rows(), cols());
+ eigen_assert(rhs.rows() == rows());
+
+ typename RhsType::PlainObject c(rhs);
+
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+ c.applyOnTheLeft(householderSequence(
+ m_qr.leftCols(rank),
+ m_hCoeffs.head(rank)).transpose()
+ );
+
+ m_qr.topLeftCorner(rank, rank)
+ .template triangularView<Upper>()
+ .solveInPlace(c.topRows(rank));
+
+ dst.topRows(rank) = c.topRows(rank);
+ dst.bottomRows(cols()-rank).setZero();
+}
+#endif
+
/** Performs the QR factorization of the given matrix \a matrix. The result of
* the factorization is stored into \c *this, and a reference to \c *this
* is returned.
@@ -353,15 +376,14 @@
* \sa class HouseholderQR, HouseholderQR(const MatrixType&)
*/
template<typename MatrixType>
-HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
+void HouseholderQR<MatrixType>::computeInPlace()
{
check_template_parameters();
- Index rows = matrix.rows();
- Index cols = matrix.cols();
+ Index rows = m_qr.rows();
+ Index cols = m_qr.cols();
Index size = (std::min)(rows,cols);
- m_qr = matrix;
m_hCoeffs.resize(size);
m_temp.resize(cols);
@@ -369,7 +391,6 @@
internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
m_isInitialized = true;
- return *this;
}
/** \return the Householder QR decomposition of \c *this.
diff --git a/Eigen/src/QR/HouseholderQR_MKL.h b/Eigen/src/QR/HouseholderQR_LAPACKE.h
similarity index 77%
rename from Eigen/src/QR/HouseholderQR_MKL.h
rename to Eigen/src/QR/HouseholderQR_LAPACKE.h
index b80f1b4..1dc7d53 100644
--- a/Eigen/src/QR/HouseholderQR_MKL.h
+++ b/Eigen/src/QR/HouseholderQR_LAPACKE.h
@@ -25,47 +25,44 @@
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************
- * Content : Eigen bindings to Intel(R) MKL
+ * Content : Eigen bindings to LAPACKe
* Householder QR decomposition of a matrix w/o pivoting based on
* LAPACKE_?geqrf function.
********************************************************************************
*/
-#ifndef EIGEN_QR_MKL_H
-#define EIGEN_QR_MKL_H
-
-#include "../Core/util/MKL_support.h"
+#ifndef EIGEN_QR_LAPACKE_H
+#define EIGEN_QR_LAPACKE_H
namespace Eigen {
- namespace internal {
+namespace internal {
- /** \internal Specialization for the data types supported by MKL */
+/** \internal Specialization for the data types supported by LAPACKe */
-#define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \
+#define EIGEN_LAPACKE_QR_NOPIV(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX) \
template<typename MatrixQR, typename HCoeffs> \
struct householder_qr_inplace_blocked<MatrixQR, HCoeffs, EIGTYPE, true> \
{ \
- static void run(MatrixQR& mat, HCoeffs& hCoeffs, \
- typename MatrixQR::Index = 32, \
+ static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index = 32, \
typename MatrixQR::Scalar* = 0) \
{ \
lapack_int m = (lapack_int) mat.rows(); \
lapack_int n = (lapack_int) mat.cols(); \
lapack_int lda = (lapack_int) mat.outerStride(); \
lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
- LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \
+ LAPACKE_##LAPACKE_PREFIX##geqrf( matrix_order, m, n, (LAPACKE_TYPE*)mat.data(), lda, (LAPACKE_TYPE*)hCoeffs.data()); \
hCoeffs.adjointInPlace(); \
} \
};
-EIGEN_MKL_QR_NOPIV(double, double, d)
-EIGEN_MKL_QR_NOPIV(float, float, s)
-EIGEN_MKL_QR_NOPIV(dcomplex, MKL_Complex16, z)
-EIGEN_MKL_QR_NOPIV(scomplex, MKL_Complex8, c)
+EIGEN_LAPACKE_QR_NOPIV(double, double, d)
+EIGEN_LAPACKE_QR_NOPIV(float, float, s)
+EIGEN_LAPACKE_QR_NOPIV(dcomplex, lapack_complex_double, z)
+EIGEN_LAPACKE_QR_NOPIV(scomplex, lapack_complex_float, c)
} // end namespace internal
} // end namespace Eigen
-#endif // EIGEN_QR_MKL_H
+#endif // EIGEN_QR_LAPACKE_H