Thanks for the bug report. I have added your bug report to the LAPACK
Errata file. See: http://www.netlib.org/lapack/Errata/ , bug0037. We will
try to have a look at some points.
On Mon, 10 Aug 2009, Pat Quillen wrote:
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.
From: Sven Hammarling [mailto:sven.hammarling@Domain.Removed]
Sent: Wednesday, August 05, 2009 5:16 AM
To: Pat Quillen
Cc: lapack@Domain.Removed; Bobby Cheng; Duncan Po; Penny Anderson; Cleve
Subject: Re: [Lapack] xLARFP and scaling.
Thank you for pointing out this behaviour. Although I was not involved
in writing this routine, I don't imagine that this behaviour was
intended. I shall try to take a look at xLARFP in the next few days
see if I can suggest anything to improve this.
As I am sure that you realise, xLARFP was written to produce a
non-negative BETA, so that in a factorization such as the QR, R would
have non-negative diagonals. (Although I notice that the comments in
xGEQRF and xGEQR2 do not actually say that the diagonals are real and