1c { dg-do compile } 2C To: egcs-bugs@cygnus.com 3C Subject: -fPIC problem showing up with fortran on x86 4C From: Dave Love <d.love@dl.ac.uk> 5C Date: 19 Dec 1997 19:31:41 +0000 6C 7C 8C This illustrates a long-standing problem noted at the end of the g77 9C `Actual Bugs' info node and thought to be in the back end. Although 10C the report is against gcc 2.7 I can reproduce it (specifically on 11C redhat 4.2) with the 971216 egcs snapshot. 12C 13C g77 version 0.5.21 14C gcc -v -fnull-version -o /tmp/gfa00415 -xf77-cpp-input /tmp/gfa00415.f -xnone 15C -lf2c -lm 16C 17 18C ------------ 19 subroutine dqage(f,a,b,epsabs,epsrel,limit,result,abserr, 20 * neval,ier,alist,blist,rlist,elist,iord,last) 21C -------------------------------------------------- 22C 23C Modified Feb 1989 by Barry W. Brown to eliminate key 24C as argument (use key=1) and to eliminate all Fortran 25C output. 26C 27C Purpose: to make this routine usable from within S. 28C 29C -------------------------------------------------- 30c***begin prologue dqage 31c***date written 800101 (yymmdd) 32c***revision date 830518 (yymmdd) 33c***category no. h2a1a1 34c***keywords automatic integrator, general-purpose, 35c integrand examinator, globally adaptive, 36c gauss-kronrod 37c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 38c de doncker,elise,appl. math. & progr. div. - k.u.leuven 39c***purpose the routine calculates an approximation result to a given 40c definite integral i = integral of f over (a,b), 41c hopefully satisfying following claim for accuracy 42c abs(i-reslt).le.max(epsabs,epsrel*abs(i)). 43c***description 44c 45c computation of a definite integral 46c standard fortran subroutine 47c double precision version 48c 49c parameters 50c on entry 51c f - double precision 52c function subprogram defining the integrand 53c function f(x). the actual name for f needs to be 54c declared e x t e r n a l in the driver program. 55c 56c a - double precision 57c lower limit of integration 58c 59c b - double precision 60c upper limit of integration 61c 62c epsabs - double precision 63c absolute accuracy requested 64c epsrel - double precision 65c relative accuracy requested 66c if epsabs.le.0 67c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 68c the routine will end with ier = 6. 69c 70c key - integer 71c key for choice of local integration rule 72c a gauss-kronrod pair is used with 73c 7 - 15 points if key.lt.2, 74c 10 - 21 points if key = 2, 75c 15 - 31 points if key = 3, 76c 20 - 41 points if key = 4, 77c 25 - 51 points if key = 5, 78c 30 - 61 points if key.gt.5. 79c 80c limit - integer 81c gives an upperbound on the number of subintervals 82c in the partition of (a,b), limit.ge.1. 83c 84c on return 85c result - double precision 86c approximation to the integral 87c 88c abserr - double precision 89c estimate of the modulus of the absolute error, 90c which should equal or exceed abs(i-result) 91c 92c neval - integer 93c number of integrand evaluations 94c 95c ier - integer 96c ier = 0 normal and reliable termination of the 97c routine. it is assumed that the requested 98c accuracy has been achieved. 99c ier.gt.0 abnormal termination of the routine 100c the estimates for result and error are 101c less reliable. it is assumed that the 102c requested accuracy has not been achieved. 103c error messages 104c ier = 1 maximum number of subdivisions allowed 105c has been achieved. one can allow more 106c subdivisions by increasing the value 107c of limit. 108c however, if this yields no improvement it 109c is rather advised to analyze the integrand 110c in order to determine the integration 111c difficulties. if the position of a local 112c difficulty can be determined(e.g. 113c singularity, discontinuity within the 114c interval) one will probably gain from 115c splitting up the interval at this point 116c and calling the integrator on the 117c subranges. if possible, an appropriate 118c special-purpose integrator should be used 119c which is designed for handling the type of 120c difficulty involved. 121c = 2 the occurrence of roundoff error is 122c detected, which prevents the requested 123c tolerance from being achieved. 124c = 3 extremely bad integrand behavior occurs 125c at some points of the integration 126c interval. 127c = 6 the input is invalid, because 128c (epsabs.le.0 and 129c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 130c result, abserr, neval, last, rlist(1) , 131c elist(1) and iord(1) are set to zero. 132c alist(1) and blist(1) are set to a and b 133c respectively. 134c 135c alist - double precision 136c vector of dimension at least limit, the first 137c last elements of which are the left 138c end points of the subintervals in the partition 139c of the given integration range (a,b) 140c 141c blist - double precision 142c vector of dimension at least limit, the first 143c last elements of which are the right 144c end points of the subintervals in the partition 145c of the given integration range (a,b) 146c 147c rlist - double precision 148c vector of dimension at least limit, the first 149c last elements of which are the 150c integral approximations on the subintervals 151c 152c elist - double precision 153c vector of dimension at least limit, the first 154c last elements of which are the moduli of the 155c absolute error estimates on the subintervals 156c 157c iord - integer 158c vector of dimension at least limit, the first k 159c elements of which are pointers to the 160c error estimates over the subintervals, 161c such that elist(iord(1)), ..., 162c elist(iord(k)) form a decreasing sequence, 163c with k = last if last.le.(limit/2+2), and 164c k = limit+1-last otherwise 165c 166c last - integer 167c number of subintervals actually produced in the 168c subdivision process 169c 170c***references (none) 171c***routines called d1mach,dqk15,dqk21,dqk31, 172c dqk41,dqk51,dqk61,dqpsrt 173c***end prologue dqage 174c 175 double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b, 176 * blist,b1,b2,dabs,defabs,defab1,defab2,dmax1,d1mach,elist,epmach, 177 * epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f, 178 * resabs,result,rlist,uflow 179 integer ier,iord,iroff1,iroff2,k,last,limit,maxerr,neval, 180 * nrmax 181c 182 dimension alist(limit),blist(limit),elist(limit),iord(limit), 183 * rlist(limit) 184c 185 external f 186c 187c list of major variables 188c ----------------------- 189c 190c alist - list of left end points of all subintervals 191c considered up to now 192c blist - list of right end points of all subintervals 193c considered up to now 194c rlist(i) - approximation to the integral over 195c (alist(i),blist(i)) 196c elist(i) - error estimate applying to rlist(i) 197c maxerr - pointer to the interval with largest 198c error estimate 199c errmax - elist(maxerr) 200c area - sum of the integrals over the subintervals 201c errsum - sum of the errors over the subintervals 202c errbnd - requested accuracy max(epsabs,epsrel* 203c abs(result)) 204c *****1 - variable for the left subinterval 205c *****2 - variable for the right subinterval 206c last - index for subdivision 207c 208c 209c machine dependent constants 210c --------------------------- 211c 212c epmach is the largest relative spacing. 213c uflow is the smallest positive magnitude. 214c 215c***first executable statement dqage 216 epmach = d1mach(4) 217 uflow = d1mach(1) 218c 219c test on validity of parameters 220c ------------------------------ 221c 222 ier = 0 223 neval = 0 224 last = 0 225 result = 0.0d+00 226 abserr = 0.0d+00 227 alist(1) = a 228 blist(1) = b 229 rlist(1) = 0.0d+00 230 elist(1) = 0.0d+00 231 iord(1) = 0 232 if(epsabs.le.0.0d+00.and. 233 * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) ier = 6 234 if(ier.eq.6) go to 999 235c 236c first approximation to the integral 237c ----------------------------------- 238c 239 neval = 0 240 call dqk15(f,a,b,result,abserr,defabs,resabs) 241 last = 1 242 rlist(1) = result 243 elist(1) = abserr 244 iord(1) = 1 245c 246c test on accuracy. 247c 248 errbnd = dmax1(epsabs,epsrel*dabs(result)) 249 if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2 250 if(limit.eq.1) ier = 1 251 if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs) 252 * .or.abserr.eq.0.0d+00) go to 60 253c 254c initialization 255c -------------- 256c 257c 258 errmax = abserr 259 maxerr = 1 260 area = result 261 errsum = abserr 262 nrmax = 1 263 iroff1 = 0 264 iroff2 = 0 265c 266c main do-loop 267c ------------ 268c 269 do 30 last = 2,limit 270c 271c bisect the subinterval with the largest error estimate. 272c 273 a1 = alist(maxerr) 274 b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) 275 a2 = b1 276 b2 = blist(maxerr) 277 call dqk15(f,a1,b1,area1,error1,resabs,defab1) 278 call dqk15(f,a2,b2,area2,error2,resabs,defab2) 279c 280c improve previous approximations to integral 281c and error and test for accuracy. 282c 283 neval = neval+1 284 area12 = area1+area2 285 erro12 = error1+error2 286 errsum = errsum+erro12-errmax 287 area = area+area12-rlist(maxerr) 288 if(defab1.eq.error1.or.defab2.eq.error2) go to 5 289 if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12) 290 * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1 291 if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1 292 5 rlist(maxerr) = area1 293 rlist(last) = area2 294 errbnd = dmax1(epsabs,epsrel*dabs(area)) 295 if(errsum.le.errbnd) go to 8 296c 297c test for roundoff error and eventually set error flag. 298c 299 if(iroff1.ge.6.or.iroff2.ge.20) ier = 2 300c 301c set error flag in the case that the number of subintervals 302c equals limit. 303c 304 if(last.eq.limit) ier = 1 305c 306c set error flag in the case of bad integrand behavior 307c at a point of the integration range. 308c 309 if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03* 310 * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3 311c 312c append the newly-created intervals to the list. 313c 314 8 if(error2.gt.error1) go to 10 315 alist(last) = a2 316 blist(maxerr) = b1 317 blist(last) = b2 318 elist(maxerr) = error1 319 elist(last) = error2 320 go to 20 321 10 alist(maxerr) = a2 322 alist(last) = a1 323 blist(last) = b1 324 rlist(maxerr) = area2 325 rlist(last) = area1 326 elist(maxerr) = error2 327 elist(last) = error1 328c 329c call subroutine dqpsrt to maintain the descending ordering 330c in the list of error estimates and select the subinterval 331c with the largest error estimate (to be bisected next). 332c 333 20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) 334c ***jump out of do-loop 335 if(ier.ne.0.or.errsum.le.errbnd) go to 40 336 30 continue 337c 338c compute final result. 339c --------------------- 340c 341 40 result = 0.0d+00 342 do 50 k=1,last 343 result = result+rlist(k) 344 50 continue 345 abserr = errsum 346 60 neval = 30*neval+15 347 999 return 348 end 349