multiroots error checking.

[Toby White <tow at theor dot ch dot cam dot ac dot uk>]

The documentation for multiroots indicates that gsl_multiroot_fd?f_iterate
functions should return error codes if the function evaluates badly:

     These functions perform a single iteration of the solver S.  If the
     iteration encounters an unexpected problem then an error code will
     be returned,

    `GSL_EBADFUNC'
          the iteration encountered a singular point where the function
          or its derivative evaluated to `Inf' or `NaN'.

As far as I can see from the source, this does not happen. This is rather
annoying since I can't stop iterating if the current guess has strayed outside
of my function's domain of validity.

The following patch adds this behaviour; if the function (or derivative) 
provided returns an error code other than GSL_SUCCESS, then 
gsl_multiroot_fd?f_iterate() will return GSL_EBADFUNC.


diff -ur multiroots.orig/broyden.c multiroots/broyden.c
--- multiroots.orig/broyden.c	Thu Nov 30 17:51:40 2000
+++ multiroots/broyden.c	Thu Nov 30 17:50:09 2000
@@ -230,6 +230,7 @@
 {
   broyden_state_t *state = (broyden_state_t *) vstate;
 
+  int status;
   double phi0, phi1, t, lambda;
 
   gsl_matrix *H = state->H;
@@ -273,7 +274,8 @@
       gsl_vector_set (x_trial, i, xi + t * pi);
     }
 
-  GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
+  status=GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
+  if (status) return GSL_EBADFUNC;
 
   phi1 = enorm (fnew);
 
@@ -313,7 +315,8 @@
           gsl_vector_set (x_trial, i, xi + t * pi);
         }
       
-      GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
+      status=GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);
+      if (status) GSL_ERROR("error evaluating function", GSL_EBADFUNC);
       
       phi1 = enorm (fnew);
     }
diff -ur multiroots.orig/dnewton.c multiroots/dnewton.c
--- multiroots.orig/dnewton.c	Thu Nov 30 17:51:40 2000
+++ multiroots/dnewton.c	Thu Nov 30 17:50:14 2000
@@ -117,6 +117,8 @@
   
   int signum ;
 
+  int status;
+
   size_t i;
 
   size_t n = function->n ;
@@ -135,8 +137,9 @@
       gsl_vector_set (x, i, y - e);
     }
   
-  GSL_MULTIROOT_FN_EVAL (function, x, f);
-  
+  status=GSL_MULTIROOT_FN_EVAL (function, x, f);
+  if (status) return GSL_EBADFUNC;
+ 
   gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->J);
 
   return GSL_SUCCESS;
diff -ur multiroots.orig/fdjac.c multiroots/fdjac.c
--- multiroots.orig/fdjac.c	Thu Nov 30 17:51:40 2000
+++ multiroots/fdjac.c	Thu Nov 30 17:56:42 2000
@@ -29,6 +29,7 @@
   const size_t m = f->size;
   const size_t n1 = jacobian->size1;
   const size_t n2 = jacobian->size2;
+  int status;
 
   if (m != n1 || n != n2)
     {
@@ -68,7 +69,8 @@
 	  }
 
 	gsl_vector_set (x1, j, xj + dx);
-	GSL_MULTIROOT_FN_EVAL (F, x1, f1);
+	status=GSL_MULTIROOT_FN_EVAL (F, x1, f1);
+        if (status) return GSL_EBADFUNC;
 	gsl_vector_set (x1, j, xj);
 
 	for (i = 0; i < m; i++)
diff -ur multiroots.orig/hybrid.c multiroots/hybrid.c
--- multiroots.orig/hybrid.c	Thu Nov 30 17:51:40 2000
+++ multiroots/hybrid.c	Thu Nov 30 17:50:27 2000
@@ -409,6 +409,7 @@
 {
   hybrid_state_t *state = (hybrid_state_t *) vstate;
 
+  int status;
   const double fnorm = state->fnorm;
 
   gsl_matrix *J = state->J;
@@ -454,8 +455,9 @@
 
   /* Evaluate function at x + p */
 
-  GSL_MULTIROOT_FN_EVAL (func, x_trial, f_trial);
-
+  status=GSL_MULTIROOT_FN_EVAL (func, x_trial, f_trial);
+  if (status) return GSL_EBADFUNC;
+  
   /* Set df = f_trial - f */
 
   compute_df (f_trial, f, df);
diff -ur multiroots.orig/hybridj.c multiroots/hybridj.c
--- multiroots.orig/hybridj.c	Thu Nov 30 17:51:40 2000
+++ multiroots/hybridj.c	Thu Nov 30 17:50:29 2000
@@ -381,6 +381,7 @@
 {
   hybridj_state_t *state = (hybridj_state_t *) vstate;
 
+  int status;
   const double fnorm = state->fnorm;
 
   gsl_matrix *q = state->q;
@@ -425,7 +426,8 @@
 
   /* Evaluate function at x + p */
 
-  GSL_MULTIROOT_FN_EVAL_F (fdf, x_trial, f_trial);
+  status=GSL_MULTIROOT_FN_EVAL_F (fdf, x_trial, f_trial);
+  if (status) return GSL_EBADFUNC;
 
   /* Set df = f_trial - f */
 
@@ -498,7 +500,8 @@
 
   if (state->ncfail == 2)
     {
-      GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J);
+      status=GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J);
+      if (status) return GSL_EBADFUNC;
 
       state->nslow2++;
 
diff -ur multiroots.orig/newton.c multiroots/newton.c
--- multiroots.orig/newton.c	Thu Nov 30 17:51:40 2000
+++ multiroots/newton.c	Thu Nov 30 17:50:39 2000
@@ -96,7 +96,7 @@
 {
   newton_state_t * state = (newton_state_t *) vstate;
   
-  int signum ;
+  int signum, status ;
 
   size_t i;
 
@@ -116,7 +116,8 @@
       gsl_vector_set (x, i, y - e);
     }
 
-  GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J);
+  status=GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J);
+  if (status) return GSL_EBADFUNC;
 
   return GSL_SUCCESS;
 }


-- 
Toby White, University Chemical Lab., Lensfield Road, Cambridge. CB2 1EW. U.K.
Email: <tow@theor.ch.cam.ac.uk> GPG Key ID: 1DE9DE75
Web: <URL:http://ket.ch.cam.ac.uk/people/tow/index.html>
Tel: +44 1223 336423
Fax: +44 1223 336362