Fortran: problem with NEWUNIT on different machines

el_duderino · 03-13-2015, 04:45 PM

Hello all,

I'm writing a program in Fortran 90 which uses random number distributions. To generate the distributions, I'm using the following program (just as a test..):

Code:

program test_rand
 implicit none

 integer :: seed, i, j, k, kk

 real :: r

 CALL init_random_seed()         ! see example of RANDOM_SEED
 
 do j = 1, 5
   kk = 0
    do k = 1, 100
       do i = 1, 1000
          kk = kk + 1
          CALL RANDOM_NUMBER(r)
          write(1000+j,*) kk, r
       enddo
    enddo
 enddo
 
 stop
end program test_rand

subroutine init_random_seed()
  implicit none
  integer, allocatable :: seed(:)
  integer :: i, n, un, istat, dt(8), pid, t(2), s
  integer(8) :: count, tms

  call random_seed(size = n)
  allocate(seed(n))
  ! First try if the OS provides a random number generator
  open(newunit=un, file="/dev/urandom", access="stream", &
       form="unformatted", action="read", status="old", iostat=istat)
  if (istat == 0) then
     read(un) seed
     close(un)
  else
     ! Fallback to XOR:ing the current time and pid. The PID is
     ! useful in case one launches multiple instances of the same
     ! program in parallel.
     call system_clock(count)
     if (count /= 0) then
        t = transfer(count, t)
     else
        call date_and_time(values=dt)
        tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 &
             + dt(2) * 31_8 * 24 * 60 * 60 * 1000 &
             + dt(3) * 24 * 60 * 60 * 60 * 1000 &
             + dt(5) * 60 * 60 * 1000 &
             + dt(6) * 60 * 1000 + dt(7) * 1000 &
             + dt(8)
        t = transfer(tms, t)
     end if
     s = ieor(t(1), t(2))
     pid = getpid() + 1099279 ! Add a prime
     s = ieor(s, pid)
     if (n >= 3) then
        seed(1) = t(1) + 36269
        seed(2) = t(2) + 72551
        seed(3) = pid
        if (n > 3) then
           seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /)
        end if
     else
        seed = s + 37 * (/ (i, i = 0, n - 1 ) /)
     end if
  end if
  call random_seed(put=seed)
end subroutine init_random_seed

Now, if I try to compile this small program on my laptop with gfortran and Ubuntu 14.10, no problem. Same as for another laptop, Ubuntu 10.10, and for a Red Hat Cluster. However, if I try to complile it on a Mac, I have the following error message:

Code:

open(newunit=un,file="/dev/urandom", access="stream", &
             1
Error: Syntax error in OPEN statement at (1)

I'm really struggling trying to understand why the Mac machine does not like the 'newunit' statement. In fact, it is also not even fonted in vim.

Based on the code posted, if I substitute 'newunit' with 'unit', I have different results if I ask to print the variables 'unit' and 'seed' in my laptop and in the Mac machine, and I do not understand why.

Could anyone please help me?

Thanks

suicidaleggroll · 03-13-2015, 05:12 PM

newunit is not a valid argument to open in fortran 90:
ftp://ftp.nag.co.uk/sc22wg5/N001-N1100/N692.pdf

It was introduced in fortran 2008:
http://www.j3-fortran.org/doc/year/10/10-007.pdf

A good reference for the Fortran standards is here (that's where those previous two links are from):
https://gcc.gnu.org/wiki/GFortranStandards

As for the different values, of course, unit is an input to open, newunit is an output. When you change newunit to unit, it means you're passing in an undefined variable for the unit number, which may have unintended side effects when you go to use it.

Just change newunit=un to unit=10 (or [almost] any other number, or set un=10 on a previous line and use unit=un).

el_duderino · 03-14-2015, 09:41 AM

Thank! Now I have similar results in all the machines.
Probably, in the Mac machine the fortran compiler is too old, while newer compiler on other machines do recognize newunit. Am I right?

However, I still don't understand why the results are so different specifying un = 10 from the case in which I do not specify it in the Mac machine. Do you have any idea?

suicidaleggroll · 03-14-2015, 10:01 AM

Yes gfortran on the Linux machine is probably newer and has been updated with the 2008 standards.

As for the different results - you were passing an uninitialized variable to the unit number input. Unit numbers can not be anything you want, if the wrong number was in that memory location (eg: stdin, stdout, stderr) it could cause open to behave erratically.

el_duderino · 03-14-2015, 11:20 AM

Very well, thanks for your help and the explainations!