[eigen] makeHouseholder() |

[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]

Dear All,
is there are specific reason for
if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol)
or could it be changed into
if(tailSqNorm <= tol && numext::abs2(c0)<=tol)
within makeHouseholder() in Householder.C, [*] see below?
I'm asking as I stumbled over the specific case of c0=1.0+0.0i and a 2x2 matrix,
where tailSqNorm==0.
The current implementation with the use the special case
tau = RealScalar(0);
beta = numext::real(c0);
essential.setZero();
while if(tailSqNorm <= tol && numext::abs2(c0)<=tol)
would select the normal path in that case,
beta = sqrt(numext::abs2(c0) + tailSqNorm);
if (numext::real(c0)>=RealScalar(0))
beta = -beta;
essential = tail / (c0 - beta);
tau = conj((beta - c0) / beta);
which is fine to evaluate.
Due to sign choice of beta, (c0 - beta) will even be further away from zero as beta itself, so numerically
that shouldn't be a problem.
I'm asking because I could pin-point, with the help of Christoph on IRC, this particular if-statement
causing problems for my own custom type.
Best regards,
Peter
[*]
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);
}
}