On 2016-09-12 20:03, Peter wrote:
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)
As discussed at the beginning of the week on IRC, I second that change.
I'll be pushing this change (if I don't get any objections within the
next minutes).
We should probably try to make some unit tests which would fail without
this change (preferably not requiring custom types).
Cheers,
Christoph
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);
}
}