Hanging in ?HSEIN if NaN in H on input

[akobotov]

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

[rodney]

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

[admin]

... and this is reported in our Errata list as bug0101 (http://www.netlib.org/lapack/Errata/index2.html)
Julie

[malxyz2]

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]

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]

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.

[rodney]

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