Re: [eigen] notes on Householder and Givens

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

hm, ok, now i can see that Golub & van Loan's code does handle
correctly the case x0=0.

What I still disagree with is their " if(sigma==0) ". Why does this
come into play ? Seems redundant to me. The code I propose has only 1
branching.

Cheers,
Benoit

2009/2/4 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> Something just occured to me.
>
> If x_0 is "zero", that is muchSmallerThan(norm(x)), then there is
> nothing wrong with this formula:
>
> v = x + norm(x) e_0
>
> Indeed, no cancellation can occur since x0 is much smaller than
> norm(x), and after that v0 is guaranteed to be approximately equal to
> norm(x).
>
> So, i'm sorry to not have thought of that earlier, but it seems that
> no pivoting is needed: i think that the following
> ei_compute_householder will work for all nonzero vectors x even if x0
> = 0:
>
> void ei_compute_householder(const Vector& x of size n, Householder *h)
> {
>  // compute natural householder vector v
>  Vector v = x;
>  if(ei_isMuchSmallerThan(x0, norm(x))
>  {
>   v0 += norm(x);
>  }
>  else
>  {
>   v0 += ( x_0 / abs(x_0) ) norm(x);
>  }
>
>  Vector v_starting_with_1 = v/v0;
>
>  h->essential = v_starting_with_1.end(n-1);
>  h->beta = 2 / v_starting_with_1.squaredNorm();
> };
>
> Notice that this is still very different from Golub's code since our
> code takes good care of the case x0=0 while they don't.
>
> Cheers,
> Benoit
>