Re: [eigen] makeHouseholder()

[Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx>]

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

--
 Dipl. Inf., Dipl. Math. Christoph Hertzberg

 Universität Bremen
 FB 3 - Mathematik und Informatik
 AG Robotik
 Robert-Hooke-Straße 1
 28359 Bremen, Germany

 Zentrale: +49 421 178 45-6611

 Besuchsadresse der Nebengeschäftsstelle:
 Robert-Hooke-Straße 5
 28359 Bremen, Germany

 Tel.:    +49 421 178 45-4021
 Empfang: +49 421 178 45-6600
 Fax:     +49 421 178 45-4150
 E-Mail:  chtz@xxxxxxxxxxxxxxxxxxxxxxxx

 Weitere Informationen: http://www.informatik.uni-bremen.de/robotik