[eigen] makeHouseholder()

[Peter <list@xxxxxxxxxxxxxxxxx>]

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);
  }
}