Dear Prof. Kahan, thanks for providing us with your ICMAX1 trick. It would be
nice to have this trick in ICMAX1 in LAPACK indeed. Cheers, Julien.
On Feb 10, 2014, at 3:28 PM, William Kahan
This responds to what, on 2/10/2014 1:20 PM, James Demmel wrote:
This is a "non-bug" report, but I wanted to report that we get these too.
I note that this routine, ICMAX1, is susceptible to an overflow error, because
the absolute value of a complex number can overflow even if the
real and imaginary parts are finite, even though the function being
computed (the index of the largest entry in absolute value) is well-defined
for all finite (or even non-NaN) inputs. So this is an example that can
return a wrong answer silently because of overflow (not NaNs).
-------- Original Message --------
Subject: Re: [Lapack] Possible bug in function icmax1
Date: Mon, 10 Feb 2014 13:19:24 -0700
From: Langou, Julien
To: slepc-maint maintenance
<lapack@Domain.Removed><mailto:lapack@Domain.Removed>, James Demmel
<demmel@Domain.Removed><mailto:demmel@Domain.Removed>, Jack Dongarra
<dongarra@Domain.Removed><mailto:dongarra@Domain.Removed>, Nicholas Higham
<nick.higham@Domain.Removed><mailto:nick.higham@Domain.Removed>, Sven J.
Dear Eloy, thanks for the bug report. There is a bug, but actually the bug is
not in the code but it is in the comments of the code. The routine ICMAX1 finds
the index of the first vector element of maximum absolute value. (And this is
not what the comment of the code is saying as you correctly points out.) There
is a historical reason for this comment mistake. (Related to us by Nick
Higham.) So the ICMAX1 routine initially was computing the index of the vector
element whose real part has maximum absolute value, (as the comments are
currently saying,) then we changed the routine so that it computes the index of
the vector element of maximum absolute value, and of course, forgot to change
the comments. Yup. A classic. Anyway. The comments are now corrected in our SVN
repository and should be out with the next LAPACK release. I am curious: do you
need a routine to the index of the vector element whose real part has maximum
absolute value (??)? I assume not and this was ju
uriosity on your part. Oh and just to be clear, there is the ICAMAX function in
the BLAS, this one works on | Re( . ) | + | Im( . ) |. The ICMAX1 function in
LAPACK works on | . | (so the ?genuine? absolute value of a complex number). We
tried to make this distinction clearer in the new commit. Good continuation
with SLEPc, this is a great package and much appreciated by the community.
On Feb 4, 2014, at 4:13 AM, Eloy Romero Alcalde
Dear LAPACK team,
The documentation of the function icmax1 says that it "finds the index
of the vector element whose real part has maximum absolute value", but I
think that it returns the index of the element with maximum absolute
value taking into account the imaginary part also. The next code tests
icmax1 with the complex array 1+99i, 2+98i, ..., 10+90i. The function
returns 1 using a LAPACK compiled with gfortran 4.6.3.
INTEGER INCX, N, R, I
INCX = 1
N = 10
DO 100 I=1,10
CX(I) = CMPLX(REAL(I), REAL(100-I))
R = ICMAX1(N, CX, INCX)
WRITE( *, * ) ' Result: ',R
My guess is that the problem is in the definition of the statement
function CABS1 in line 112 in icmax1.f:
* NEXT LINE IS THE ONLY MODIFICATION.
CABS1( ZDUM ) = ABS( ZDUM )
If it is replaced by the next one, the above code returns 10.
* NEXT LINE IS THE ONLY MODIFICATION.
CABS1( ZDUM ) = ABS( REAL(ZDUM) )
Possibly, there is a similar bug in the sister function izmax1.
Lapack mailing list
This concerns the comparison of cabs(x + i*y) vs. cabs(a + i*b) without
contamination by over/underflow, and without the computation (slowly) of
square roots. This problem arose about 62 years ago on IBM's 7090/7094
at the Univ. of Toronto. I do not have at hand the program we wrote to cope
with this problem, so I shall have to rely on an old man's unreliable memory.
Because I had rewritten IBM's Floating-Point Trap Handler, we could enable
and disable Gradual Underflow in less than a multiply time, and could also
compute an extended product p = (q1 + r1)*(q2 + r2)*(...)*(qN + rN) quickly
with NO over/underflow by returning both K and p/2.0**(128*K) . This product
capability is probably not available to you yet, so I shall do without it in
Preprocessing can ensure that x >= y >= 0 and a >= b >= 0 . Then
cabs(x + i*y) vs. cabs(a + i*b)
compares the same way (absent roundoff) as does
(x - a) vs. (b - y)*((0.5*b +
0.5*y)/(0.5*a + 0.5*x)) .
The latter expressions cannot suffer from overflow, nor from gradual underflow.
by more than roundoff, which interferes noticeably only when cabs(x + i*y)
and cabs(a + i*b) are too nearly equal, but this rarely matters.
I hope you find the foregoing recollection useful. For more such stuff see
-------------- next part --------------
An HTML attachment was scrubbed...