LAPACK Archives

[Lapack] Possible bug in function icmax1

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 
<wkahan@Domain.Removed<mailto:wkahan@Domain.Removed>> wrote:

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).
Jim



-------- Original Message --------
Subject:        Re: [Lapack] Possible bug in function icmax1
Date:   Mon, 10 Feb 2014 13:19:24 -0700
From:   Langou, Julien 
<Julien.Langou@Domain.Removed><mailto:Julien.Langou@Domain.Removed>
To:     slepc-maint maintenance 
<slepc-maint@Domain.Removed><mailto:slepc-maint@Domain.Removed>
CC:     lapack@Domain.Removed<mailto:lapack@Domain.Removed> 
<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. 
Hammarling <sven@Domain.Removed><mailto:sven@Domain.Removed>



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
st c
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. 
Cheers, Julien.



On Feb 4, 2014, at 4:13 AM, Eloy Romero Alcalde 
<elroal@Domain.Removed><mailto:elroal@Domain.Removed> wrote:

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 ICMAX1
      EXTERNAL ICMAX1

      INTEGER INCX, N, R, I
      COMPLEX CX(10)
      INCX = 1
      N = 10
      DO 100 I=1,10
        CX(I) = CMPLX(REAL(I), REAL(100-I))
 100  CONTINUE

      R = ICMAX1(N, CX, INCX)
      WRITE( *, * ) ' Result: ',R
      END

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.

Regards,
Eloy Romero
_______________________________________________
Lapack mailing list
Lapack@Domain.Removed<mailto:Lapack@Domain.Removed>
http://lists.eecs.utk.edu/mailman/listinfo/lapack


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 
what
follows.

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
www.eecs.berkeley.edu/~wkahan/7094II.pdf<http://www.eecs.berkeley.edu/~wkahan/7094II.pdf>
                              W. Kahan


-------------- next part --------------
An HTML attachment was scrubbed...
URL: 
<http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20140214/2eb2677c/attachment.html>

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or