This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: Additional multiroots functions.


Brian Gough <bjg@network-theory.co.uk> writes:

> It probably makes sense for alloc and set to become orthogonal, so
> that alloc only allocates memory, which would also solve your problem
> I think.

Please find below a patch which does this. I've also updated test.c
to reflect the change in semantics, and adjusted doc/multiroots.texi
to correctly describe the functions.

Thanks,
Toby

diff -ur gsl.orig/doc/multiroots.texi gsl/doc/multiroots.texi
--- gsl.orig/doc/multiroots.texi	Sun May 14 21:10:13 2000
+++ gsl/doc/multiroots.texi	Wed Feb 14 12:15:16 2001
@@ -94,15 +94,15 @@
 @node Initializing the Multidimensional Solver
 @section Initializing the Solver
 
-@deftypefun {gsl_multiroot_fsolver *} gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * @var{T}, gsl_multiroot_function * @var{f}, gsl_vector * @var{x})
+@deftypefun {gsl_multiroot_fsolver *} gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * @var{T}, size_t @var{n})
 This function returns a pointer to a a newly allocated instance of a
-solver of type @var{T} for the function @var{f}, with an initial bracket
-on the root of @var{x}.  For example, the following code creates an
-instance of a Brent solver,
+solver of type @var{T} for a system of @var{n} dimensions.
+For example, the following code creates an instance of a Brent solver, 
+to solve a 3-dimensional system of equations.
 
 @example
 gsl_multiroot_fsolver * s = 
-    gsl_multiroot_fsolver_alloc (gsl_multiroot_fsolver_brent, f, x);
+    gsl_multiroot_fsolver_alloc (gsl_multiroot_fsolver_brent, 3);
 @end example
 
 If there is insufficient memory to create the solver then the function
@@ -110,15 +110,15 @@
 code of @code{GSL_ENOMEM}.
 @end deftypefun
 
-@deftypefun {gsl_multiroot_fdfsolver *} gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * @var{T}, gsl_multiroot_function_fdf * @var{fdf}, gsl_vector * @var{x})
+@deftypefun {gsl_multiroot_fdfsolver *} gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * @var{T}, size_t @var{n})
 This function returns a pointer to a a newly allocated instance of a
-derivative solver of type @var{T} for the function @var{fdf}, with an
-initial guess for the root of @var{x}.  For example, the following code
-creates an instance of a Newton-Raphson solver,
+derivative solver of type @var{T} for a system of size_t @var{n} dimensions.
+For example, the following code creates an instance of a Newton-Raphson solver,
+for a 2-dimensional system of equations.
 
 @example
 gsl_multiroot_fdfsolver * s = 
-    gsl_multiroot_fdfsolver_alloc (gsl_multiroot_fdfsolver_newton, fdf, x);
+    gsl_multiroot_fdfsolver_alloc (gsl_multiroot_fdfsolver_newton, 2);
 @end example
 
 If there is insufficient memory to create the solver then the function
@@ -127,12 +127,12 @@
 @end deftypefun
 
 @deftypefun int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * @var{s}, gsl_multiroot_function * @var{f}, gsl_vector * @var{x})
-This function reinitializes an existing solver @var{s} to use the
+This function sets, or resets, an existing solver @var{s} to use the
 function @var{f} and the initial guess @var{x}.
 @end deftypefun
 
 @deftypefun int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * @var{s}, gsl_function_fdf * @var{fdf}, gsl_vector * @var{x})
-This function reinitializes an existing solver @var{s} to use the
+This function sets, or resets, an existing solver @var{s} to use the
 function and derivative @var{fdf} and the initial guess @var{x}.
 @end deftypefun
 
@@ -778,7 +778,8 @@
   gsl_vector_set (x, 1, x_init[1]);
 
   T = gsl_multiroot_fsolver_hybrids;
-  s = gsl_multiroot_fsolver_alloc (T, &f, x);
+  s = gsl_multiroot_fsolver_alloc (T, 2);
+  gsl_multiroot_fsolver_set (s, &f, x);
 
   print_state (iter, s);
 
@@ -919,7 +920,8 @@
   gsl_vector_set (x, 1, x_init[1]);
 
   T = gsl_multiroot_fdfsolver_gnewton;
-  s = gsl_multiroot_fdfsolver_alloc (T, &f, x);
+  s = gsl_multiroot_fdfsolver_alloc (T, n);
+  gsl_multiroot_fdfsolver_alloc (s, &f, x);
 
   print_state (iter, s);
 
diff -ur gsl.orig/multiroots/fdfsolver.c gsl/multiroots/fdfsolver.c
--- gsl.orig/multiroots/fdfsolver.c	Mon Jun  5 20:27:37 2000
+++ gsl/multiroots/fdfsolver.c	Wed Feb 14 12:07:12 2001
@@ -25,21 +25,12 @@
 
 gsl_multiroot_fdfsolver *
 gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, 
-                               gsl_multiroot_function_fdf * f, 
-                               gsl_vector * x)
+                                       size_t n)
 {
   int status;
 
   gsl_multiroot_fdfsolver * s;
 
-  const size_t n = f->n ;
-
-  if (x->size != n) 
-    {
-      GSL_ERROR_VAL ("vector length not compatible with function", 
-                        GSL_EBADLEN, 0);
-    }
-
   s = (gsl_multiroot_fdfsolver *) malloc (sizeof (gsl_multiroot_fdfsolver));
 
   if (s == 0)
@@ -116,21 +107,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 +117,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 +137,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 gsl.orig/multiroots/fsolver.c gsl/multiroots/fsolver.c
--- gsl.orig/multiroots/fsolver.c	Mon Jun  5 20:27:37 2000
+++ gsl/multiroots/fsolver.c	Wed Feb 14 12:06:49 2001
@@ -23,22 +23,14 @@
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_multiroots.h>
 
-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 (const gsl_multiroot_fsolver_type * T,
+                                     size_t n)
 {
   int status;
 
   gsl_multiroot_fsolver * s;
 
-  const size_t n = f->n ;
-
-  if (x->size != n) 
-    {
-      GSL_ERROR_VAL ("vector length not compatible with function", 
-                        GSL_EBADLEN, 0);
-    }
-
   s = (gsl_multiroot_fsolver *) malloc (sizeof (gsl_multiroot_fsolver));
 
   if (s == 0)
@@ -102,20 +94,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 +105,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 +125,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 gsl.orig/multiroots/gsl_multiroots.h gsl/multiroots/gsl_multiroots.h
--- gsl.orig/multiroots/gsl_multiroots.h	Thu May  4 12:25:05 2000
+++ gsl/multiroots/gsl_multiroots.h	Wed Feb 14 12:05:53 2001
@@ -77,9 +77,10 @@
   }
 gsl_multiroot_fsolver;
 
-gsl_multiroot_fsolver *
+gsl_multiroot_fsolver * 
 gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, 
-                             gsl_multiroot_function * f, gsl_vector * x);
+                                     size_t n); 
+
 void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s);
 
 int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s, 
@@ -131,11 +132,9 @@
   }
 gsl_multiroot_fdfsolver;
 
-
 gsl_multiroot_fdfsolver *
-gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, 
-                               gsl_multiroot_function_fdf * fdf,
-                               gsl_vector * x);
+gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T,
+                                      size_t n);
 
 int
 gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s, 
diff -ur gsl.orig/multiroots/test.c gsl/multiroots/test.c
--- gsl.orig/multiroots/test.c	Thu May  4 12:25:05 2000
+++ gsl/multiroots/test.c	Wed Feb 14 12:07:58 2001
@@ -134,7 +134,8 @@
 
   if (factor != 1.0) scale(x, factor);
 
-  s = gsl_multiroot_fdfsolver_alloc (T, function, x);
+  s = gsl_multiroot_fdfsolver_alloc (T, n);
+  gsl_multiroot_fdfsolver_set (s, function, x);
  
   do
     {
@@ -220,7 +221,8 @@
 
   if (factor != 1.0) scale(x, factor);
 
-  s = gsl_multiroot_fsolver_alloc (T, &function, x);
+  s = gsl_multiroot_fsolver_alloc (T, n);
+  gsl_multiroot_fsolver_set (s, &function, x);
 
 /*   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); */
 /*   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); */


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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]