htswanalysis/MACS/lib/gsl/gsl-1.11/multiroots/test.c

   1 /* multiroots/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <stdio.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_multiroots.h>
  26
  27 #include <gsl/gsl_ieee_utils.h>
  28
  29 #include "test_funcs.h"
  30 int test_fdf (const char * desc, gsl_multiroot_function_fdf * function, initpt_function initpt, double factor, const gsl_multiroot_fdfsolver_type * T);
  31 int test_f (const char * desc, gsl_multiroot_function_fdf * fdf, initpt_function initpt, double factor, const gsl_multiroot_fsolver_type * T);
  32
  33
  34 int
  35 main (void)
  36 {
  37   const gsl_multiroot_fsolver_type * fsolvers[5] ;
  38   const gsl_multiroot_fsolver_type ** T1 ;
  39
  40   const gsl_multiroot_fdfsolver_type * fdfsolvers[5] ;
  41   const gsl_multiroot_fdfsolver_type ** T2 ;
  42
  43   double f;
  44
  45   fsolvers[0] = gsl_multiroot_fsolver_dnewton;
  46   fsolvers[1] = gsl_multiroot_fsolver_broyden;
  47   fsolvers[2] = gsl_multiroot_fsolver_hybrid;
  48   fsolvers[3] = gsl_multiroot_fsolver_hybrids;
  49   fsolvers[4] = 0;
  50
  51   fdfsolvers[0] = gsl_multiroot_fdfsolver_newton;
  52   fdfsolvers[1] = gsl_multiroot_fdfsolver_gnewton;
  53   fdfsolvers[2] = gsl_multiroot_fdfsolver_hybridj;
  54   fdfsolvers[3] = gsl_multiroot_fdfsolver_hybridsj;
  55   fdfsolvers[4] = 0;
  56
  57   gsl_ieee_env_setup();
  58
  59
  60   f = 1.0 ;
  61
  62   T1 = fsolvers ;
  63
  64   while (*T1 != 0)
  65     {
  66       test_f ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T1);
  67       test_f ("Roth", &roth, roth_initpt, f, *T1);
  68       test_f ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T1);
  69       test_f ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T1);
  70       test_f ("Powell singular", &powellsing, powellsing_initpt, f, *T1);
  71       test_f ("Wood", &wood, wood_initpt, f, *T1);
  72       test_f ("Helical", &helical, helical_initpt, f, *T1);
  73       test_f ("Discrete BVP", &dbv, dbv_initpt, f, *T1);
  74       test_f ("Trig", &trig, trig_initpt, f, *T1);
  75       T1++;
  76     }
  77
  78   T2 = fdfsolvers ;
  79
  80   while (*T2 != 0)
  81     {
  82       test_fdf ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T2);
  83       test_fdf ("Roth", &roth, roth_initpt, f, *T2);
  84       test_fdf ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T2);
  85       test_fdf ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T2);
  86       test_fdf ("Powell singular", &powellsing, powellsing_initpt, f, *T2);
  87       test_fdf ("Wood", &wood, wood_initpt, f, *T2);
  88       test_fdf ("Helical", &helical, helical_initpt, f, *T2);
  89       test_fdf ("Discrete BVP", &dbv, dbv_initpt, f, *T2);
  90       test_fdf ("Trig", &trig, trig_initpt, f, *T2);
  91       T2++;
  92     }
  93
  94   exit (gsl_test_summary ());
  95 }
  96
  97 void scale (gsl_vector * x, double factor);
  98
  99 void
 100 scale (gsl_vector * x, double factor)
 101 {
 102   size_t i, n = x->size;
 103
 104   if (gsl_vector_isnull(x))
 105     {
 106       for (i = 0; i < n; i++)
 107         {
 108           gsl_vector_set (x, i, factor);
 109         }
 110     }
 111   else
 112     {
 113       for (i = 0; i < n; i++)
 114         {
 115           double xi = gsl_vector_get(x, i);
 116           gsl_vector_set(x, i, factor * xi);
 117         }
 118     }
 119 }
 120
 121 int
 122 test_fdf (const char * desc, gsl_multiroot_function_fdf * function,
 123           initpt_function initpt, double factor,
 124           const gsl_multiroot_fdfsolver_type * T)
 125 {
 126   int status;
 127   double residual = 0;
 128   size_t i, n = function->n, iter = 0;
 129
 130   gsl_vector *x = gsl_vector_alloc (n);
 131   gsl_matrix *J = gsl_matrix_alloc (n, n);
 132
 133   gsl_multiroot_fdfsolver *s;
 134
 135   (*initpt) (x);
 136
 137   if (factor != 1.0) scale(x, factor);
 138
 139   s = gsl_multiroot_fdfsolver_alloc (T, n);
 140   gsl_multiroot_fdfsolver_set (s, function, x);
 141
 142   do
 143     {
 144       iter++;
 145       status = gsl_multiroot_fdfsolver_iterate (s);
 146
 147       if (status)
 148         break ;
 149
 150       status = gsl_multiroot_test_residual (s->f, 0.0000001);
 151     }
 152   while (status == GSL_CONTINUE && iter < 1000);
 153
 154 #ifdef DEBUG
 155   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
 156   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
 157 #endif
 158
 159
 160 #ifdef TEST_JACOBIAN
 161  {
 162     double r,sum; size_t j;
 163
 164     gsl_multiroot_function f1 ;
 165     f1.f = function->f ;
 166     f1.n = function->n ;
 167     f1.params = function->params ;
 168
 169     gsl_multiroot_fdjacobian (&f1, s->x, s->f, GSL_SQRT_DBL_EPSILON, J);
 170
 171     /* compare J and s->J */
 172
 173     r=0;sum=0;
 174     for (i = 0; i < n; i++)
 175       for (j = 0; j< n ; j++)
 176         {
 177           double u = gsl_matrix_get(J,i,j);
 178           double su = gsl_matrix_get(s->J, i, j);
 179           r = fabs(u - su)/(1e-6 + 1e-6 * fabs(u)); sum+=r;
 180           if (fabs(u - su) > 1e-6 + 1e-6 * fabs(u))
 181             printf("broken jacobian %g\n", r);
 182         }
 183     printf("avg r = %g\n", sum/(n*n));
 184   }
 185 #endif
 186
 187   for (i = 0; i < n ; i++)
 188     {
 189       residual += fabs(gsl_vector_get(s->f, i));
 190     }
 191
 192   gsl_multiroot_fdfsolver_free (s);
 193   gsl_matrix_free(J);
 194   gsl_vector_free(x);
 195
 196   gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);
 197
 198   return status;
 199 }
 200
 201
 202 int
 203 test_f (const char * desc, gsl_multiroot_function_fdf * fdf,
 204         initpt_function initpt, double factor,
 205         const gsl_multiroot_fsolver_type * T)
 206 {
 207   int status;
 208   size_t i, n = fdf->n, iter = 0;
 209   double residual = 0;
 210
 211   gsl_vector *x;
 212
 213   gsl_multiroot_fsolver *s;
 214   gsl_multiroot_function function;
 215
 216   function.f = fdf->f;
 217   function.params = fdf->params;
 218   function.n = n ;
 219
 220   x = gsl_vector_alloc (n);
 221
 222   (*initpt) (x);
 223
 224   if (factor != 1.0) scale(x, factor);
 225
 226   s = gsl_multiroot_fsolver_alloc (T, n);
 227   gsl_multiroot_fsolver_set (s, &function, x);
 228
 229 /*   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); */
 230 /*   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); */
 231
 232   do
 233     {
 234       iter++;
 235       status = gsl_multiroot_fsolver_iterate (s);
 236
 237       if (status)
 238         break ;
 239
 240       status = gsl_multiroot_test_residual (s->f, 0.0000001);
 241     }
 242   while (status == GSL_CONTINUE && iter < 1000);
 243
 244 #ifdef DEBUG
 245   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n");
 246   printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n");
 247 #endif
 248
 249   for (i = 0; i < n ; i++)
 250     {
 251       residual += fabs(gsl_vector_get(s->f, i));
 252     }
 253
 254   gsl_multiroot_fsolver_free (s);
 255   gsl_vector_free(x);
 256
 257   gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual);
 258
 259   return status;
 260 }