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

   1 /* multiroots/fdjac.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 <gsl/gsl_multiroots.h>
  22
  23 int
  24 gsl_multiroot_fdjacobian (gsl_multiroot_function * F,
  25                            const gsl_vector * x, const gsl_vector * f,
  26                            double epsrel, gsl_matrix * jacobian)
  27 {
  28   const size_t n = x->size;
  29   const size_t m = f->size;
  30   const size_t n1 = jacobian->size1;
  31   const size_t n2 = jacobian->size2;
  32   int status = 0;
  33
  34   if (m != n1 || n != n2)
  35     {
  36       GSL_ERROR ("function and jacobian are not conformant", GSL_EBADLEN);
  37     }
  38
  39   {
  40     size_t i,j;
  41     gsl_vector *x1, *f1;
  42
  43     x1 = gsl_vector_alloc (n);
  44
  45     if (x1 == 0)
  46       {
  47         GSL_ERROR ("failed to allocate space for x1 workspace", GSL_ENOMEM);
  48       }
  49
  50     f1 = gsl_vector_alloc (m);
  51
  52     if (f1 == 0)
  53       {
  54         gsl_vector_free (x1);
  55
  56         GSL_ERROR ("failed to allocate space for f1 workspace", GSL_ENOMEM);
  57       }
  58
  59     gsl_vector_memcpy (x1, x);  /* copy x into x1 */
  60
  61     for (j = 0; j < n; j++)
  62       {
  63         double xj = gsl_vector_get (x, j);
  64         double dx = epsrel * fabs (xj);
  65
  66         if (dx == 0)
  67           {
  68             dx = epsrel;
  69           }
  70
  71         gsl_vector_set (x1, j, xj + dx);
  72
  73         {
  74           int f_stat = GSL_MULTIROOT_FN_EVAL (F, x1, f1);
  75
  76           if (f_stat != GSL_SUCCESS)
  77             {
  78               status = GSL_EBADFUNC;
  79               break; /* n.b. avoid memory leak for x1,f1 */
  80             }
  81         }
  82
  83         gsl_vector_set (x1, j, xj);
  84
  85         for (i = 0; i < m; i++)
  86           {
  87             double g1 = gsl_vector_get (f1, i);
  88             double g0 = gsl_vector_get (f, i);
  89             gsl_matrix_set (jacobian, i, j, (g1 - g0) / dx);
  90           }
  91
  92         {
  93           gsl_vector_view col = gsl_matrix_column (jacobian, j);
  94           int null_col = gsl_vector_isnull (&col.vector);
  95           /* if column is null, return an error - this may be due to
  96              dx being too small. Try increasing epsrel */
  97           if (null_col) {
  98             status = GSL_ESING;
  99           }
 100         }
 101       }
 102
 103     gsl_vector_free (x1);
 104     gsl_vector_free (f1);
 105   }
 106
 107   if (status)
 108     return status;
 109   else
 110     return GSL_SUCCESS;
 111 }