Date: Mon, 10 Aug 2009 08:50:16 -0600
From: Pat Quillen
To: Sven Hammarling
Cc: Penny Anderson, Duncan Po, Cleve Moler, Bobby Cheng, lapack
Subject: Re: [Lapack] xLARFP and scaling.

Dear Sven,

Thanks for your note.  I had realized that the intent of xLARFP was to
construct TAU and V such that (I - TAU*V*V')(ALPHA; X')' = (BETA; 0')' with
BETA > 0 and I had read LAWN 203 discussing these changes.  In addition to my
previously proposed change (which constructed bounded TAU but at the price of
needing to store V(1) separately) I have thought of another way to change
xLARFP, which I believe preserves the desired non-negative diagonal behavior.

Essentially, one could maintain exactly the construction of xLARFG, however if
BETA comes out negative, set it to be positive and mark TAU to be negative.
Then, when applying the reflector, if TAU < 0, apply |TAU|*V*V' - I instead of
the usual I - TAU*V*V'.  That is, if xLARFG wants to construct a reflector to
make BETA negative, construct that reflector (call it H) but apply -H instead,
using the sign of TAU as a sentinel to indicate this behavior.  Also, if I
understand correctly, one could do the same thing for the complex versions of
these codes by keying off the real part of TAU.

An annoyance with this approach is that xGER only allows A = A + a*x*y'.   I'd
really like to be able to say A = b*A + a*x*y', where in this case b = -1.
Without such a call, one needs to make an extra pass over the memory to scale
by -1, which is bothersome.

I've attached some code (dlarfp3.f, dlarf2.f) based on LAPACK 3.2.1 which
follows this approach.  Again, all comments and advice are welcome.  We are
getting close to a foundation change freeze (Monday 17 August) here at The
MathWorks, and we'd like to have some resolution to these xLARFP issues.

Thanks.
Pat.