aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
authorGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2009-11-08 10:21:26 -0500
committerGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2009-11-08 10:21:26 -0500
commitba7bfe110cf9a2df84b2691dd19f1cfe13d2356c (patch)
treeea678d9b0c02a06b426ecdcd0a72586c3fa2e003 /Eigen
parent68210b03c17915b49655dfa4a13c28cc31a59092 (diff)
port the qr module to ei_solve_xxx.
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/LU/PartialPivLU.h2
-rw-r--r--Eigen/src/QR/ColPivHouseholderQR.h131
-rw-r--r--Eigen/src/QR/FullPivHouseholderQR.h141
-rw-r--r--Eigen/src/QR/HouseholderQR.h80
-rw-r--r--Eigen/src/SVD/SVD.h1
5 files changed, 181 insertions, 174 deletions
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 8f3b7dfc1..eeec3533f 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -198,8 +198,6 @@ PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
compute(matrix);
}
-
-
/** This is the blocked version of ei_fullpivlu_unblocked() */
template<typename Scalar, int StorageOrder>
struct ei_partial_lu_impl
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 05287ff3c..a774fdd73 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -42,17 +42,17 @@
*
* \sa MatrixBase::colPivHouseholderQr()
*/
-template<typename MatrixType> class ColPivHouseholderQR
+template<typename _MatrixType> class ColPivHouseholderQR
{
public:
+ typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
-
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@@ -83,22 +83,27 @@ template<typename MatrixType> class ColPivHouseholderQR
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
* *this is the QR decomposition, if any exists.
*
- * \returns \c true if a solution exists, \c false if no solution exists.
- *
* \param b the right-hand-side of the equation to solve.
*
- * \param result a pointer to the vector/matrix in which to store the solution, if any exists.
- * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
- * If no solution exists, *result is left with undefined coefficients.
+ * \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
+ * \note_about_checking_solutions
+ *
+ * \note_about_arbitrary_choice_of_solution
+ *
* Example: \include ColPivHouseholderQR_solve.cpp
* Output: \verbinclude ColPivHouseholderQR_solve.out
*/
- template<typename OtherDerived, typename ResultType>
- bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+ template<typename Rhs>
+ inline const ei_solve_return_value<ColPivHouseholderQR, Rhs>
+ solve(const MatrixBase<Rhs>& b) const
+ {
+ ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
+ return ei_solve_return_value<ColPivHouseholderQR, Rhs>(*this, b.derived());
+ }
HouseholderSequenceType matrixQ(void) const;
@@ -204,36 +209,24 @@ template<typename MatrixType> class ColPivHouseholderQR
return isInjective() && isSurjective();
}
- /** Computes the inverse of the matrix of which *this is the QR decomposition.
- *
- * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
- *
- * \note If this matrix is not invertible, *result is left with undefined coefficients.
- * Use isInvertible() to first determine whether this matrix is invertible.
- *
- * \sa inverse()
- */
- inline void computeInverse(MatrixType *result) const
- {
- ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
- solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
- }
-
/** \returns the inverse of the matrix of which *this is the QR decomposition.
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
- *
- * \sa computeInverse()
*/
- inline MatrixType inverse() const
+ inline const
+ ei_solve_return_value<ColPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
+ inverse() const
{
- MatrixType result;
- computeInverse(&result);
- return result;
+ ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
+ return ei_solve_return_value<ColPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
+ (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
}
+ inline int rows() const { return m_qr.rows(); }
+ inline int cols() const { return m_qr.cols(); }
+ const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@@ -331,50 +324,56 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
return *this;
}
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-bool ColPivHouseholderQR<MatrixType>::solve(
- const MatrixBase<OtherDerived>& b,
- ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<ColPivHouseholderQR<MatrixType>, Rhs, Dest>
+ : ei_solve_return_value<ColPivHouseholderQR<MatrixType>, Rhs>
{
- ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- result->resize(m_qr.cols(), b.cols());
- if(m_rank==0)
+ void evalTo(Dest& dst) const
{
- if(b.squaredNorm() == RealScalar(0))
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ const ColPivHouseholderQR<MatrixType>& dec = this->m_dec;
+ const Rhs& rhs = this->m_rhs;
+ const int rows = dec.rows(), cols = dec.cols();
+ dst.resize(cols, rhs.cols());
+ ei_assert(rhs.rows() == rows);
+
+ // FIXME introduce nonzeroPivots() and use it here. and more generally,
+ // make the same improvements in this dec as in FullPivLU.
+ if(dec.rank()==0)
{
- result->setZero();
- return true;
+ dst.setZero();
+ return;
}
- else return false;
- }
- const int rows = m_qr.rows();
- ei_assert(b.rows() == rows);
+ typename Rhs::PlainMatrixType c(rhs);
- typename OtherDerived::PlainMatrixType c(b);
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+ c.applyOnTheLeft(makeHouseholderSequence(
+ dec.matrixQR().corner(TopLeft,rows,dec.rank()),
+ dec.hCoeffs().start(dec.rank())).transpose()
+ );
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
- c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose());
-
- if(!isSurjective())
- {
- // is c is in the image of R ?
- RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
- RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
- if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
- return false;
- }
+ if(!dec.isSurjective())
+ {
+ // is c is in the image of R ?
+ RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
+ RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
+ // FIXME brain dead
+ const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
+ if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
+ return;
+ }
- m_qr.corner(TopLeft, m_rank, m_rank)
- .template triangularView<UpperTriangular>()
- .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
+ dec.matrixQR()
+ .corner(TopLeft, dec.rank(), dec.rank())
+ .template triangularView<UpperTriangular>()
+ .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
- for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
- for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
- return true;
-}
+ for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
+ for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
+ }
+};
/** \returns the matrix Q as a sequence of householder transformations */
template<typename MatrixType>
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 07ec343a5..36ec71b95 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -42,17 +42,17 @@
*
* \sa MatrixBase::fullPivHouseholderQr()
*/
-template<typename MatrixType> class FullPivHouseholderQR
+template<typename _MatrixType> class FullPivHouseholderQR
{
public:
+ typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
-
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@@ -78,22 +78,27 @@ template<typename MatrixType> class FullPivHouseholderQR
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
* *this is the QR decomposition, if any exists.
*
- * \returns \c true if a solution exists, \c false if no solution exists.
- *
* \param b the right-hand-side of the equation to solve.
*
- * \param result a pointer to the vector/matrix in which to store the solution, if any exists.
- * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
- * If no solution exists, *result is left with undefined coefficients.
+ * \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
+ * \note_about_checking_solutions
+ *
+ * \note_about_arbitrary_choice_of_solution
+ *
* Example: \include FullPivHouseholderQR_solve.cpp
* Output: \verbinclude FullPivHouseholderQR_solve.out
*/
- template<typename OtherDerived, typename ResultType>
- bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+ template<typename Rhs>
+ inline const ei_solve_return_value<FullPivHouseholderQR, Rhs>
+ solve(const MatrixBase<Rhs>& b) const
+ {
+ ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
+ return ei_solve_return_value<FullPivHouseholderQR, Rhs>(*this, b.derived());
+ }
MatrixQType matrixQ(void) const;
@@ -205,36 +210,23 @@ template<typename MatrixType> class FullPivHouseholderQR
return isInjective() && isSurjective();
}
- /** Computes the inverse of the matrix of which *this is the QR decomposition.
- *
- * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
- *
- * \note If this matrix is not invertible, *result is left with undefined coefficients.
- * Use isInvertible() to first determine whether this matrix is invertible.
- *
- * \sa inverse()
- */
- inline void computeInverse(MatrixType *result) const
- {
- ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
- ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
- solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
- }
-
/** \returns the inverse of the matrix of which *this is the QR decomposition.
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
- *
- * \sa computeInverse()
- */
- inline MatrixType inverse() const
+ */ inline const
+ ei_solve_return_value<FullPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
+ inverse() const
{
- MatrixType result;
- computeInverse(&result);
- return result;
+ ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
+ return ei_solve_return_value<FullPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
+ (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
}
+ inline int rows() const { return m_qr.rows(); }
+ inline int cols() const { return m_qr.cols(); }
+ const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@@ -340,56 +332,59 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
return *this;
}
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-bool FullPivHouseholderQR<MatrixType>::solve(
- const MatrixBase<OtherDerived>& b,
- ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<FullPivHouseholderQR<MatrixType>, Rhs, Dest>
+ : ei_solve_return_value<FullPivHouseholderQR<MatrixType>, Rhs>
{
- ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
- result->resize(m_qr.cols(), b.cols());
- if(m_rank==0)
+ void evalTo(Dest& dst) const
{
- if(b.squaredNorm() == RealScalar(0))
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ const FullPivHouseholderQR<MatrixType>& dec = this->m_dec;
+ const Rhs& rhs = this->m_rhs;
+ const int rows = dec.rows(), cols = dec.cols();
+ dst.resize(cols, rhs.cols());
+ ei_assert(rhs.rows() == rows);
+
+ // FIXME introduce nonzeroPivots() and use it here. and more generally,
+ // make the same improvements in this dec as in FullPivLU.
+ if(dec.rank()==0)
{
- result->setZero();
- return true;
+ dst.setZero();
+ return;
}
- else return false;
- }
- const int rows = m_qr.rows();
- const int cols = b.cols();
- ei_assert(b.rows() == rows);
+ typename Rhs::PlainMatrixType c(rhs);
- typename OtherDerived::PlainMatrixType c(b);
+ Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs.cols());
+ for (int k = 0; k < dec.rank(); ++k)
+ {
+ int remainingSize = rows-k;
+ c.row(k).swap(c.row(dec.rowsTranspositions().coeff(k)));
+ c.corner(BottomRight, remainingSize, rhs.cols())
+ .applyHouseholderOnTheLeft(dec.matrixQR().col(k).end(remainingSize-1),
+ dec.hCoeffs().coeff(k), &temp.coeffRef(0));
+ }
- Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
- for (int k = 0; k < m_rank; ++k)
- {
- int remainingSize = rows-k;
- c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
- c.corner(BottomRight, remainingSize, cols)
- .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
- }
+ if(!dec.isSurjective())
+ {
+ // is c is in the image of R ?
+ RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
+ RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
+ // FIXME brain dead
+ const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
+ if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
+ return;
+ }
+ dec.matrixQR()
+ .corner(TopLeft, dec.rank(), dec.rank())
+ .template triangularView<UpperTriangular>()
+ .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
- if(!isSurjective())
- {
- // is c is in the image of R ?
- RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
- RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
- if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
- return false;
+ for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
+ for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
}
- m_qr.corner(TopLeft, m_rank, m_rank)
- .template triangularView<UpperTriangular>()
- .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
-
- for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
- for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
- return true;
-}
+};
/** \returns the matrix Q */
template<typename MatrixType>
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index a32aa4eaf..6db0411d9 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -46,17 +46,17 @@
*
* \sa MatrixBase::householderQr()
*/
-template<typename MatrixType> class HouseholderQR
+template<typename _MatrixType> class HouseholderQR
{
public:
+ typedef _MatrixType MatrixType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
-
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
@@ -85,19 +85,26 @@ template<typename MatrixType> class HouseholderQR
*
* \param b the right-hand-side of the equation to solve.
*
- * \param result a pointer to the vector/matrix in which to store the solution, if any exists.
- * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
- * If no solution exists, *result is left with undefined coefficients.
+ * \returns a solution.
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
+ * \note_about_checking_solutions
+ *
+ * \note_about_arbitrary_choice_of_solution
+ *
* Example: \include HouseholderQR_solve.cpp
* Output: \verbinclude HouseholderQR_solve.out
*/
- template<typename OtherDerived, typename ResultType>
- void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
-
+ template<typename Rhs>
+ inline const ei_solve_return_value<HouseholderQR, Rhs>
+ solve(const MatrixBase<Rhs>& b) const
+ {
+ ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+ return ei_solve_return_value<HouseholderQR, Rhs>(*this, b.derived());
+ }
+
MatrixQType matrixQ() const;
HouseholderSequenceType matrixQAsHouseholderSequence() const
@@ -145,6 +152,10 @@ template<typename MatrixType> class HouseholderQR
*/
typename MatrixType::RealScalar logAbsDeterminant() const;
+ inline int rows() const { return m_qr.rows(); }
+ inline int cols() const { return m_qr.cols(); }
+ const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@@ -198,31 +209,36 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
return *this;
}
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-void HouseholderQR<MatrixType>::solve(
- const MatrixBase<OtherDerived>& b,
- ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<HouseholderQR<MatrixType>, Rhs, Dest>
+ : ei_solve_return_value<HouseholderQR<MatrixType>, Rhs>
{
- ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
- result->derived().resize(m_qr.cols(), b.cols());
- const int rows = m_qr.rows();
- const int rank = std::min(m_qr.rows(), m_qr.cols());
- ei_assert(b.rows() == rows);
-
- typename OtherDerived::PlainMatrixType c(b);
-
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
- c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose());
-
- m_qr.corner(TopLeft, rank, rank)
- .template triangularView<UpperTriangular>()
- .solveInPlace(c.corner(TopLeft, rank, c.cols()));
-
- result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
- result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
-}
+ void evalTo(Dest& dst) const
+ {
+ const HouseholderQR<MatrixType>& dec = this->m_dec;
+ const Rhs& rhs = this->m_rhs;
+ const int rows = dec.rows(), cols = dec.cols();
+ dst.resize(cols, rhs.cols());
+ const int rank = std::min(rows, cols);
+ ei_assert(rhs.rows() == rows);
+
+ typename Rhs::PlainMatrixType c(rhs);
+
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+ c.applyOnTheLeft(makeHouseholderSequence(
+ dec.matrixQR().corner(TopLeft,rows,rank),
+ dec.hCoeffs().start(rank)).transpose()
+ );
+
+ dec.matrixQR()
+ .corner(TopLeft, rank, rank)
+ .template triangularView<UpperTriangular>()
+ .solveInPlace(c.corner(TopLeft, rank, c.cols()));
+
+ dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols());
+ dst.corner(BottomLeft, cols-rank, c.cols()).setZero();
+ }
+};
/** \returns the matrix Q */
template<typename MatrixType>
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index b43123384..8ca425525 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -86,7 +86,6 @@ template<typename _MatrixType> class SVD
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
- * \note_about_using_kernel_to_study_multiple_solutions
*
* \sa MatrixBase::svd(),
*/