lapackers

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

Regarding underflow, does it cause loss of relative accuracy even when 
the final computed singular value
is far from underflow? The scaling strategy that we have used in 
analogous routines is
(1) scale the matrix to have as large a norm as possible that overflow 
is not a consideration
(eg norm =  sqrt(overflow), if we multiply 2 matrix entries somewhere in 
the code)
(2) Use a "mixed" relative/absolute stopping criterion that guarantees 
tiny relative
error for large enough numbers (say sqrt(underflow)) but uses an 
absolute tolerance
(say macheps*sqrt(underflow)) for tinier ones.
(3) scale the result back, and warn the user that the code "only" guarantees
high relative accuracy for tiny singular values that are at least 
underflow*norm(A).
This would be more than satisfactory. If the code multiplies 3 matrix 
entries,
the sqrt() above might change to cube-root, and underflow*norm(A) to
underflow^(2/3)*norm(A), but this is not a big deal.

Jim

Paul Willems wrote:
> Hello Prof. Demmel,
>
> well, I am not graduated yet, there is the small matter of the oral exam
> still missing ...
>
> You can download the thesis from
>
> http://www-ai.math.uni-wuppertal.de/~willems/thesis
>
> This is the version I handed in, not yet the published one. Be warned:
> It is too long.
>
> Indeed, we found some ways to improve the tridiagonal MRRR and they are
> described (Chap 2).
> Blocks can greatly help concerning robustness, also for the tridiagonal
> case (Chap 4),
> but I cannot really evaluate their impact on efficiency until we have an
> optimized code.
> The one I am working on right now will definitely use blocks.
>
> There is one issue that might be of special interest to lapackers in
> general (i.e., "the other folks cc-ed on this list"):
> There are cases where as-yet unregarded *underflow* in the
> qd-transformations can lose relative accuracy
> (basically: at some points we multiply 3 or more numbers). However, I do
> not see an easy scaling to
> avoid it always. This affects the blocked- and non-blocked versions.
> It does not happen often, but when it does it can completely ruin your day.
>     The solution I will probably use is based on IEEE-exception flags:
> Run fast loop first, check flags, rerun
> safe version if under- or overflow or invalid ops occured. The fast loop
> is then actually faster than before
> (no pivmin-check for the pivots).
>     Now, lacking Fortran2003 I plan to make a detour over C99-routines
> to get to the flags.
> This might spoil LAPACK-integrability somewhat (?).
>
> Best regards,
> Paul.
>
>
>
> James Demmel schrieb:
>   
> > Congratulations on graduating! Where can we see a copy of your thesis?
> >
> > Thanks for the update on MRRR for the bidiagonal SVD (bSVD). Does your
> > thesis describe the improvements you mention for both the tridiagonal
> > and bidiagonal cases? Are you saying that the block diagonal LDL'
> > factorization that you need for stability in the bSVD-MRRR is also useful
> > just for the usual tridiagonal case?
> >
> > Let me take the liberty of attaching your Monterrey slide talks for the
> > the other folks cc-ed on this list.
> >
> > Jim
> >
> >
> >
> > Paul Willems wrote:
> >     
> > > Hello Prof. Demmel,
> > >
> > > thanks for your email and interest. I hope you had nice holidays and a
> > > good transition into the new year.
> > >
> > > Since Monterey two things have happened:
> > > 1) Some parts of the codebase were redesigned to improve robustness,
> > > so that
> > > at the moment we have codes for the tridiagonal and bidiagonal case that
> > > are working quite satisfactory. The problems that were still there in
> > > november are resolved
> > > now.
> > > 2) I finished my thesis (finally) and handed it in a couple of days
> > > before christmas.
> > >
> > > Concerning the code & LAPACK:
> > > Well, some time ago I decided to do my  prototypes in C++, to have more
> > > flexibility in experimenting with different subcomponents, and also for
> > > easier code reuse.
> > > Furthermore the code is heavily instrumented for statistics, logging &
> > > debugging.
> > > The downside is that we do not have an optimized implementation, just a
> > > prototype.
> > >
> > > Thus it would have to be rewritten in Fortran from scratch (although now
> > > we do know how,
> > > so no experimentation but just straight coding). I would love to get the
> > > final code into LAPACK,
> > > of course, but I cannot foresee if I will have the time to do it myself,
> > > since I don't know how long
> > > I will be staying in Wuppertal (just started looking at the job market).
> > >     At the moment I am paid within a research project for large scalable
> > > eigensolvers (ELPA: http://elpa.rzg.mpg.de).
> > > For this we first need our improvements to MRRR in (parallel) optimized
> > > form---this might become a conteder for Sca/LAPACK-MRRR,
> > > but the SVD code has lower priority right now.
> > > [Otoh: Once MRRR is done, many parts can be reused for bSVD, so this
> > > *should* then be not *that* much work]
> > >
> > > I will keep you updated.
> > >
> > > Best regards,
> > > Paul.
> > >
> > > P.S.: Just for interest: Do you already have a concrete schedule for the
> > > next release of LAPACK ?
> > >
> > >
> > > James Demmel schrieb:
> > >  
> > >       
> > > > Just wondering how you're doing, and what the status of
> > > > the SVD code is. We're still interested in putting it in
> > > > Sca/LAPACK!
> > > > Regards,
> > > > Jim
> > > >
> > > > Paul Willems wrote:
> > > >    
> > > >         
> > > > > Hello Prof. Demmel,
> > > > >
> > > > > thanks for your interest in my talk, attached you will find the slides
> > > > > of my talk as presented in Monterey.
> > > > > As the slides alone do not convey that much of information, I will
> > > > > give
> > > > > some comments on the contents.
> > > > >
> > > > > The test results shown are from our prototype implementation over a
> > > > > representative testset (including
> > > > > practical as well as synthetic and glued test matrices).
> > > > >
> > > > > Things are looking up (finally), the Block-Factorizations really were
> > > > > the key.
> > > > > However, I don't want to hide from you that there are some problematic
> > > > > cases where the coupling approach does not work well yet (just BFs on
> > > > > TGK does, though).  I studied these in  detail and  it always boils
> > > > > down
> > > > > to one of two problems:
> > > > >   a)  An employed (non-blocked) representation is not relatively
> > > > > robust,
> > > > > although our a priori criteria (element growth)
> > > > >         make us think so.
> > > > >   b)  Underflow issues in dstqds/dqds that were not resolved yet.
> > > > > Both are actually problems in MRRR/STEGR itself, which is why I did
> > > > > not
> > > > > want to elaborate  on  them in my talk.
> > > > > (In fact I have a couple of cases where the current STEGR doesn't look
> > > > > that well either, due to the same problems...)
> > > > > In my opinion these should really be fixed in a next version of MRRR,
> > > > > and also to get the bSVD working.
> > > > > I have a couple of ideas for a) and b) and already talked with
> > > > > B.Parlett
> > > > > about some of them, we will have to see what
> > > > > comes out of this.
> > > > >
> > > > > Please feel free to ask for more information or if anything of the
> > > > > above
> > > > > was unclear.
> > > > > I will keep you posted about (major) new developments.
> > > > >
> > > > > Best regards,
> > > > > Paul.
> > > > >
> > > > >
> > > > >         
> > > > >