Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran 95: WRITE to screen changes output to WRITE to file

  1. Jul 30, 2010 #1
    I have a short code that calculates the running median for three columns of data in an input file and writes the running medians out to three columns of data in an output file. I use Numerical Recipe subroutines and functions to compute the running median. If I include a WRITE(*,*) directly after the computation of each running median and then write all three running medians to the output file using WRITE(UNIT=outputfile,FMT=*), then I get the expected result. If I comment out the WRITE(*,*) commands, then the output repeats the running median of the first column in the other two columns. I'm running Mac OS10.6, using gfortran to compile and I can reproduce the effect on three different machines. Has anyone seen this before? Is there a better workaround than writing to screen?
    PROGRAM RunningMedianFinder

    IMPLICIT NONE

    INTEGER,PARAMETER :: inputfile=10, outputfile=20

    INTEGER :: rint,i,k,ios

    REAL :: f(1025),h1(1025),h2(1025),h3(1025), &
    h1m,h2m,h3m,fm,selip,medians(3)

    CHARACTER(LEN=100) :: inputname,outputname

    ! First we find the name of the input file (raw data) and
    ! the output file (running median), and then we open them.

    WRITE(*,*) 'What is the name of the input file (raw data)?'
    READ(*,*) inputname

    WRITE(*,*) 'What is the name of the output file (running medians)?'
    READ(*,*) outputname

    OPEN(UNIT=inputfile,FILE=inputname,POSITION='REWIND')
    OPEN(UNIT=outputfile,FILE=outputname,POSITION='REWIND')

    ! Here, we determine the length of the running median -- called rint
    ! -- for Running INTerval.

    WRITE(*,*) 'What is the length of the interval for the running median?'
    READ(*,*) rint

    ! First, initialize the arrays.

    DO i=1,1025,1
    f(i)=0
    h1(i)=0
    h2(i)=0
    h3(i)=0
    END DO

    ! Read in the first interval for the running meadian.

    DO i=1,rint,1
    READ(UNIT=inputfile,FMT=*) f(i),h1(i),h2(i),h3(i)
    END DO

    ! Because we are finding a median, the value of k used in selip will
    ! be half the running interval.

    k = rint/2

    ! We will assign the value of the median to be at the frequency in
    ! the middle of the running interval.


    ! Now we start the loop over the entire file and record the running
    ! medians. We use the ios variable to determine when we have reached
    ! the end of the inputfile

    ios = 0

    DO WHILE (ios.ge.0)

    h1m = selip(k,rint,h1)
    WRITE(*,*) h1m
    h2m = selip(k,rint,h2)
    WRITE(*,*) h2m
    h3m = selip(k,rint,h3)
    WRITE(*,*) h3m

    fm = f(k)

    WRITE(UNIT=outputfile,FMT=*) fm,h1m,h2m,h3m

    DO i=2,rint,1
    f(i-1)=f(i)
    h1(i-1)=h1(i)
    h2(i-1)=h2(i)
    h3(i-1)=h3(i)
    END DO

    READ(UNIT=inputfile,FMT=*,IOSTAT=ios) f(rint),h1(rint),h2(rint),h3(rint)
    END DO

    CLOSE(UNIT=inputfile)
    CLOSE(UNIT=outputfile)

    END PROGRAM

    !!!!!!SUBROUTINES!!!!!!!!!

    SUBROUTINE shell(n,s)

    IMPLICIT NONE

    INTEGER :: n
    REAL :: s(n)

    INTEGER :: i,j,inc
    REAL :: v

    inc=1

    DO WHILE(inc.le.n)
    inc=3*inc+1
    END DO
    DO WHILE(inc.gt.1)
    inc=inc/3
    DO i=inc+1,n
    v=s(i)
    j=i
    DO WHILE(s(j-inc).gt.v.and.j.gt.inc)
    s(j)=s(j-inc)
    j=j-inc
    END DO
    s(j)=v
    END DO
    END DO

    RETURN

    END SUBROUTINE


    !!!!!!!FUNCTIONS!!!!!!!!!!

    FUNCTION selip(k,n,arr)

    IMPLICIT NONE

    INTEGER,PARAMETER :: M=64
    INTEGER :: k,n
    REAL,PARAMETER :: BIG=1.e30
    REAL :: selip,arr(n)

    INTEGER :: i,j,jl,jm,ju,kk,mm,nlo,nxtmm,isel(M+2),a
    REAL :: ahi,alo,sum,sel(M+2)

    kk=k
    ahi=BIG
    alo=-BIG
    DO WHILE(a.ne.1)
    a=1
    mm=0
    nlo=0
    sum=0.
    nxtmm=M+1
    DO i=1,n
    IF(arr(i).ge.alo.and.arr(i).le.ahi) THEN
    mm=mm+1
    IF(arr(i).eq.alo) nlo=nlo+1
    IF(mm.le.M) THEN
    sel(mm)=arr(i)
    ELSE IF(mm.eq.nxtmm) THEN
    nxtmm=mm+mm/M
    sel(1+mod(i+mm+kk,M))=arr(i)
    END IF
    sum=sum+arr(i)
    END IF
    END DO
    IF(kk.le.nlo) THEN
    selip=alo
    RETURN
    ELSE IF(mm.le.M) THEN
    CALL shell(mm,sel)
    selip=sel(kk)
    RETURN
    END IF
    sel(M+1)=sum/mm
    CALL shell(M+1,sel)
    sel(M+2)=ahi
    DO j=1,M+2
    isel(j)=0
    END DO
    DO i=1,n
    IF(arr(i).ge.alo.and.arr(i).le.ahi) THEN
    jl=0
    ju=M+2
    DO WHILE(ju-jl.gt.1)
    jm=(ju+jl)/2
    IF(arr(i).ge.sel(jm)) THEN
    jl=jm
    ELSE
    ju=jm
    END IF
    END DO
    isel(ju)=isel(ju)+1
    END IF
    END DO
    j=1
    DO WHILE(kk.gt.isel(j))
    alo=sel(j)
    kk=kk-isel(j)
    j=j+1
    END DO
    ahi=sel(j)
    a=0
    END DO

    END FUNCTION
     
  2. jcsd
  3. Jul 30, 2010 #2

    Mark44

    Staff: Mentor

    So are you saying that if you comment out the three WRITE (*,*) statements in this DO loop
    Code (Text):

    DO WHILE (ios.ge.0)

    h1m = selip(k,rint,h1)
    WRITE(*,*) h1m
    h2m = selip(k,rint,h2)
    WRITE(*,*) h2m
    h3m = selip(k,rint,h3)
    WRITE(*,*) h3m

    fm = f(k)

    WRITE(UNIT=outputfile,FMT=*) fm,h1m,h2m,h3m

    DO i=2,rint,1
    f(i-1)=f(i)
    h1(i-1)=h1(i)
    h2(i-1)=h2(i)
    h3(i-1)=h3(i)
    END DO

    READ(UNIT=inputfile,FMT=*,IOSTAT=ios) f(rint),h1(rint),h2(rint),h3(rint)
    END DO
     
    then the WRITE(UNIT=outputfile, FMT=* ) statement produces the screwy results you described? If that's the case, you might consider using an actual format statement for that last WRITE statement.
     
  4. Jul 30, 2010 #3
    Thanks for the suggestion. What you describe is precisely the problem. Unfortunately, I tried an actual format statement (using FMT='(E7.6)' ) and I still get the screwy results unless I comment out the three WRITE(*,*) statements. Perhaps this is a bug in the gfortan compiler?
     
  5. Jul 30, 2010 #4

    alphysicist

    User Avatar
    Homework Helper

    Hi benacquista,

    I believe your problem is in the selip function:


    I believe your next line is the problem; you are using the variable a without it being initialized. If right here you will insert:

    Code (Text):
    a=0.
    I think that will fix your problem.

     
  6. Jul 31, 2010 #5

    Mark44

    Staff: Mentor

    This is the opposite of what you said earlier.

    In the first quote above, you're now saying that you don't get the screwy results if you comment out the three WRITE statements.
     
  7. Jul 31, 2010 #6

    Mark44

    Staff: Mentor

    alphysicist, I agree that the selip function should initialize a before testing it in the DO WHILE loop, but I don't see how selip returning a bogus value would affect the formatting of a WRITE statement.
     
  8. Jul 31, 2010 #7
    Sorry! This has got me so confused, I'm confusing others. My first statement is correct. If the WRITE statements are there, then it works. If I comment out the WRITE statements, then I get weird results.
     
  9. Jul 31, 2010 #8
    I tried alphysicist's suggestion and it worked. I'm trying to figure out why, and my only guess is that the WRITE statement causes the stored value of 'a' to change ... otherwise, it remains equal to 1 (it's value when selip returns a number). If it remains equal to 1, then the following calls to selip do nothing and simply return the previous value of selip.

    Thank you very much!
     
  10. Jul 31, 2010 #9

    alphysicist

    User Avatar
    Homework Helper

    I don't know anything about compilers, so I can't say I know "why" it works. However, years of painful experience has taught me that when something weird is happening, it quite often is an unitialized variable using whatever garbage happens to be at the memory location.

    That's the right idea, but remember there is not really a "stored value" for a between calls to selip. After the first call to selip is over, the memory location for a had a value of 1. On the second call to selip (when the function calls are immediately after each other), a is apparently assigned the identical memory location, which has not been disturbed. Since you have not initialized a, it uses whatever garbage is already there, which happens to be your previous value.

    And when you output something between calls, my guess is that a is assigned a different memory location, and now the garbage there is not equal to 1, and your function works again. (But as I said above, I don't know anything about the executables created by the fortran compiler; if anyone knows more details about this I would like to hear it.)
     
  11. Aug 3, 2010 #10

    DrDu

    User Avatar
    Science Advisor

    This problem sounds all to familiar to me. Is there any compiler option to warn about un-initialized variables?
     
  12. Aug 4, 2010 #11

    alphysicist

    User Avatar
    Homework Helper

    Hi DrDu,

    I don't use gfortran, but with g95 you can use the -Wall switch to catch these simple cases of uninitialized variables. gfortran might use the same option.
     
  13. Aug 4, 2010 #12
    I tried this by commenting out the place where I initialized the variable, and compiled with gfortran using -Wall. It didn't catch the uninitalized variable, but it flagged every one of my tabs as a "Nonconforming tab character". I'm not quite sure what that means, but since I tab liberally in order to make it easy for me to follow the code structure, it flagged every line in the code. Probably the best solution is simply to be diligent about initializing variables, and targeting that as the most likely problem when debugging.
     
  14. Aug 4, 2010 #13

    alphysicist

    User Avatar
    Homework Helper

    You can stop the warning about tabs by adding -Wtabs as an extra compiler option (in addition to the -Wall).

    But was -Wall catching uninitialized variables for gfortran? I've recently installed Ubuntu on one of my machines and was planning on using gfortran.
     
  15. Aug 4, 2010 #14
    Thanks, after using -Wtabs (and commenting out the initialization of the variable) I could see the one other warning. Unfortunately, gfortran did not catch the uninitialized variable. It did catch an unused variable in the declarations of the main program.
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook