Fortran program freezes

[Li__Jun]

Hi,

I had an issue running a simulation program. The program run in loops and, within each loop, uses OpenML in a parallel way. The program can run more than 100 loops without any problem, but it freezes completely at a loop 136 and certain thread, without any error message.

After digging into the program, it appears that it freezes at the Poisson process step:

CALL S_RNPOI(LAMBDA, IR)

for that thread. The lambda value was 3.74, so it was a reasonable value.

My question is:

1) Do someone know why this happens, and if there is any direct way of resolving this?

2) Is there any indirect way to get around this? For example, try to run this line of code, but it takes more than 5 seconds, it runs an alternative code.

Thank you in advance!

Jun

[JohnNichols]

You should provide a sample that fails as observed.  You did not provide enough information.  

[Li__Jun]

Thanks, JohnNichols.

I just googled and found this post: https://community.intel.com/t5/Intel-Fortran-Compiler/Fortran-frozen-after-iterations-in-a-loop/td-p/941535

It appears that this is a bug in the function RNPOI in IMSL that trigger nonexit loop forever. The error depends on random seed and Poisson intensity, so it could be difficult to replicate on your end.

There is no way I can change the IMSL function, so I wonder if there is a way to do something like 2) I mentioned above?

[mecej4 · Accepted solution]

The IMSL error that you encountered was reported about nine years ago in this forum,  also 13 years ago in comp.lang.fortran, and pertains to an old version of IMSL (some platform distributions of versions 4, 5, 6 and early releases of 7). Which version of IMSL are you using?

You can find an extended discussion of this bug here . Note that this is an IMSL bug that was fixed many years ago, and has little to do with the Intel Fortran Compiler. Intel no longer provides a bundled version of IMSL, and the bug was no longer present in the last IMSL version (7.0.1) that was available from Intel as an optional add on to the Fortran compiler package.

The MKL library provides the same functionality. For details, read the MKL documentation for virngpoisson . After some initial set up, you can use the following lines of code to fill an array with the desired random numbers as follows:

errcode=vslnewstream( stream, brng, iseed ) 
do j = 1, nr 
   errcode = virngpoisson( method, stream, 1, y(j), lambda(j) ) 
end do

[Li__Jun]

Thank you so much, mecej4！ I followed your suggestion, managed to install MKL, and used the RNG there. Now it works well.

I am not sure the version of IMSL I have. The folder name is w_fcompxe_imsl_p_15.0.1.032. It came with my purchase of "Intel® Visual Fortran Composer XE for Windows* with IMSL*" about 10 years ago. Since the thread was back in 2013, so I think the version of IMSL has not fixed that issue.

[mecej4]

You have a workaround now, but I think that if you actually have installed IMSL from w_fcompxe_imsl_p_15.0.1.032 you should have version of 7.0.1 of IMSL and with it the program should not hang. I wonder if you have an older version of IMSL that gets linked because of your environmental setup.

Here is a short test program that you can run to help identify the version of IMSL that gets used. Please report the program output and also tell us if the program freezes.

!https://comp.lang.fortran.narkive.com/71xJ5w7A/why-does-the-program-stop-when-random-numbers-from-a-poisson-distribution-are-generated
!Hangs with CVF 6.6C + IMSL4
!OK with IFort 14 + IMSL 7.0.1
PROGRAM TEST
IMPLICIT NONE
INTEGER NR
PARAMETER (NR=7)
INTEGER IR(NR), ISEED
REAL MEAN
INTEGER    ISELCT, NOUT
CHARACTER  STRING(4)*80, TEMP*80, VERML*80
EXTERNAL   VERML
!
STRING(1) = '(" IMSL MATH/LIBRARY Version Number:  ", A)'
STRING(3) = '(" Fortran Compiler Version Number:  ", A)'
!                                Print the versions and numbers.
DO ISELCT=1, 3, 2
   TEMP = VERML(ISELCT)
   WRITE (*,STRING(ISELCT)) trim(TEMP)
END DO
MEAN = Z'41224925'
ISEED = 1013973114
print *,'Calling RNSET(',iseed,')'
CALL RNSET (ISEED)
print *,'Calling RNPOI'
CALL RNPOI(NR,MEAN,IR)
END PROGRAM

[Li__Jun]

Thanks, mecej4.

I just tried it. Somehow, I got an error message in running STRING(3), so I modified the code slightly so that it does not print out the fortran compiler version number. Below is what I got:

*************

IMSL MATH/LIBRARY Version Number: IMSL Fortran Numerical Library, Version 7.0.1
Calling RNSET( 1013973114 )
Calling RNPOI
Press any key to continue . . .

*************

It appears that it is version 7.0.1, and the program did finish.

But my original program does freeze even when I use different parameter values, seeds, and it stops at the function RNPOI. I guess  the problem remains at version 7.0.1 at different circumstances.

[Li__Jun]

Thanks, mercej4!

Below is what I got. I had to mute the STRING(3) as I got an error message. But it appears that the version is 7.0.1, and the program did finish. I guess the issue remains unfixed at different occasions.

IMSL MATH/LIBRARY Version Number: IMSL Fortran Numerical Library, Version 7.0.1
Calling RNSET( 1013973114 )
Calling RNPOI
Press any key to continue . . .

[Li__Jun]

Thanks, mercej4!

Below is what I got. I had to mute the STRING(3) as I got an error message. But it appears that the version is 7.0.1, and the program did finish. I guess the issue remains unfixed at different occasions.

IMSL MATH/LIBRARY Version Number: IMSL Fortran Numerical Library, Version 7.0.1
Calling RNSET( 1013973114 )
Calling RNPOI
Press any key to continue . . .

[Li__Jun]

I tried to post my reply several times without success. Did I violate some forum policy?

[nordmay0]

perhaps you tried to write an incorrect answer?

[Li__Jun]

Thanks, mercej4!

I had to mute the STRING(3) as I got an error message. But based on the outputs, it appears that the version is 7.0.1 and the program did finish without being frozen. I guess the issue remains unfixed at different occasions.

[mecej4]

The bug in the IMSL routine has since been fixed, but the reason for the failure is worth noting. The IMSL function RNPOI calls a uniform-distribution RNG, and solves a nonlinear equation to find the corresponding random number in a Poisson distribution. The solution fails to converge for certain values when the calculations are performed in single precision. Here is a simple demonstration program, which does not depend on having IMSL or MKL available. 

In each iteration, variables u and fp are updated, and the iterations are stopped when u falls below fp. In the old version of the IMSL library, the variable u was single-precision, and that caused the failure. To see this, change the declaration of u to single precision and run the program again.

In the old forum thread (posted in 2013) that Li__Jun linked to, there were some codes that were originally available for testing related issues, but the revamping of the forum may have caused those to be removed. 

program xrnp
   implicit none
   real :: ymu = 8.0, fp
   integer :: ix
   real*8 :: u !single precision causes failure
   integer :: i , j, iter
   real, external :: rnunf
!
   i=1
   ix = 0
   fp = exp(-ymu)
   u = 0.9999999    !value very near 1.0
   iter = 0
   do
      if( u <= fp )exit
      ix = ix + 1
      u = u - fp
      fp = fp*ymu/ix
      iter = iter+1
      print '(2i4,2x,2ES12.5)',iter,ix,u,fp
      if(iter <= 30)cycle  ! place a limit of 30 on iterations
      stop 'Failed'        ! u never went below fp
   end do
   stop
end program

[JohnNichols]

A brief look at the results with EXCEL for Real* 4 shows the problem. 

 program xrnp
   implicit none
   real :: ymu = 8.0, fp
   real :: ix
   real :: u !single precision causes failure
   integer :: i , j, iter
   real, external :: rnunf
!
   i=1
   ix = 0.0D0
   fp = exp(-ymu)
   u = 0.9999999    !value very near 1.0
   iter = 0
   do
      if( abs(u- fp) < 0.0000001 )exit
      ix = ix + 1.0D0
      u = u - fp
      fp = fp*ymu/ix
      iter = iter+1
      write(*, '(i4,2x,4ES12.5)'),iter,ymu,ix,u,fp
      if(iter <= 100)cycle  ! place a limit of 30 on iterations
      stop 'Failed'        ! u never went below fp
   end do
   stop
end program

Does solve the problem and is the way my very old and probably dead lecturers taught us how to do it.  

Conte and de Boore page 82 is the sort of thing I was taught. 

ERROR = ABS(B-A)
DO WHILE ERROR .GT. XTOL

[Li__Jun]

Yes, I can replicate the failure (with a small change, u = 0.99999999, not u = 0.9999999), but I have no clue about the bug in IMSL. I noticed this freeze more a year ago, but never doubted the problem was due to the RNPOI function until recently.

[Li__Jun]

Thanks, mercej4! Somehow I tried to post my reply several times without success.

Below is what I got. I had to mute the STRING(3) as I got an error message. But it appears that the version is 7.0.1, and the program did finish. I guess the issue remains unfixed at different occasions.

IMSL MATH/LIBRARY Version Number: IMSL Fortran Numerical Library, Version 7.0.1
Calling RNSET( 1013973114 )
Calling RNPOI
Press any key to continue .

[JohnNichols]

In my humble opinion, a human cannot write a RNG that is truly  random, it may pass the spectral tests, but that is also human.  

Even measuring some physical property that is supposedly random, is dependent on the generator.  

So randomness is a neat mathematical idea, it is purely an idea.  

[Steve_Lionel]

Just understand that Intel can't help you with IMSL anymore. Visual Numerics withdrew all rights to IMSL from Intel (Intel chose not to renew their contract.)

I'd definitely recommend an MKL solution if that works for you.

[nordmay0]

most likely this is due to an error in the RNPOI function