I'm not going to get into registering to submit a bug in the formal way,
I am trying urgently to get some quick and dirty results.
But the gfortran compiler told me to report a bug, so I am doing so
Command line
[EMAIL PROTECTED] QDPfig]# gfortran -o contour_fun
contour_fun.f90
Error message
contour_fun.f90: In function 'findval':
contour_fun.f90:91: internal compiler error: in gfc_conv_variable, at
fortran/trans-expr.c:403
Please submit a full bug report,
with preprocessed source if appropriate.
See <URL:http://gcc.gnu.org/bugs.html> for instructions.
Version
[EMAIL PROTECTED] QDPfig]# gfortran -v
Using built-in specs.
Target: i386-apple-darwin8
Configured with: ../gcc-4.2-20060610/configure --prefix=/sw
--prefix=/sw/lib/gcc4 --enable-languages=c,c++,fortran,objc,java
--with-gmp=/sw --with-included-gettext --host=i386-apple-darwin8
--with-as=/sw/lib/odcctools/bin/as --with-ld=/sw/lib/odcctools/bin/ld
--with-nm=/sw/lib/odcctools/bin/nm --mandir=/sw/lib/gcc4/share/man
--infodir=/sw/lib/gcc4/share/info --disable-multilib
Thread model: posix
gcc version 4.2.0 20060610 (experimental)
Code
Attached
Regards
Gerry Skinner
PROGRAM countour_fun
!
! write a qdp file to draw contour line of the likelihood and chisqd functions
! for the few photon coded mask case
!
IMPLICIT NONE
!
INTEGER, PARAMETER :: ncntr = 4
REAL :: cntr_levels(ncntr)=(/16.0, 25.0, 36.0, 49.0/)
INTEGER, PARAMETER :: nf=2
REAL :: f_values(2) = (/0.4,0.6/)
REAL :: c1, c0, f
REAL :: c1_min=1.0, c1_max=100.0
REAL :: c0_min=1.0, c0_max=100.0, c0step = 1.0
REAL :: fmin, fmax
INTEGER :: n, fun, k
!
REAL, EXTERNAL :: findval, funget
!
!---
!
OPEN ( 1, FILE='contour_fun.qdp', STATUS='unknown')
WRITE (1,'(A)') ' skip on'
f=f_values(1)
c1=c1_min
DO WHILE (c1 < c1_max)
write (2,*) c1 , funget ( 1, c1, c0_min, f )
c1 = c1 + 1.0
end DO
DO n = 1,2
DO fun = 1,nf
f = f_values(fun)
DO k = 1, ncntr
c0 = c0_min
DO WHILE (c0 < c0_max )
! seek c1 which gives the required level of function f with current c0
c1 = findval( fun, c0, f, c1_min, c1_max, cntr_levels(k))
IF (c1 > 0.0) THEN
WRITE (1,'(2F10.4)') c1, c0
END IF
c0 = c0 + c0step
END DO
WRITE (1,'(A)') ' no no no'
END DO
END DO ! next function
f = fmax
END DO ! next f value
END PROGRAM countour_fun
REAL FUNCTION findval( fun, c0, f, c1_min, c1_max, contour)
!
! finds point on c1 axis at which function number fun
! has value contour when other parameter are c0, f
IMPLICIT NONE
!
INTEGER, INTENT(in) :: fun
REAL, INTENT(in) :: c0
REAL, INTENT(in) :: f
REAL, INTENT(in) :: c1_min, c1_max
REAL, INTENT(in) :: contour
!
INTEGER, PARAMETER :: maxloop=20
REAL, PARAMETER :: tol = 1.0E-3
!
INTEGER :: n
REAL :: c1a, c1b, c
REAL :: v_1, v_2, v
!
REAL, EXTERNAL :: funget
!
c1a = c1_min
c1b = c1_max
DO n = 1, maxloop
v_1 = funget ( fun, c1a, c0, f )
v_2 = funget ( fun, c1b, c0, f )
IF ( ABS (v_2 - v_1) < tol ) THEN
findval = (c1a+c1b)/2.0
RETURN
END IF
IF (( v_1 > contour).OR. (v_2< contour) ) THEN
! write (*,'(i6,9F12.4)') n, contour, c1a, c, c1b, v_1, v, v_2
! WRITE (*,*) ' outside range', fun, c0, f
STOP
funget = -99.9
RETURN
END IF
IF ( v_1 > v_2 ) THEN
write (*,'(i6,9F12.4)') n, contour, c1a, c, c1b, v_1, v, v_2
WRITE (*,*) ' reverse slope', fun, c0, f
STOP
END IF
c = ( c1a + c1b ) /2.0
v = funget ( fun, c, c0, f )
write (*,'(i6,9F12.4)') n, contour, c1a, c, c1b, v_1, v, v_2
IF ( v >= contour ) THEN ! in first half
c1b = c
ELSE
c1a = c
END IF
END DO
END FUNCTION findval
REAL FUNCTION funget( fun, c1, c0, f )
!
! returns function number fun
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: fun
REAL, INTENT(in) :: c1, c0
REAL, INTENT(in) :: f
!
REAL, EXTERNAL :: cash, chi2
!
!---
!
IF (fun == 1 ) THEN
funget = cash( c1, c0, f )
ELSE IF (fun == 2 ) THEN
funget = chi2( c1, c0, f )
ELSE
WRITE (*,*) ' Bad fun'
STOP
END IF
END FUNCTION funget
REAL FUNCTION cash( c1, c0, f )
!
! returns cash function
!
IMPLICIT NONE
!
REAL, INTENT(in) :: c1, c0
REAL, INTENT(in) :: f
!
REAL :: b, s
!
!---
!
b = c0 /(1.0-f)
s = (c1/f) - b
IF (s < 0.0 ) THEN
cash = 0.0
ELSE
cash = (2*(c1*log(c1/f)+c0*log(c0/(1.0-f))-(c1+c0)*log(c1+c0)))
END IF
END FUNCTION cash
REAL FUNCTION chi2( c1, c0, f )
!
! returns chisqd function
!
IMPLICIT NONE
!
REAL, INTENT(in) :: c1, c0
REAL, INTENT(in) :: f
!
REAL :: e0, e1, b, s
!---
!
b = c1 /(1.0-f)
s = (c0/f) - b
IF (s < 0.0 ) THEN
chi2 = 0.0
ELSE
e1 = f * (c0+c1)
e0 = (1-f) * (c0+c1)
chi2 = (((e1-c1)**2.0)/e1) + (((e0-c0)**2.0)/e0)
END IF
END FUNCTION chi2