Skip to content

Instantly share code, notes, and snippets.

@yixuan

yixuan/Householder.h Secret

Last active April 17, 2017 19:34
Download ZIP

Revisions

  1. yixuan revised this gist Apr 17, 2017. 1 changed file with 125 additions and 9 deletions.
    134 changes: 125 additions & 9 deletions Householder.h
    View file
    Original file line number Original file line Diff line number Diff line change
    @@ -11,7 +11,7 @@
    #ifndef EIGEN_HOUSEHOLDER_H #ifndef EIGEN_HOUSEHOLDER_H
    #define EIGEN_HOUSEHOLDER_H #define EIGEN_HOUSEHOLDER_H


    namespace Eigen { namespace Eigen {


    namespace internal { namespace internal {
    template<int n> struct decrement_size template<int n> struct decrement_size
    @@ -30,7 +30,7 @@ template<int n> struct decrement_size
    * \f$ v^T = [1 essential^T] \f$ * \f$ v^T = [1 essential^T] \f$
    * *
    * The essential part of the vector \c v is stored in *this. * The essential part of the vector \c v is stored in *this.
    * *
    * On output: * On output:
    * \param tau the scaling factor of the Householder transformation * \param tau the scaling factor of the Householder transformation
    * \param beta the result of H * \c *this * \param beta the result of H * \c *this
    @@ -69,10 +69,10 @@ void MatrixBase<Derived>::makeHouseholder(
    { {
    using std::sqrt; using std::sqrt;
    using numext::conj; using numext::conj;

    EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
    VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);

    RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
    Scalar c0 = coeff(0); Scalar c0 = coeff(0);
    const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
    @@ -105,7 +105,7 @@ void MatrixBase<Derived>::makeHouseholder(
    * \param workspace a pointer to working space with at least * \param workspace a pointer to working space with at least
    * this->cols() * essential.size() entries * this->cols() * essential.size() entries
    * *
    * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
    * MatrixBase::applyHouseholderOnTheRight() * MatrixBase::applyHouseholderOnTheRight()
    */ */
    template<typename Derived> template<typename Derived>
    @@ -115,11 +115,70 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
    const Scalar& tau, const Scalar& tau,
    Scalar* workspace) Scalar* workspace)
    { {
    using numext::conj;
    using numext::abs2;

    if(rows() == 1) if(rows() == 1)
    { {
    *this *= Scalar(1)-tau; *this *= Scalar(1)-tau;
    return;
    }

    // Frequently used in RealSchur
    if(EssentialPart::SizeAtCompileTime == 1)
    {
    // H = 1 - tau * v * v^*
    // v = [1 ], vv^* = [1 , _ess_ ], H = [1 - tau , -tau * _ess_ ] = [h00, h01]
    // [ess] [ess, |ess|^2] [-tau * ess, 1 - tau * |ess|^2] [h10, h11]
    const Scalar h00 = Scalar(1) - tau;
    const Scalar h10 = -tau * essential.coeff(0);
    const Scalar h01 = -tau * conj(essential.coeff(0));
    const Scalar h11 = Scalar(1) - tau * abs2(essential.coeff(0));

    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols());
    tmp.noalias() = h00 * this->row(0) + h01 * this->row(1);
    this->row(1).noalias() = h10 * this->row(0) + h11 * this->row(1);
    this->row(0).noalias() = tmp;

    return;
    }

    // Frequently used in RealSchur
    if(EssentialPart::SizeAtCompileTime == 2)
    {
    // H = 1 - tau * v * v^*
    // H * X = X - tau * v * (v^* * X)
    const Scalar ess1 = essential.coeff(0);
    const Scalar ess2 = essential.coeff(1);
    const Scalar essconj1 = conj(ess1);
    const Scalar essconj2 = conj(ess2);

    // If the storage order of this matrix is row-major, we use vectorized operations
    if(Flags & RowMajorBit)
    {
    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols());
    tmp.noalias() = tau * (this->row(0) + essconj1 * this->row(1) + essconj2 * this->row(2));
    this->row(0) -= tmp;
    this->row(1) -= ess1 * tmp;
    this->row(2) -= ess2 * tmp;
    }
    else // Otherwise, we use only one loop to reduce overhead
    {
    const Index ncol = cols();
    for(Index i = 0; i < ncol; ++i)
    {
    const Scalar tmp = tau * (coeff(0, i) + essconj1 * coeff(1, i) + essconj2 * coeff(2, i));
    coeffRef(0, i) -= tmp;
    coeffRef(1, i) -= ess1 * tmp;
    coeffRef(2, i) -= ess2 * tmp;
    }
    }

    return;
    } }
    else if(tau!=Scalar(0))
    // General case
    if(tau!=Scalar(0))
    { {
    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols()); Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
    Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
    @@ -140,9 +199,9 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
    * \param essential the essential part of the vector \c v * \param essential the essential part of the vector \c v
    * \param tau the scaling factor of the Householder transformation * \param tau the scaling factor of the Householder transformation
    * \param workspace a pointer to working space with at least * \param workspace a pointer to working space with at least
    * this->cols() * essential.size() entries * this->rows() * essential.size() entries
    * *
    * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
    * MatrixBase::applyHouseholderOnTheLeft() * MatrixBase::applyHouseholderOnTheLeft()
    */ */
    template<typename Derived> template<typename Derived>
    @@ -152,11 +211,68 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
    const Scalar& tau, const Scalar& tau,
    Scalar* workspace) Scalar* workspace)
    { {
    using numext::conj;
    using numext::abs2;

    if(cols() == 1) if(cols() == 1)
    { {
    *this *= Scalar(1)-tau; *this *= Scalar(1)-tau;
    return;
    } }
    else if(tau!=Scalar(0))
    // Frequently used in RealSchur
    if(EssentialPart::SizeAtCompileTime == 1)
    {
    // H = 1 - tau * v * v^*
    // v = [1 ], vv^* = [1 , _ess_ ], H = [1 - tau , -tau * _ess_ ] = [h00, h01]
    // [ess] [ess, |ess|^2] [-tau * ess, 1 - tau * |ess|^2] [h10, h11]
    const Scalar h00 = Scalar(1) - tau;
    const Scalar h10 = -tau * essential.coeff(0);
    const Scalar h01 = -tau * conj(essential.coeff(0));
    const Scalar h11 = Scalar(1) - tau * abs2(essential.coeff(0));

    Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace, rows());
    tmp.noalias() = h00 * this->col(0) + h10 * this->col(1);
    this->col(1).noalias() = h01 * this->col(0) + h11 * this->col(1);
    this->col(0).noalias() = tmp;

    return;
    }

    if(EssentialPart::SizeAtCompileTime == 2)
    {
    // H = 1 - tau * v * v^*
    // X * H = X - tau * (X * v) * v^*
    const Scalar ess1 = essential.coeff(0);
    const Scalar ess2 = essential.coeff(1);
    const Scalar essconj1 = conj(ess1);
    const Scalar essconj2 = conj(ess2);

    // If the storage order of this matrix is row-major, use only one loop to reduce overhead
    if(Flags & RowMajorBit)
    {
    const Index nrow = rows();
    for(Index i = 0; i < nrow; ++i)
    {
    const Scalar tmp = tau * (coeff(i, 0) + ess1 * coeff(i, 1) + ess2 * coeff(i, 2));
    coeffRef(i, 0) -= tmp;
    coeffRef(i, 1) -= essconj1 * tmp;
    coeffRef(i, 2) -= essconj2 * tmp;
    }
    }
    else // Otherwise, we use vectorized operations
    {
    Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace, rows());
    tmp.noalias() = tau * (this->col(0) + ess1 * this->col(1) + ess2 * this->col(2));
    this->col(0) -= tmp;
    this->col(1) -= essconj1 * tmp;
    this->col(2) -= essconj2 * tmp;
    }

    return;
    }

    if(tau!=Scalar(0))
    { {
    Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows()); Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
    Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
  2. yixuan created this gist Apr 17, 2017.
    172 changes: 172 additions & 0 deletions Householder.h
    View file
    Original file line number Original file line Diff line number Diff line change
    @@ -0,0 +1,172 @@
    // This file is part of Eigen, a lightweight C++ template library
    // for linear algebra.
    //
    // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
    // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
    //
    // 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_HOUSEHOLDER_H
    #define EIGEN_HOUSEHOLDER_H

    namespace Eigen {

    namespace internal {
    template<int n> struct decrement_size
    {
    enum {
    ret = n==Dynamic ? n : n-1
    };
    };
    }

    /** Computes the elementary reflector H such that:
    * \f$ H *this = [ beta 0 ... 0]^T \f$
    * where the transformation H is:
    * \f$ H = I - tau v v^*\f$
    * and the vector v is:
    * \f$ v^T = [1 essential^T] \f$
    *
    * The essential part of the vector \c v is stored in *this.
    *
    * On output:
    * \param tau the scaling factor of the Householder transformation
    * \param beta the result of H * \c *this
    *
    * \sa MatrixBase::makeHouseholder(), MatrixBase::applyHouseholderOnTheLeft(),
    * MatrixBase::applyHouseholderOnTheRight()
    */
    template<typename Derived>
    void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
    {
    VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
    makeHouseholder(essentialPart, tau, beta);
    }

    /** Computes the elementary reflector H such that:
    * \f$ H *this = [ beta 0 ... 0]^T \f$
    * where the transformation H is:
    * \f$ H = I - tau v v^*\f$
    * and the vector v is:
    * \f$ v^T = [1 essential^T] \f$
    *
    * On output:
    * \param essential the essential part of the vector \c v
    * \param tau the scaling factor of the Householder transformation
    * \param beta the result of H * \c *this
    *
    * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(),
    * MatrixBase::applyHouseholderOnTheRight()
    */
    template<typename Derived>
    template<typename EssentialPart>
    void MatrixBase<Derived>::makeHouseholder(
    EssentialPart& essential,
    Scalar& tau,
    RealScalar& beta) const
    {
    using std::sqrt;
    using numext::conj;

    EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
    VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);

    RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
    Scalar c0 = coeff(0);
    const RealScalar tol = (std::numeric_limits<RealScalar>::min)();

    if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol)
    {
    tau = RealScalar(0);
    beta = numext::real(c0);
    essential.setZero();
    }
    else
    {
    beta = sqrt(numext::abs2(c0) + tailSqNorm);
    if (numext::real(c0)>=RealScalar(0))
    beta = -beta;
    essential = tail / (c0 - beta);
    tau = conj((beta - c0) / beta);
    }
    }

    /** Apply the elementary reflector H given by
    * \f$ H = I - tau v v^*\f$
    * with
    * \f$ v^T = [1 essential^T] \f$
    * from the left to a vector or matrix.
    *
    * On input:
    * \param essential the essential part of the vector \c v
    * \param tau the scaling factor of the Householder transformation
    * \param workspace a pointer to working space with at least
    * this->cols() * essential.size() entries
    *
    * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
    * MatrixBase::applyHouseholderOnTheRight()
    */
    template<typename Derived>
    template<typename EssentialPart>
    void MatrixBase<Derived>::applyHouseholderOnTheLeft(
    const EssentialPart& essential,
    const Scalar& tau,
    Scalar* workspace)
    {
    if(rows() == 1)
    {
    *this *= Scalar(1)-tau;
    }
    else if(tau!=Scalar(0))
    {
    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
    Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
    tmp.noalias() = essential.adjoint() * bottom;
    tmp += this->row(0);
    this->row(0) -= tau * tmp;
    bottom.noalias() -= tau * essential * tmp;
    }
    }

    /** Apply the elementary reflector H given by
    * \f$ H = I - tau v v^*\f$
    * with
    * \f$ v^T = [1 essential^T] \f$
    * from the right to a vector or matrix.
    *
    * On input:
    * \param essential the essential part of the vector \c v
    * \param tau the scaling factor of the Householder transformation
    * \param workspace a pointer to working space with at least
    * this->cols() * essential.size() entries
    *
    * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
    * MatrixBase::applyHouseholderOnTheLeft()
    */
    template<typename Derived>
    template<typename EssentialPart>
    void MatrixBase<Derived>::applyHouseholderOnTheRight(
    const EssentialPart& essential,
    const Scalar& tau,
    Scalar* workspace)
    {
    if(cols() == 1)
    {
    *this *= Scalar(1)-tau;
    }
    else if(tau!=Scalar(0))
    {
    Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
    Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
    tmp.noalias() = right * essential.conjugate();
    tmp += this->col(0);
    this->col(0) -= tau * tmp;
    right.noalias() -= tau * tmp * essential.transpose();
    }
    }

    } // end namespace Eigen

    #endif // EIGEN_HOUSEHOLDER_H