# HG changeset patch
# User David Bateman
# Date 1351199704 7200
# Node ID f9ac079d152f018fe21cf94d0c14dae6cba63d69
# Parent 48b1312bb0fe497218be21dcd53599d6527bf909
Fix early exit from quadgk for infinite intervals
* quadgk.m : Move test for small intervals from the main loop to
__quadgk_eval__ and perform the test on the abscissa rather than the
limits of the subintervals
diff git a/scripts/general/quadgk.m b/scripts/general/quadgk.m
 a/scripts/general/quadgk.m
+++ b/scripts/general/quadgk.m
@@ 303,28 +303,19 @@
## Singularity will cause divide by zero warnings
warning ("off", "Octave:dividebyzero");
 ## Initial evaluation of the integrand on the subintervals
 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
 q0 = sum (q_subs);
 err0 = sum (q_errs);

if (isa (a, "single")  isa (b, "single")  isa (waypoints, "single"))
myeps = eps ("single");
else
myeps = eps;
endif
+ ## Initial evaluation of the integrand on the subintervals
+ [q_subs, q_errs] = __quadgk_eval__ (f, subs, myeps);
+ q0 = sum (q_subs);
+ err0 = sum (q_errs);
+
first = true;
while (true)
 ## Check for subintervals that are too small. Test must be
 ## performed in untransformed subintervals. What is a good
 ## value for this test. Shampine suggests 100*eps
 if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
 q = q0;
 err = err0;
 break;
 endif

## Quit if any evaluations are not finite (Inf or NaN)
if (any (! isfinite (q_subs)))
warning ("quadgk: non finite integrand encountered");
@@ 379,7 +370,16 @@
endif
## Evaluation of the integrand on the remaining subintervals
 [q_subs, q_errs] = __quadgk_eval__ (f, subs);
+ [q_subs, q_errs, tooclose] = __quadgk_eval__ (f, subs, myeps);
+
+ if (tooclose)
+ ## Quit as __quadgk_eval__ flagged that the
+ ## subintervals were too close.
+ q = q0;
+ err = err0;
+ break;
+ endif
+
endwhile
if (err > max (abstol, reltol * abs (q)))
@@ 393,7 +393,7 @@
endif
endfunction
function [q, err] = __quadgk_eval__ (f, subs)
+function [q, err, tooclose] = __quadgk_eval__ (f, subs, myeps)
## A (15,7) point pair of GaussKonrod quadrature rules. The abscissa
## and weights are copied directly from dqk15w.f from quadpack
@@ 423,8 +423,20 @@
0.1294849661688697e+00]);
halfwidth = diff (subs, [], 2) ./ 2;
 center = sum (subs, 2) ./ 2;;
+ center = sum (subs, 2) ./ 2;
x = bsxfun (@plus, halfwidth * abscissa, center);
+
+ ## Check for subintervals that are too small. Test must be
+ ## performed in untransformed subintervals. What is a good
+ ## value for this test. Shampine suggests 100*eps
+ if (any (abs (bsxfun (@(y,z) y ./ z, diff (x, [], 2), (x(:,end)x(:,1))))(:)
+ < 100 * myeps))
+ tooclose = true;
+ return;
+ else
+ tooclose = false;
+ endif
+
y = reshape (f (x(:)), size (x));
## This is faster than using bsxfun as the * operator can use a