Re: multidimensional root-finding problem :-(

["Peter S. Christopher" <peterc at midway dot uchicago dot edu>]

Hi there,

	I've had experience with the Multi-Min package. My best guess is
that some of the routines use ONLY fdf and some of them use f, df and fdf.
What I do to check is to put asserts in my f,df functions -- to make sure
that they aren't being called. So, for example, replace resenbrock{_f,_df}
with the following
> > int
> > rosenbrock_f (const gsl_vector * x, void *params,
> >               gsl_vector * f)
> > {
      assert(0);
> >   return GSL_SUCCESS;
> > }
> >
> > int
> > rosenbrock_df (const gsl_vector * x, void *params,
> >                gsl_matrix * J)
> > {
      assert(0);
> >   return GSL_SUCCESS;
> > }

If these functions are called you'll get termination and you'll be able to
tell that these functions were called. The initial x is probably remaining
the same because these functions are returning GSL_SUCCESS and the
algorithm thinks that they are doing something. 

Best of luck,
Pete

On Fri, 23 May 2003, Ivo Alxneit wrote:

> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
> 
> On Friday 23 May 2003 04:34, Slaven Peles wrote:
> > Hi,
> >
> > I was working with the example in GSL Reference Manual for multidimensional
> > root finding using Newton's method:
> > http://sources.redhat.com/gsl/ref/gsl-ref_33.html#SEC445
> > I wanted to have both, equations and the Jacobi matrix, in the same
> > function in order to optimize calculation, so I left functions rosenbrock_f
> > and rosenbrock_df empty, and put everything into rosenbrock_fdf (see
> > below). According to the Manual this should work, however when I run this
> > program my initial value for x remains unchanged with every iteration. If
> > anyone on the list has experience with GSL multidimensional root-finding
> > functions please help. My code is below.
> >
> > Many thanks,
> > Slaven
> >
> > /**************************************************************************
> >******/ #include <stdlib.h>
> > #include <stdio.h>
> > #include <gsl/gsl_vector.h>
> > #include <gsl/gsl_multiroots.h>
> >
> > struct rparams
> >   {
> >     double a;
> >     double b;
> >   };
> >
> > int
> > rosenbrock_f (const gsl_vector * x, void *params,
> >               gsl_vector * f)
> > {
> >   return GSL_SUCCESS;
> > }
> >
> > int
> > rosenbrock_df (const gsl_vector * x, void *params,
> >                gsl_matrix * J)
> > {
> >   return GSL_SUCCESS;
> > }
> >
> > int
> > rosenbrock_fdf (const gsl_vector * x, void *params,
> >                 gsl_vector * f, gsl_matrix * J)
> > {
> >   double a = ((struct rparams *) params)->a;
> >   double b = ((struct rparams *) params)->b;
> >
> >   double x0 = gsl_vector_get (x, 0);
> >   double x1 = gsl_vector_get (x, 1);
> >
> >   double y0 = a * (1 - x0);
> >   double y1 = b * (x1 - x0 * x0);
> >
> >   gsl_vector_set (f, 0, y0);
> >   gsl_vector_set (f, 1, y1);
> >
> >   gsl_matrix_set (J, 0, 0, -a);
> >   gsl_matrix_set (J, 0, 1, 0);
> >   gsl_matrix_set (J, 1, 0,  -2 * b  * x0);
> >   gsl_matrix_set (J, 1, 1, b);
> >
> >   return GSL_SUCCESS;
> > }
> >
> >
> > int
> > print_state (size_t iter, gsl_multiroot_fdfsolver * s)
> > {
> >   printf ("iter = %3u x = % .3f % .3f "
> >           "f(x) = % .3e % .3e\n",
> >           iter,
> >           gsl_vector_get (s->x, 0),
> >           gsl_vector_get (s->x, 1),
> >           gsl_vector_get (s->f, 0),
> >           gsl_vector_get (s->f, 1));
> > }
> >
> >
> > int
> > main (void)
> > {
> >   const gsl_multiroot_fdfsolver_type *T;
> >   gsl_multiroot_fdfsolver *s;
> >
> >   int status;
> >   size_t iter = 0;
> >
> >   const size_t n = 2;
> >   struct rparams p = {1.0, 10.0};
> >   gsl_multiroot_function_fdf f = {&rosenbrock_f,
> >                                   &rosenbrock_df,
> >                                   &rosenbrock_fdf,
> >                                   n, &p};
> >
> >   double x_init[2] = {-10.0, -5.0};
> >   gsl_vector *x = gsl_vector_alloc (n);
> >
> >   gsl_vector_set (x, 0, x_init[0]);
> >   gsl_vector_set (x, 1, x_init[1]);
> >
> >   T = gsl_multiroot_fdfsolver_gnewton;
> >   s = gsl_multiroot_fdfsolver_alloc (T, n);
> >   gsl_multiroot_fdfsolver_set (s, &f, x);
> >
> >   print_state (iter, s);
> >
> >   do
> >     {
> >       iter++;
> >
> >       status = gsl_multiroot_fdfsolver_iterate (s);
> >
> >       print_state (iter, s);
> >
> >       if (status)
> >         break;
> >
> >       status = gsl_multiroot_test_residual (s->f, 1e-7);
> >     }
> >   while (status == GSL_CONTINUE && iter < 1000);
> >
> >   printf ("status = %s\n", gsl_strerror (status));
> >
> >   gsl_multiroot_fdfsolver_free (s);
> >   gsl_vector_free (x);
> >   return 0;
> > }
> i think there is a missconception here
> you need all three (f, df, and fdf) perform what they are supposed to do. your 
> iteration fails because your f is just a dummy function.
> the "additional" function fdf is used in cases where both f and df are needed 
> at the same time. you often can code it to be much faster that separate calls 
> to f and df. so the way it is in the manual is actually correct. 
> - -- 
> Dr. Ivo Alxneit
> Laboratory for Solar Technology   phone: +41 56 310 4092
> Paul Scherrer Institute             fax: +41 56 310 2624
> CH-5232 Villigen                   http://solar.web.psi.ch
> Switzerland                   gnupg key: 0x515E30C7
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.2.0 (GNU/Linux)
> 
> iD8DBQE+zifeAd7CE1FeMMcRAiS6AJ9+ZPXznDOuAgmKQh1FCO4Jj7+NLACgrn17
> FK6c5DE++DVjAA/BdKNFkbI=
> =1Q/q
> -----END PGP SIGNATURE-----
> 
>