changeset 15579:5fb80374c881

quadgk.m: cleanup unwind_protect cruft
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Fri, 02 Nov 2012 15:18:49 -0400
parents 79083c78eac9
children 673889d49256
files scripts/general/quadgk.m
diffstat 1 files changed, 78 insertions(+), 89 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/general/quadgk.m	Sun Oct 28 07:36:13 2012 +0100
+++ b/scripts/general/quadgk.m	Fri Nov 02 15:18:49 2012 -0400
@@ -301,110 +301,99 @@
     endwhile
     subs = [subs(1:end-1), subs(2:end)];
 
-    # Not needed anmoyre
-    #warn_state = warning ("query", "Octave:divide-by-zero");
     # Set divide-by-zero warning off locally
     warning ("off", "Octave:divide-by-zero", "local");
 
     warn_msg   = "Octave:quadgk:warning-termination";
-    unwind_protect
-      ## Singularity will cause divide by zero warnings
-      warning ("off", "Octave:divide-by-zero");
+
+    ## Singularity will cause divide by zero warnings
+    warning ("off", "Octave:divide-by-zero");
+
+    ## 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);
-      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
 
-      if (isa (a, "single") || isa (b, "single") || isa (waypoints, "single"))
-        myeps = eps ("single");
-      else
-        myeps = eps;
+      ## Quit if any evaluations are not finite (Inf or NaN)
+      if (any (! isfinite (q_subs)))
+        warning (warn_msg, "quadgk: non finite integrand encountered");
+        q = q0;
+        err = err0;
+        break;
+      endif
+
+      tol = max (abstol, reltol .* abs (q0));
+
+      ## If the global error estimate is meet exit
+      if (err0 < tol)
+        q = q0;
+        err = err0;
+        break;
       endif
 
-      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 (warn_msg, "quadgk: non finite integrand encountered");
-          q = q0;
-          err = err0;
-          break;
-        endif
-
-        tol = max (abstol, reltol .* abs (q0));
+      ## Accept the subintervals that meet the convergence criteria
+      idx = find (abs (q_errs) < tol .* abs (diff (subs, [], 2)) ./ h);
+      if (first)
+        q = sum (q_subs (idx));
+        err = sum (q_errs(idx));
+        first = false;
+      else
+        q0 = q + sum (q_subs);
+        err0 = err + sum (q_errs);
+        q += sum (q_subs (idx));
+        err += sum (q_errs(idx));
+      endif
+      subs(idx,:) = [];
 
-        ## If the global error estimate is meet exit
-        if (err0 < tol)
-          q = q0;
-          err = err0;
-          break;
-        endif
+      ## If no remaining subintervals exit
+      if (rows (subs) == 0)
+        break;
+      endif
 
-        ## Accept the subintervals that meet the convergence criteria
-        idx = find (abs (q_errs) < tol .* abs (diff (subs, [], 2)) ./ h);
-        if (first)
-          q = sum (q_subs (idx));
-          err = sum (q_errs(idx));
-          first = false;
-        else
-          q0 = q + sum (q_subs);
-          err0 = err + sum (q_errs);
-          q += sum (q_subs (idx));
-          err += sum (q_errs(idx));
-        endif
-        subs(idx,:) = [];
+      if (trace)
+        disp ([rows(subs), err, q0]);
+      endif
 
-        ## If no remaining subintervals exit
-        if (rows (subs) == 0)
-          break;
-        endif
-
-        if (trace)
-          disp ([rows(subs), err, q0]);
-        endif
-
-        ## Split remaining subintervals in two
-        mid = (subs(:,2) + subs(:,1)) ./ 2;
-        subs = [subs(:,1), mid; mid, subs(:,2)];
+      ## Split remaining subintervals in two
+      mid = (subs(:,2) + subs(:,1)) ./ 2;
+      subs = [subs(:,1), mid; mid, subs(:,2)];
 
-        ## If the maximum subinterval count is met accept remaining
-        ## subinterval and exit
-        if (rows (subs) > maxint)
-          warning (warn_msg,
-                             "quadgk: maximum interval count (%d) met", maxint);
-          q += sum (q_subs);
-          err += sum (q_errs);
-          break;
-        endif
-
-        ## Evaluation of the integrand on the remaining subintervals
-        [q_subs, q_errs] = __quadgk_eval__ (f, subs);
-      endwhile
+      ## If the maximum subinterval count is met accept remaining
+      ## subinterval and exit
+      if (rows (subs) > maxint)
+        warning (warn_msg,
+                 "quadgk: maximum interval count (%d) met", maxint);
+        q += sum (q_subs);
+        err += sum (q_errs);
+        break;
+      endif
 
-      if (err > max (abstol, reltol * abs (q)))
-        warning (warn_msg,
-                    "quadgk: Error tolerance not met. Estimated error %g", err);
-      endif
-    unwind_protect_cleanup
+      ## Evaluation of the integrand on the remaining subintervals
+      [q_subs, q_errs] = __quadgk_eval__ (f, subs);
+    endwhile
 
-      # not needed aynomre, used local off of warnings
-      %{
-      if (strcmp (warn_state.state, "on"))
-        warning ("on", "Octave:divide-by-zero");
-      endif
-      %}
+    if (err > max (abstol, reltol * abs (q)))
+      warning (warn_msg,
+               "quadgk: Error tolerance not met. Estimated error %g", err);
+    endif
 
-    end_unwind_protect
   endif
 endfunction