lapackers

[James Demmel <demmel@cs.berkeley.edu>]

Julien, thanks for the response.

Your idea of keeping slarfp for QR alone, and dropping it in eig, SVD, etc.,
is interesting. But I wonder about the case of the SVD of a tall skinny 
matrix,
which starts by doing QR to reduce to a smaller matrix for bidiagonal 
reduction.
It seems that we'd want the potential speed benefit of the new QR, if 
the original
tall-skinny matrix also has a narrow profile (this could make the 
overall SVD
much cheaper). I also wonder if we can easily separate which calls should
be to the old slarfg and which to slarfp, it depends on the call trees 
of all the
routines.

So there is a set of design choices. We will be having a conference call 
with the
Mathworks in about 30 minutes, and maybe get some feedback.

Jim


Julien Langou wrote:
> Hi Jim,
>
> Thanks for the very clear update. I do not have much strong feelings.
>
> - I am in favor of "as simple as" possible codes. I did not consider 
> the patch 1 necessarilly simple actually! But if patch 3 and 4 are 
> even more complex ... having an extra routine to compute the norm 
> diferently might be cleaner at the end.
>
> - I am not worried about a factor of 4x, 5x in loss of orthogonality 
> additional to the O(n) anyhow.
>
> - We could also
>   - restore xlarfg in all the routines we do not know why we are using 
> xlarfp ... (eg: eig, svd, etc.) then for QR factorization I do not 
> know maybe having the Alan Edleman QR factorization with xlarfp.
>
> best wishes, --julien.
>
>
>
>
> On Thu, 6 May 2010, James Demmel wrote:
>
> > I attach some interesting data Guillaume has collected from our work
> > on debugging slarfp. We are trying to come to a conclusion before
> > reporting progress to the Mathworks tomorrow. But depending on
> > our ultimate goal (cost vs accuracy), the data comparing various 
> > possible
> > patches for slarfp is not conclusive, so we'd like your opinion.
> >
> > Briefly, there are several patches of different costs and worst-case
> > accuracies that are much better than the buggy version, and probably
> > "good enough". But all are less accurate than the original slarfg, by
> > about 4x or 5x (eyeball metric from the plots; more exact data is
> > available).
> > Whether this lower accuracy is worth worrying about is a judgment
> > call about which I'd like your opinion.
> >
> > To remind you of notation, the input is the vector [alpha;x] and the
> > output is the orthogonal matrix Q = I - tau*[1;v]*[1,v'] and beta >=0
> > such that Q*[alpha;x] = [beta;0].
> >
> > The attached pdf has 3 pages, showing histograms of errors over
> > certain "torture test" suites. The error shown is loss of orthogonality,
> > measured as norm(Q'*Q-I), computed in double. (There is also
> > data for the residual norm(Q*[alpha;x] - [beta;0])/beta, but
> > this is enough data for now). More precisely, we histogram
> >   log10(error/macheps)
> > i.e. how many decimal digits we lose over the "best" error of
> > macheps = 2^(-24). Actually, a large number of test cases had
> > smaller errors than this, and are not shown, but counted as
> > "not shown" in the caption to each figure.
> >
> > There are 3 pages of plots corresponding to 3 kinds of tests.
> > The input vector is of the form [alpha;x], and the dimension
> > of x can be 2, 100 and 1000. When the dimension is 2,
> > each component alpha and x(i) is chosen from a large set
> > of floating point numbers spanning the range from 0,
> > through the denorms, and from underflow to overflow threshholds,
> > with various mantissas and signs. When the dimension
> > of x is 100 or 1000, alpha and x(1) are chosen similarly generally,
> > but x(2:end) are chosen very tiny compared to x(1) to maximize
> > the relative error in computing snrm2(x) (i.e. x(2:end) are all
> > about sqrt(macheps)*x(1)).
> >
> > For each page of tests, 5 plots are shown:
> >   (a) slarfg: original implementation, before our "improvements"
> >   (b) slarfp (3.2, buggy): the version from LAPACK 3.2 that
> >         caused the bug report
> >   (c) slarfp (patch 1): the version that Julien and I worked on
> >   (d) slarfp (patch 2): not a practical patch, rather it uses a bit of
> >             double precision in a subset of the code that we identified
> >             as problematic from our delta-debugging work
> >   (e) slarfp (patch 3): a suggestion of Velvel, where the output tau
> >             is explicitly computed from its definition 2/(1+norm(v,2)^2)
> >             and so costing an additional O(n) work
> >   (f) slarfp (patch 4): a suggestion of Guillaume, that also recomputes
> >             a vector norm, but differently.
> >
> > The main features of these plots are as follows:
> >    The bug in slarfp  (3.2) shows up clearly.
> >    For dimension n = 100 and n=1000, the worst case error in slarfg
> >       and all patched versions of slarfp grows like O(n*macheps) as
> > expected,
> >       since this is the error from snrm2. We have lived with this 
> > situation
> >       (in slarfg) for a long time, and it dominates the differences
> > between the
> >       different patches. We could decide that because of this, the
> > difference
> >       between the patches is too small to matter, and we should pick the
> >       simplest and cheapest (patch 1).
> >    For dimension n=2, the worst case error for slarfg and slarfp 
> > (patch 3)
> >       are about the same, and about 4x smaller than slarfp (patches 1
> > and 2)
> >       and 20x smaller than patch 4. (The 4x and 20x are "eyeball"
> > values, Guillaume
> >       has the exact data.)
> >   For dimension n=100, slarfg is about 5x (eyeball) better than
> >       patches 1, 2 and 3, which are about the same, but only 2x 
> > better than
> >       patch 4.
> >   For dimension n = 1000, the situation is similar, but now patch 4 is
> > about as good
> >      as slarfg.
> >
> > So to repeat my question: do we care that the simplest patch, and
> > sometimes the
> > fancier ones, are as much as 4x to 5x worse than the original slarfg
> > (epending
> > on dimension)?
> >
> > Jim