- #1
benacquista
- 6
- 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
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