Hanging in ?HSEIN if NaN in H on input

Post here if you want to report a bug to the LAPACK team

Hanging in ?HSEIN if NaN in H on input

Postby akobotov » Wed Oct 03, 2012 4:59 am

Hi,

?HSEIN functions could hang if there is a NaN on input in the matrix H and eigenvalues in W much larger than SMLNUM
The hanging happens in loop 70 and, since the NaN causes the value of HNORM to be also NaN, the value of EPS3 is set to SMLNUM at lines 402-405. If eigenvalues are large enough the operation at line 417 WK = WK + EPS3 doesn't alter the eigenvalue since WK .EQ. WK + EPS3 just because EPS3 is below significant digits for WK. And the loop infinitely iterates over the same value of WK.

As I see the fix is simple given we have the new norm function with NaN checkers from LAPACK 3.4.2. Just check the received HNORM value for NaN and report error for the input parameter H. It could be something like this:
Code: Select all
-      LOGICAL            LSAME
+      LOGICAL            LSAME, SISNAN
       REAL               CLANHS, SLAMCH
-      EXTERNAL           LSAME, CLANHS, SLAMCH
+      EXTERNAL           LSAME, CLANHS, SLAMCH, SISNAN

                HNORM = CLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
                IF( HNORM.GT.RZERO ) THEN
                   EPS3 = HNORM*ULP
+               ELSE IF( SISNAN( HNORM ) ) THEN
+                  INFO = -6
+                  CALL XERBLA( 'CHSEIN', -INFO )
+                  RETURN
                ELSE
                   EPS3 = SMLNUM
                END IF


Thanks,
Alexander
akobotov
 
Posts: 8
Joined: Wed Feb 03, 2010 7:38 am
Location: Intel Corp., Russia, Novosibirsk

Re: Hanging in ?HSEIN if NaN in H on input

Postby rodney » Thu Oct 11, 2012 2:23 pm

Alexander,

I have added a check for NaN in the HNORM in revision 1350 of the LAPACK svn repository. The routine now will exit with INFO=-6, but no call to XERBLA, if a NaN is detected in the norm of H. Hopefully this will fix the problem.

Rodney
rodney
 
Posts: 49
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College

Re: Hanging in ?HSEIN if NaN in H on input

Postby admin » Thu Oct 11, 2012 5:46 pm

... and this is reported in our Errata list as bug0101 (http://www.netlib.org/lapack/Errata/index2.html)
Julie
admin
Site Admin
 
Posts: 504
Joined: Wed Dec 08, 2004 7:07 pm

Re: Hanging in ?HSEIN if NaN in H on input

Postby malxyz2 » Thu Dec 27, 2012 9:13 am

I think you also need to take care about positive & negative infinity.
The code example below would loop forever not only for NaN
but also for positive infinity.
So recent lapack changes (NaN checking) is not sufficient to avoid infinite loops.
The best way is to have explicit integer limit for such loops.

java example. FORTRAN results are identical

class k {
public static void main(String [] args)
{
double x=Double.POSITIVE_INFINITY;
//double x=Double.NaN;
ll:for(;;)
{
System.err.println("kkkk x="+x);
if(x<=0) break ll;
x=x/2;
}
}
}


kkkk x=Infinity
kkkk x=Infinity
kkkk x=Infinity
kkkk x=Infinity
............
malxyz2
 
Posts: 9
Joined: Tue Dec 25, 2012 2:28 pm

Re: Hanging in ?HSEIN if NaN in H on input

Postby malxyz2 » Thu Dec 27, 2012 9:23 am

fortran example of infinite loop with infinity


program f
double precision x
integer i
x=1.0

do i=1,100
x=x*1e9
end do
write(*,*) "x=",x

do while (.true.)
write(*,*) "lx=",x
x=x/2
if(x.le.0) exit
end do


end
malxyz2
 
Posts: 9
Joined: Tue Dec 25, 2012 2:28 pm

Re: Hanging in ?HSEIN if NaN in H on input

Postby malxyz2 » Thu Dec 27, 2012 9:36 am

also, to minimize lapack changes,
just the note

the loop

do while (.true.)
....
if(abs(x).le.eps) goto 10
end do

will loop forever if x id NaN or infinity
in the same time the loop

do while (.true.)
....
if(.not.(abs(x).gt.eps)) goto 10
end do

will exit on goto 10 when x is NaN of infinity
most lapack code is of type 1 loops.
rewriting it as in the second example would prevent NaN or infinity problem.

But again, I think that maximal integer limit is a better option.
malxyz2
 
Posts: 9
Joined: Tue Dec 25, 2012 2:28 pm

Re: Hanging in ?HSEIN if NaN in H on input

Postby rodney » Thu Dec 27, 2012 10:16 am

But again, I think that maximal integer limit is a better option.


I agree, we should just put in a maximum iteration count to prevent infinite loops instead of trying to detect NaN's and Inf's.

--Rodney
rodney
 
Posts: 49
Joined: Thu Feb 10, 2011 8:20 pm
Location: Colorado College


Return to Bug report

Who is online

Users browsing this forum: No registered users and 1 guest