1 SUBROUTINE dqagpe(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,
2 * ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,PTS,IORD,LEVEL,NDIN,
192 DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
193 * A2,B,BLIST,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,DMIN1,
194 * dres,
d1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,
195 * errmax,error1,erro12,error2,errsum,ertest,oflow,points,pts,
196 * resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp,uflow
197 INTEGER I,ID,IER,IERRO,IND1,IND2,IORD,IP1,IROFF1,IROFF2,IROFF3,J,
198 * JLOW,JUPBND,K,KSGN,KTMIN,LAST,LEVCUR,LEVEL,LEVMAX,LIMIT,MAXERR,
199 * ndin,neval,nint,nintp1,npts,npts2,nres,nrmax,numrl2
203 dimension alist(limit),blist(limit),elist(limit),iord(limit),
204 * level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3),
205 * rlist(limit),rlist2(52)
278 IF(npts2.LT.2.OR.limit.LE.npts.OR.(epsabs.LE.0.0d+00.AND.
279 * epsrel.LT.dmax1(0.5d+02*epmach,0.5d-28))) ier = 6
280 IF(ier.EQ.6)
GO TO 999
286 IF(a.GT.b) sign = -1.0d+00
288 IF(npts.EQ.0)
GO TO 15
292 15 pts(npts+2) = dmax1(a,b)
295 IF(npts.EQ.0)
GO TO 40
300 IF(pts(i).LE.pts(j))
GO TO 20
305 IF(pts(1).NE.dmin1(a,b).OR.pts(nintp1).NE.dmax1(a,b)) ier = 6
306 IF(ier.EQ.6)
GO TO 999
314 CALL dqk21(f,a1,b1,area1,error1,defabs,resa,ier)
315 IF (ier .LT. 0)
RETURN 316 abserr = abserr+error1
317 result = result+area1
319 IF(error1.EQ.resa.AND.error1.NE.0.0d+00) ndin(i) = 1
320 resabs = resabs+defabs
331 IF(ndin(i).EQ.1) elist(i) = abserr
332 errsum = errsum+elist(i)
340 errbnd = dmax1(epsabs,epsrel*dres)
341 IF(abserr.LE.0.1d+03*epmach*resabs.AND.abserr.GT.errbnd) ier = 2
342 IF(nint.EQ.1)
GO TO 80
348 IF(elist(ind1).GT.elist(ind2))
GO TO 60
352 IF(ind1.EQ.iord(i))
GO TO 70
356 IF(limit.LT.npts2) ier = 1
357 80
IF(ier.NE.0.OR.abserr.LE.errbnd)
GO TO 210
364 errmax = elist(maxerr)
383 IF(dres.GE.(0.1d+01-0.5d+02*epmach)*resabs) ksgn = 1
388 DO 160 last = npts2,limit
393 levcur = level(maxerr)+1
395 b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
399 CALL dqk21(f,a1,b1,area1,error1,resa,defab1,ier)
400 IF (ier .LT. 0)
RETURN 401 CALL dqk21(f,a2,b2,area2,error2,resa,defab2,ier)
402 IF (ier .LT. 0)
RETURN 409 erro12 = error1+error2
410 errsum = errsum+erro12-errmax
411 area = area+area12-rlist(maxerr)
412 IF(defab1.EQ.error1.OR.defab2.EQ.error2)
GO TO 95
413 IF(dabs(rlist(maxerr)-area12).GT.0.1d-04*dabs(area12)
414 * .OR.erro12.LT.0.99d+00*errmax)
GO TO 90
415 IF(extrap) iroff2 = iroff2+1
416 IF(.NOT.extrap) iroff1 = iroff1+1
417 90
IF(last.GT.10.AND.erro12.GT.errmax) iroff3 = iroff3+1
418 95 level(maxerr) = levcur
420 rlist(maxerr) = area1
422 errbnd = dmax1(epsabs,epsrel*dabs(area))
426 IF(iroff1+iroff2.GE.10.OR.iroff3.GE.20) ier = 2
427 IF(iroff2.GE.5) ierro = 3
432 IF(last.EQ.limit) ier = 1
437 IF(dmax1(dabs(a1),dabs(b2)).LE.(0.1d+01+0.1d+03*epmach)*
438 * (dabs(a2)+0.1d+04*uflow)) ier = 4
442 IF(error2.GT.error1)
GO TO 100
446 elist(maxerr) = error1
449 100 alist(maxerr) = a2
452 rlist(maxerr) = area2
454 elist(maxerr) = error2
461 110
CALL dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
463 IF(errsum.LE.errbnd)
GO TO 190
465 IF(ier.NE.0)
GO TO 170
467 erlarg = erlarg-erlast
468 IF(levcur+1.LE.levmax) erlarg = erlarg+erro12
474 IF(level(maxerr)+1.LE.levmax)
GO TO 160
477 120
IF(ierro.EQ.3.OR.erlarg.LE.ertest)
GO TO 140
485 IF(last.GT.(2+limit/2)) jupbnd = limit+3-last
488 errmax = elist(maxerr)
490 IF(level(maxerr)+1.LE.levmax)
GO TO 160
496 140 numrl2 = numrl2+1
497 rlist2(numrl2) = area
498 IF(numrl2.LE.2)
GO TO 155
499 CALL dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
501 IF(ktmin.GT.5.AND.abserr.LT.0.1d-02*errsum) ier = 5
502 IF(abseps.GE.abserr)
GO TO 150
507 ertest = dmax1(epsabs,epsrel*dabs(reseps))
509 IF(abserr.LT.ertest)
GO TO 170
513 150
IF(numrl2.EQ.1) noext = .true.
514 IF(ier.GE.5)
GO TO 170
516 errmax = elist(maxerr)
527 170
IF(abserr.EQ.oflow)
GO TO 190
528 IF((ier+ierro).EQ.0)
GO TO 180
529 IF(ierro.EQ.3) abserr = abserr+correc
531 IF(result.NE.0.0d+00.AND.area.NE.0.0d+00)
GO TO 175
532 IF(abserr.GT.errsum)
GO TO 190
533 IF(area.EQ.0.0d+00)
GO TO 210
535 175
IF(abserr/dabs(result).GT.errsum/dabs(area))
GO TO 190
539 180
IF(ksgn.EQ.(-1).AND.dmax1(dabs(result),dabs(area)).LE.
540 * resabs*0.1d-01)
GO TO 210
541 IF(0.1d-01.GT.(result/area).OR.(result/area).GT.0.1d+03.OR.
542 * errsum.GT.dabs(area)) ier = 6
549 result = result+rlist(k)
552 210
IF(ier.GT.2) ier = ier-1
double precision function d1mach(i)
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)