devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add block Householder right-side application for HouseholderSequence

Browse Source 
libeigen/eigen!2342

Closes #3057

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-03-27 19:56:08 -07:00
parent 79d7d280a5
commit cf508c096b
6 changed files with 448 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

45
Eigen/src/Householder/BlockHouseholder.h
View File

@@ -36,14 +36,8 @@ void make_block_householder_triangular_factor(TriangularFactorType& triFactor, c
triFactor.row(i).tail(rt).noalias() = -hCoeffs(i) * vectors.col(i).tail(rs).adjoint() *
vectors.bottomRightCorner(rs, rt).template triangularView<UnitLower>();
// FIXME: use the following line with .noalias() once triangular product supports in-place operation.
// triFactor.row(i).tail(rt) = triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt,rt).template
// triangularView<Upper>();
for (Index j = nbVecs - 1; j > i; --j) {
typename TriangularFactorType::Scalar z = triFactor(i, j);
triFactor(i, j) = z * triFactor(j, j);
if (nbVecs - j - 1 > 0) triFactor.row(i).tail(nbVecs - j - 1) += z * triFactor.row(j).tail(nbVecs - j - 1);
}
triFactor.row(i).tail(rt) =
(triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt, rt).template triangularView<Upper>()).eval();
}
triFactor(i, i) = hCoeffs(i);
}
@@ -71,14 +65,43 @@ void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vec
(VectorsType::MaxColsAtCompileTime == 1 && MatrixType::MaxColsAtCompileTime != 1) ? RowMajor : ColMajor,
VectorsType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime>
tmp = V.adjoint() * mat;
// FIXME: add .noalias() once triangular product supports in-place operation.
if (forward)
tmp = T.template triangularView<Upper>() * tmp;
tmp = (T.template triangularView<Upper>() * tmp).eval();
else
tmp = T.template triangularView<Upper>().adjoint() * tmp;
tmp = (T.template triangularView<Upper>().adjoint() * tmp).eval();
mat.noalias() -= V * tmp;
}
/** \internal
* if forward then perform mat = mat * H0 * H1 * H2
* otherwise perform mat = mat * H2 * H1 * H0
*/
template <typename MatrixType, typename VectorsType, typename CoeffsType>
void apply_block_householder_on_the_right(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs,
bool forward) {
enum { TFactorSize = VectorsType::ColsAtCompileTime };
Index nbVecs = vectors.cols();
Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, RowMajor> T(nbVecs, nbVecs);
if (forward)
make_block_householder_triangular_factor(T, vectors, hCoeffs);
else
make_block_householder_triangular_factor(T, vectors, hCoeffs.conjugate());
const TriangularView<const VectorsType, UnitLower> V(vectors);
// A -= (A * V) * T * V^* (forward)
// A -= (A * V) * T^* * V^* (backward)
Matrix<typename MatrixType::Scalar, MatrixType::RowsAtCompileTime, VectorsType::ColsAtCompileTime,
(MatrixType::MaxRowsAtCompileTime == 1 && VectorsType::MaxColsAtCompileTime != 1) ? ColMajor : RowMajor,
MatrixType::MaxRowsAtCompileTime, VectorsType::MaxColsAtCompileTime>
tmp = mat * V;
if (forward)
tmp = (tmp * T.template triangularView<Upper>()).eval();
else
tmp = (tmp * T.template triangularView<Upper>().adjoint()).eval();
mat.noalias() -= tmp * V.adjoint();
}
} // end namespace internal
} // end namespace Eigen

50
Eigen/src/Householder/HouseholderSequence.h
View File

@@ -281,10 +281,7 @@ class HouseholderSequence : public EigenBase<HouseholderSequence<VectorsType, Co
for (Index k = 0; k < cols() - vecs; ++k) dst.col(k).tail(rows() - k - 1).setZero();
} else if (m_length > BlockSize) {
dst.setIdentity(rows(), rows());
if (m_reverse)
applyThisOnTheLeft(dst, workspace, true);
else
applyThisOnTheLeft(dst, workspace, true);
applyThisOnTheLeft(dst, workspace, true);
} else {
dst.setIdentity(rows(), rows());
for (Index k = vecs - 1; k >= 0; --k) {
@@ -309,14 +306,49 @@ class HouseholderSequence : public EigenBase<HouseholderSequence<VectorsType, Co
/** \internal */
template <typename Dest, typename Workspace>
inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const {
workspace.resize(dst.rows());
for (Index k = 0; k < m_length; ++k) {
Index actual_k = m_reverse ? m_length - k - 1 : k;
dst.rightCols(rows() - m_shift - actual_k)
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
// Use the block path when the reflectors are long enough for GEMM to outperform GEMV.
// The threshold on rows() (the reflector length) is higher than for the left-side path because
// the right-side block application has more overhead from the tmp = mat * V product layout.
if (m_length >= BlockSize && rows() - m_shift >= 4 * BlockSize) {
applyBlockOnTheRight(dst);
} else {
workspace.resize(dst.rows());
for (Index k = 0; k < m_length; ++k) {
Index actual_k = m_reverse ? m_length - k - 1 : k;
dst.rightCols(rows() - m_shift - actual_k)
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
}
}
}
private:
/** \internal Block Householder application on the right, kept out-of-line
* to avoid template bloat pessimizing the scalar path above. */
template <typename Dest>
EIGEN_DONT_INLINE void applyBlockOnTheRight(Dest& dst) const {
// Make sure we have at least 2 useful blocks, otherwise it is point-less:
Index blockSize = m_length < Index(2 * BlockSize) ? (m_length + 1) / 2 : Index(BlockSize);
for (Index i = 0; i < m_length; i += blockSize) {
// Right-side application processes blocks in opposite order to left-side:
// forward (m_reverse=false): first block first; reversed: last block first.
Index end = m_reverse ? m_length - i : (std::min)(m_length, i + blockSize);
Index k = m_reverse ? (std::max)(Index(0), end - blockSize) : i;
Index bs = end - k;
Index start = k + m_shift;
typedef Block<internal::remove_all_t<VectorsType>, Dynamic, Dynamic> SubVectorsType;
SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side == OnTheRight ? k : start,
Side == OnTheRight ? start : k, Side == OnTheRight ? bs : m_vectors.rows() - start,
Side == OnTheRight ? m_vectors.cols() - start : bs);
std::conditional_t<Side == OnTheRight, Transpose<SubVectorsType>, SubVectorsType&> sub_vecs(sub_vecs1);
Index dstCols = rows() - m_shift - k;
auto sub_dst = dst.rightCols(dstCols);
internal::apply_block_householder_on_the_right(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
}
}
public:
/** \internal */
template <typename Dest>
inline void applyThisOnTheLeft(Dest& dst, bool inputIsIdentity = false) const {