gcc/testsuite/g77.f-torture/compile/980310-4.f

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