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

AI Thread Summary
The discussion centers around a Fortran program designed to calculate running medians for three columns of data from an input file and output these medians to an output file. The issue arises when the WRITE statements used for debugging are commented out, leading to incorrect output where the median values for the second and third columns repeat the first column's median. This behavior is reproducible across multiple machines using gfortran on Mac OS 10.6. Participants suggest that the problem may stem from an uninitialized variable in the selip function, specifically the variable 'a', which is used in a loop condition. Initializing 'a' resolves the issue, indicating that the WRITE statements inadvertently affect the program's flow by altering the memory state. Compiler options like -Wall are discussed for catching uninitialized variables, but it appears that gfortran does not always flag such issues. The conversation emphasizes the importance of diligent variable initialization and debugging practices in programming.
benacquista
Messages
6
Reaction score
0
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
 
Technology news on Phys.org
So are you saying that if you comment out the three WRITE (*,*) statements in this DO loop
Code:
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.
 
Mark44 said:
So are you saying that if you comment out the three WRITE (*,*) statements in this DO loop
Code:
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.

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?
 
Hi benacquista,

I believe your problem is in the selip function:


benacquista said:
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

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:
a=0.

I think that will fix your problem.

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
 
benacquista said:
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?
This is the opposite of what you said earlier.

benacquista said:
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.
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.
 
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.
 
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.
 
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!
 
Mark44 said:
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.

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.

benacquista said:
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!

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.)
 
  • #10
This problem sounds all to familiar to me. Is there any compiler option to warn about un-initialized variables?
 
  • #11
Hi DrDu,

DrDu said:
This problem sounds all to familiar to me. Is there any compiler option to warn about un-initialized variables?

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.
 
  • #12
alphysicist said:
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.

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.
 
  • #13
benacquista said:
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.

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.
 
  • #14
alphysicist said:
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.

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.
 

Similar threads

Replies
12
Views
3K
Replies
6
Views
3K
Replies
9
Views
1K
Replies
5
Views
5K
Replies
2
Views
2K
Replies
3
Views
3K
Replies
11
Views
3K
Back
Top