Additional multiroots functions.

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

As it stands, to initialise a multiroots solver, there are two
functions available;

gsl_multiroot_fsolver_alloc()

which allocates *and sets* a solver, and

gsl_multiroot_fsolver_set()

which sets an already allocated solver.

------

However, in several situations it is useful to simply allocate
a solver, without yet providing a function and initial guess.

The following patch provides the following additional functions:

gsl_multiroot_fsolver * 
gsl_multiroot_fsolver_alloc_no_func (const gsl_multiroot_fsolver_type * T, 
                                     size_t n); 

which allocates space for an n-dimensional fsolver.

gsl_multiroot_fdfsolver *
gsl_multiroot_fdfsolver_alloc_no_func(const gsl_multiroot_fdfsolver_type * T,
                                      size_t n);

which allocates space for an n-dimensional fdfsolver.

It also adds some error checking to the _iterate functions to ensure
that the solver has been set before use; and some further error-checking
to the _set functions to ensure that the solver, function, and initial
guess are all of the same dimensionality.

Personally, I find it immensely useful. Any chance of it being it put
in gsl?

------

diff -ur multiroots.orig/fdfsolver.c multiroots/fdfsolver.c
--- multiroots.orig/fdfsolver.c	Mon Feb 12 16:20:56 2001
+++ multiroots/fdfsolver.c	Mon Feb 12 16:39:47 2001
@@ -29,17 +29,37 @@
                                gsl_vector * x)
 {
   int status;
+  
+  gsl_multiroot_fdfsolver * s = NULL;
 
-  gsl_multiroot_fdfsolver * s;
-
-  const size_t n = f->n ;
+  s = gsl_multiroot_fdfsolver_alloc_no_func (T, f->n);
 
-  if (x->size != n) 
+  status = gsl_multiroot_fdfsolver_set (s, f, x); /* seed the generator */
+  
+  if (status != GSL_SUCCESS)
     {
-      GSL_ERROR_VAL ("vector length not compatible with function", 
-                        GSL_EBADLEN, 0);
+      (s->type->free)(s->state);
+      free (s->state);
+      gsl_vector_free (s->dx);
+      gsl_vector_free (s->x);
+      gsl_vector_free (s->f);
+      gsl_matrix_free (s->J);
+      free (s);		/* exception in constructor, avoid memory leak */
+      
+      GSL_ERROR_VAL ("failed to set solver", status, 0);
     }
 
+  return s;
+}
+
+gsl_multiroot_fdfsolver *
+gsl_multiroot_fdfsolver_alloc_no_func (const gsl_multiroot_fdfsolver_type * T, 
+                                       size_t n)
+{
+  int status;
+
+  gsl_multiroot_fdfsolver * s;
+
   s = (gsl_multiroot_fdfsolver *) malloc (sizeof (gsl_multiroot_fdfsolver));
 
   if (s == 0)
@@ -116,21 +136,8 @@
       GSL_ERROR_VAL ("failed to set solver", status, 0);
     }
   
-  status = gsl_multiroot_fdfsolver_set (s, f, x); /* seed the generator */
+  s->fdf = NULL;
   
-  if (status != GSL_SUCCESS)
-    {
-      (s->type->free)(s->state);
-      free (s->state);
-      gsl_vector_free (s->dx);
-      gsl_vector_free (s->x);
-      gsl_vector_free (s->f);
-      gsl_matrix_free (s->J);
-      free (s);		/* exception in constructor, avoid memory leak */
-      
-      GSL_ERROR_VAL ("failed to set solver", status, 0);
-    }
-
   return s;
 }
 
@@ -139,6 +146,17 @@
                              gsl_multiroot_function_fdf * f, 
                              gsl_vector * x)
 {
+  if (s->x->size != f->n)
+  {
+    GSL_ERROR("function incompatible with solver size", GSL_EBADLEN);
+  }
+    
+  if (x->size != f->n) 
+  {
+    GSL_ERROR_VAL ("vector length not compatible with function", 
+                        GSL_EBADLEN, 0);
+  }  
+    
   s->fdf = f;
   gsl_vector_memcpy(s->x,x);
   
@@ -148,6 +166,11 @@
 int
 gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * s)
 {
+  if (!s->fdf)
+  {
+    GSL_ERROR ("fdfsolver not set", GSL_EINVAL);
+  }
+    
   return (s->type->iterate) (s->state, s->fdf, s->x, s->f, s->J, s->dx);
 }
 
diff -ur multiroots.orig/fsolver.c multiroots/fsolver.c
--- multiroots.orig/fsolver.c	Mon Feb 12 16:20:56 2001
+++ multiroots/fsolver.c	Mon Feb 12 16:39:16 2001
@@ -31,14 +31,32 @@
 
   gsl_multiroot_fsolver * s;
 
-  const size_t n = f->n ;
-
-  if (x->size != n) 
+  s = gsl_multiroot_fsolver_alloc_no_func (T, f->n);
+  
+  status = gsl_multiroot_fsolver_set (s, f, x); /* seed the generator */
+  
+  if (status != GSL_SUCCESS)
     {
-      GSL_ERROR_VAL ("vector length not compatible with function", 
-                        GSL_EBADLEN, 0);
+      free (s->state);
+      gsl_vector_free (s->dx);
+      gsl_vector_free (s->x);
+      gsl_vector_free (s->f);
+      free (s);		/* exception in constructor, avoid memory leak */
+      
+      GSL_ERROR_VAL ("failed to set solver", status, 0);
     }
 
+  return s;
+}
+
+gsl_multiroot_fsolver * 
+gsl_multiroot_fsolver_alloc_no_func (const gsl_multiroot_fsolver_type * T,
+                                     size_t n)
+{
+  int status;
+
+  gsl_multiroot_fsolver * s;
+
   s = (gsl_multiroot_fsolver *) malloc (sizeof (gsl_multiroot_fsolver));
 
   if (s == 0)
@@ -102,20 +120,9 @@
       
       GSL_ERROR_VAL ("failed to set solver", status, 0);
     }
-  
-  status = gsl_multiroot_fsolver_set (s, f, x); /* seed the generator */
-  
-  if (status != GSL_SUCCESS)
-    {
-      free (s->state);
-      gsl_vector_free (s->dx);
-      gsl_vector_free (s->x);
-      gsl_vector_free (s->f);
-      free (s);		/* exception in constructor, avoid memory leak */
-      
-      GSL_ERROR_VAL ("failed to set solver", status, 0);
-    }
 
+  s->function = NULL;
+  
   return s;
 }
 
@@ -124,6 +131,17 @@
                            gsl_multiroot_function * f, 
                            gsl_vector * x)
 {
+  if (s->x->size != f->n)
+  {
+    GSL_ERROR("function incompatible with solver size", GSL_EBADLEN);
+  }
+    
+  if (x->size != f->n)
+  {
+    GSL_ERROR_VAL ("vector length not compatible with function", 
+                        GSL_EBADLEN, 0);
+  }  
+    
   s->function = f;
   gsl_vector_memcpy(s->x,x);
   
@@ -133,6 +151,11 @@
 int
 gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * s)
 {
+  if (!s->function)
+  {
+    GSL_ERROR ("fsolver not set", GSL_EINVAL);
+  }
+            
   return (s->type->iterate) (s->state, s->function, s->x, s->f, s->dx);
 }
 
diff -ur multiroots.orig/gsl_multiroots.h multiroots/gsl_multiroots.h
--- multiroots.orig/gsl_multiroots.h	Mon Feb 12 16:20:56 2001
+++ multiroots/gsl_multiroots.h	Mon Feb 12 15:59:10 2001
@@ -80,6 +80,11 @@
 gsl_multiroot_fsolver *
 gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, 
                              gsl_multiroot_function * f, gsl_vector * x);
+
+gsl_multiroot_fsolver * 
+gsl_multiroot_fsolver_alloc_no_func (const gsl_multiroot_fsolver_type * T, 
+                                     size_t n); 
+
 void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s);
 
 int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s, 
@@ -136,6 +141,10 @@
 gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, 
                                gsl_multiroot_function_fdf * fdf,
                                gsl_vector * x);
+
+gsl_multiroot_fdfsolver *
+gsl_multiroot_fdfsolver_alloc_no_func(const gsl_multiroot_fdfsolver_type * T,
+                                      size_t n);
 
 int
 gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s, 

------ 

Toby

-- 
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