| c { dg-do compile } |
| C To: egcs-bugs@cygnus.com |
| C Subject: -fPIC problem showing up with fortran on x86 |
| C From: Dave Love <d.love@dl.ac.uk> |
| C Date: 19 Dec 1997 19:31:41 +0000 |
| C |
| C |
| C This illustrates a long-standing problem noted at the end of the g77 |
| C `Actual Bugs' info node and thought to be in the back end. Although |
| C the report is against gcc 2.7 I can reproduce it (specifically on |
| C redhat 4.2) with the 971216 egcs snapshot. |
| C |
| C g77 version 0.5.21 |
| C gcc -v -fnull-version -o /tmp/gfa00415 -xf77-cpp-input /tmp/gfa00415.f -xnone |
| C -lf2c -lm |
| C |
| |
| C ------------ |
| subroutine dqage(f,a,b,epsabs,epsrel,limit,result,abserr, |
| * neval,ier,alist,blist,rlist,elist,iord,last) |
| C -------------------------------------------------- |
| C |
| C Modified Feb 1989 by Barry W. Brown to eliminate key |
| C as argument (use key=1) and to eliminate all Fortran |
| C output. |
| C |
| C Purpose: to make this routine usable from within S. |
| C |
| C -------------------------------------------------- |
| c***begin prologue dqage |
| c***date written 800101 (yymmdd) |
| c***revision date 830518 (yymmdd) |
| c***category no. h2a1a1 |
| c***keywords automatic integrator, general-purpose, |
| c integrand examinator, globally adaptive, |
| c gauss-kronrod |
| c***author piessens,robert,appl. math. & progr. div. - k.u.leuven |
| c de doncker,elise,appl. math. & progr. div. - k.u.leuven |
| c***purpose the routine calculates an approximation result to a given |
| c definite integral i = integral of f over (a,b), |
| c hopefully satisfying following claim for accuracy |
| c abs(i-reslt).le.max(epsabs,epsrel*abs(i)). |
| c***description |
| c |
| c computation of a definite integral |
| c standard fortran subroutine |
| c double precision version |
| c |
| c parameters |
| c on entry |
| c f - double precision |
| c function subprogram defining the integrand |
| c function f(x). the actual name for f needs to be |
| c declared e x t e r n a l in the driver program. |
| c |
| c a - double precision |
| c lower limit of integration |
| c |
| c b - double precision |
| c upper limit of integration |
| c |
| c epsabs - double precision |
| c absolute accuracy requested |
| c epsrel - double precision |
| c relative accuracy requested |
| c if epsabs.le.0 |
| c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), |
| c the routine will end with ier = 6. |
| c |
| c key - integer |
| c key for choice of local integration rule |
| c a gauss-kronrod pair is used with |
| c 7 - 15 points if key.lt.2, |
| c 10 - 21 points if key = 2, |
| c 15 - 31 points if key = 3, |
| c 20 - 41 points if key = 4, |
| c 25 - 51 points if key = 5, |
| c 30 - 61 points if key.gt.5. |
| c |
| c limit - integer |
| c gives an upperbound on the number of subintervals |
| c in the partition of (a,b), limit.ge.1. |
| c |
| c on return |
| c result - double precision |
| c approximation to the integral |
| c |
| c abserr - double precision |
| c estimate of the modulus of the absolute error, |
| c which should equal or exceed abs(i-result) |
| c |
| c neval - integer |
| c number of integrand evaluations |
| c |
| c ier - integer |
| c ier = 0 normal and reliable termination of the |
| c routine. it is assumed that the requested |
| c accuracy has been achieved. |
| c ier.gt.0 abnormal termination of the routine |
| c the estimates for result and error are |
| c less reliable. it is assumed that the |
| c requested accuracy has not been achieved. |
| c error messages |
| c ier = 1 maximum number of subdivisions allowed |
| c has been achieved. one can allow more |
| c subdivisions by increasing the value |
| c of limit. |
| c however, if this yields no improvement it |
| c is rather advised to analyze the integrand |
| c in order to determine the integration |
| c difficulties. if the position of a local |
| c difficulty can be determined(e.g. |
| c singularity, discontinuity within the |
| c interval) one will probably gain from |
| c splitting up the interval at this point |
| c and calling the integrator on the |
| c subranges. if possible, an appropriate |
| c special-purpose integrator should be used |
| c which is designed for handling the type of |
| c difficulty involved. |
| c = 2 the occurrence of roundoff error is |
| c detected, which prevents the requested |
| c tolerance from being achieved. |
| c = 3 extremely bad integrand behavior occurs |
| c at some points of the integration |
| c interval. |
| c = 6 the input is invalid, because |
| c (epsabs.le.0 and |
| c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), |
| c result, abserr, neval, last, rlist(1) , |
| c elist(1) and iord(1) are set to zero. |
| c alist(1) and blist(1) are set to a and b |
| c respectively. |
| c |
| c alist - double precision |
| c vector of dimension at least limit, the first |
| c last elements of which are the left |
| c end points of the subintervals in the partition |
| c of the given integration range (a,b) |
| c |
| c blist - double precision |
| c vector of dimension at least limit, the first |
| c last elements of which are the right |
| c end points of the subintervals in the partition |
| c of the given integration range (a,b) |
| c |
| c rlist - double precision |
| c vector of dimension at least limit, the first |
| c last elements of which are the |
| c integral approximations on the subintervals |
| c |
| c elist - double precision |
| c vector of dimension at least limit, the first |
| c last elements of which are the moduli of the |
| c absolute error estimates on the subintervals |
| c |
| c iord - integer |
| c vector of dimension at least limit, the first k |
| c elements of which are pointers to the |
| c error estimates over the subintervals, |
| c such that elist(iord(1)), ..., |
| c elist(iord(k)) form a decreasing sequence, |
| c with k = last if last.le.(limit/2+2), and |
| c k = limit+1-last otherwise |
| c |
| c last - integer |
| c number of subintervals actually produced in the |
| c subdivision process |
| c |
| c***references (none) |
| c***routines called d1mach,dqk15,dqk21,dqk31, |
| c dqk41,dqk51,dqk61,dqpsrt |
| c***end prologue dqage |
| c |
| double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b, |
| * blist,b1,b2,dabs,defabs,defab1,defab2,dmax1,d1mach,elist,epmach, |
| * epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f, |
| * resabs,result,rlist,uflow |
| integer ier,iord,iroff1,iroff2,k,last,limit,maxerr,neval, |
| * nrmax |
| c |
| dimension alist(limit),blist(limit),elist(limit),iord(limit), |
| * rlist(limit) |
| c |
| external f |
| c |
| c list of major variables |
| c ----------------------- |
| c |
| c alist - list of left end points of all subintervals |
| c considered up to now |
| c blist - list of right end points of all subintervals |
| c considered up to now |
| c rlist(i) - approximation to the integral over |
| c (alist(i),blist(i)) |
| c elist(i) - error estimate applying to rlist(i) |
| c maxerr - pointer to the interval with largest |
| c error estimate |
| c errmax - elist(maxerr) |
| c area - sum of the integrals over the subintervals |
| c errsum - sum of the errors over the subintervals |
| c errbnd - requested accuracy max(epsabs,epsrel* |
| c abs(result)) |
| c *****1 - variable for the left subinterval |
| c *****2 - variable for the right subinterval |
| c last - index for subdivision |
| c |
| c |
| c machine dependent constants |
| c --------------------------- |
| c |
| c epmach is the largest relative spacing. |
| c uflow is the smallest positive magnitude. |
| c |
| c***first executable statement dqage |
| epmach = d1mach(4) |
| uflow = d1mach(1) |
| c |
| c test on validity of parameters |
| c ------------------------------ |
| c |
| ier = 0 |
| neval = 0 |
| last = 0 |
| result = 0.0d+00 |
| abserr = 0.0d+00 |
| alist(1) = a |
| blist(1) = b |
| rlist(1) = 0.0d+00 |
| elist(1) = 0.0d+00 |
| iord(1) = 0 |
| if(epsabs.le.0.0d+00.and. |
| * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) ier = 6 |
| if(ier.eq.6) go to 999 |
| c |
| c first approximation to the integral |
| c ----------------------------------- |
| c |
| neval = 0 |
| call dqk15(f,a,b,result,abserr,defabs,resabs) |
| last = 1 |
| rlist(1) = result |
| elist(1) = abserr |
| iord(1) = 1 |
| c |
| c test on accuracy. |
| c |
| errbnd = dmax1(epsabs,epsrel*dabs(result)) |
| if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2 |
| if(limit.eq.1) ier = 1 |
| if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs) |
| * .or.abserr.eq.0.0d+00) go to 60 |
| c |
| c initialization |
| c -------------- |
| c |
| c |
| errmax = abserr |
| maxerr = 1 |
| area = result |
| errsum = abserr |
| nrmax = 1 |
| iroff1 = 0 |
| iroff2 = 0 |
| c |
| c main do-loop |
| c ------------ |
| c |
| do 30 last = 2,limit |
| c |
| c bisect the subinterval with the largest error estimate. |
| c |
| a1 = alist(maxerr) |
| b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) |
| a2 = b1 |
| b2 = blist(maxerr) |
| call dqk15(f,a1,b1,area1,error1,resabs,defab1) |
| call dqk15(f,a2,b2,area2,error2,resabs,defab2) |
| c |
| c improve previous approximations to integral |
| c and error and test for accuracy. |
| c |
| neval = neval+1 |
| area12 = area1+area2 |
| erro12 = error1+error2 |
| errsum = errsum+erro12-errmax |
| area = area+area12-rlist(maxerr) |
| if(defab1.eq.error1.or.defab2.eq.error2) go to 5 |
| if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12) |
| * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1 |
| if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1 |
| 5 rlist(maxerr) = area1 |
| rlist(last) = area2 |
| errbnd = dmax1(epsabs,epsrel*dabs(area)) |
| if(errsum.le.errbnd) go to 8 |
| c |
| c test for roundoff error and eventually set error flag. |
| c |
| if(iroff1.ge.6.or.iroff2.ge.20) ier = 2 |
| c |
| c set error flag in the case that the number of subintervals |
| c equals limit. |
| c |
| if(last.eq.limit) ier = 1 |
| c |
| c set error flag in the case of bad integrand behavior |
| c at a point of the integration range. |
| c |
| if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03* |
| * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3 |
| c |
| c append the newly-created intervals to the list. |
| c |
| 8 if(error2.gt.error1) go to 10 |
| alist(last) = a2 |
| blist(maxerr) = b1 |
| blist(last) = b2 |
| elist(maxerr) = error1 |
| elist(last) = error2 |
| go to 20 |
| 10 alist(maxerr) = a2 |
| alist(last) = a1 |
| blist(last) = b1 |
| rlist(maxerr) = area2 |
| rlist(last) = area1 |
| elist(maxerr) = error2 |
| elist(last) = error1 |
| c |
| c call subroutine dqpsrt to maintain the descending ordering |
| c in the list of error estimates and select the subinterval |
| c with the largest error estimate (to be bisected next). |
| c |
| 20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) |
| c ***jump out of do-loop |
| if(ier.ne.0.or.errsum.le.errbnd) go to 40 |
| 30 continue |
| c |
| c compute final result. |
| c --------------------- |
| c |
| 40 result = 0.0d+00 |
| do 50 k=1,last |
| result = result+rlist(k) |
| 50 continue |
| abserr = errsum |
| 60 neval = 30*neval+15 |
| 999 return |
| end |
