htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/mathieu_workspace.c

   1 /* specfunc/mathieu_workspace.c
   2  *
   3  * Copyright (C) 2003 Lowell Johnson
   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., 675 Mass Ave, Cambridge, MA 02139, USA.
  18  */
  19
  20 /* Author:  L. Johnson */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_mathieu.h>
  27
  28
  29 gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,
  30                                                const double qq)
  31 {
  32   gsl_sf_mathieu_workspace *workspace;
  33   unsigned int even_order = nn/2 + 1, odd_order = (nn + 1)/2,
  34                extra_values;
  35
  36   /* Compute the maximum number of extra terms required for 10^-18 root
  37      accuracy for a given value of q (contributed by Brian Gladman). */
  38   extra_values = (int)(2.1*pow(fabs(qq), 0.37)) + 9;
  39
  40   if (nn + 1 == 0)
  41   {
  42       GSL_ERROR_NULL("matrix dimension must be positive integer", GSL_EINVAL);
  43   }
  44
  45   workspace =
  46          (gsl_sf_mathieu_workspace *)malloc(sizeof(gsl_sf_mathieu_workspace));
  47   if (workspace == NULL)
  48   {
  49       GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);
  50   }
  51
  52   /* Extend matrices to ensure accuracy. */
  53   even_order += extra_values;
  54   odd_order += extra_values;
  55
  56   workspace->size = nn;
  57   workspace->even_order = even_order;
  58   workspace->odd_order = odd_order;
  59   workspace->extra_values = extra_values;
  60
  61   /* Allocate space for the characteristic values. */
  62   workspace->aa = (double *)malloc((nn+1)*sizeof(double));
  63   if (workspace->aa == NULL)
  64   {
  65       free(workspace);
  66       GSL_ERROR_NULL("Error allocating memory for characteristic a values",
  67                      GSL_ENOMEM);
  68   }
  69
  70   workspace->bb = (double *)malloc((nn+1)*sizeof(double));
  71   if (workspace->bb == NULL)
  72   {
  73       free(workspace->aa);
  74       free(workspace);
  75       GSL_ERROR_NULL("Error allocating memory for characteristic b values",
  76                      GSL_ENOMEM);
  77   }
  78
  79   /* Since even_order is always >= odd_order, dimension the arrays for
  80      even_order. */
  81
  82   workspace->dd = (double *)malloc(even_order*sizeof(double));
  83   if (workspace->dd == NULL)
  84   {
  85       free(workspace->aa);
  86       free(workspace->bb);
  87       free(workspace);
  88       GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);
  89   }
  90
  91   workspace->ee = (double *)malloc(even_order*sizeof(double));
  92   if (workspace->ee == NULL)
  93   {
  94       free(workspace->dd);
  95       free(workspace->aa);
  96       free(workspace->bb);
  97       free(workspace);
  98       GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);
  99   }
 100
 101   workspace->tt = (double *)malloc(3*even_order*sizeof(double));
 102   if (workspace->tt == NULL)
 103   {
 104       free(workspace->ee);
 105       free(workspace->dd);
 106       free(workspace->aa);
 107       free(workspace->bb);
 108       free(workspace);
 109       GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);
 110   }
 111
 112   workspace->e2 = (double *)malloc(even_order*sizeof(double));
 113   if (workspace->e2 == NULL)
 114   {
 115       free(workspace->tt);
 116       free(workspace->ee);
 117       free(workspace->dd);
 118       free(workspace->aa);
 119       free(workspace->bb);
 120       free(workspace);
 121       GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);
 122   }
 123
 124   workspace->zz = (double *)malloc(even_order*even_order*sizeof(double));
 125   if (workspace->zz == NULL)
 126   {
 127       free(workspace->e2);
 128       free(workspace->tt);
 129       free(workspace->ee);
 130       free(workspace->dd);
 131       free(workspace->aa);
 132       free(workspace->bb);
 133       free(workspace);
 134       GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);
 135   }
 136
 137   workspace->eval = gsl_vector_alloc(even_order);
 138
 139   if (workspace->eval == NULL)
 140     {
 141       free(workspace->zz);
 142       free(workspace->e2);
 143       free(workspace->tt);
 144       free(workspace->ee);
 145       free(workspace->dd);
 146       free(workspace->aa);
 147       free(workspace->bb);
 148       free(workspace);
 149       GSL_ERROR_NULL("failed to allocate space for eval", GSL_ENOMEM);
 150     }
 151
 152   workspace->evec = gsl_matrix_alloc(even_order, even_order);
 153
 154   if (workspace->evec == NULL)
 155     {
 156       gsl_vector_free (workspace->eval);
 157       free(workspace->zz);
 158       free(workspace->e2);
 159       free(workspace->tt);
 160       free(workspace->ee);
 161       free(workspace->dd);
 162       free(workspace->aa);
 163       free(workspace->bb);
 164       free(workspace);
 165       GSL_ERROR_NULL("failed to allocate space for evec", GSL_ENOMEM);
 166     }
 167
 168   workspace->wmat = gsl_eigen_symmv_alloc(even_order);
 169
 170   if (workspace->wmat == NULL)
 171     {
 172       gsl_matrix_free (workspace->evec);
 173       gsl_vector_free (workspace->eval);
 174       free(workspace->zz);
 175       free(workspace->e2);
 176       free(workspace->tt);
 177       free(workspace->ee);
 178       free(workspace->dd);
 179       free(workspace->aa);
 180       free(workspace->bb);
 181       free(workspace);
 182       GSL_ERROR_NULL("failed to allocate space for wmat", GSL_ENOMEM);
 183     }
 184
 185   return workspace;
 186 }
 187
 188
 189 void gsl_sf_mathieu_free(gsl_sf_mathieu_workspace *workspace)
 190 {
 191   gsl_vector_free(workspace->eval);
 192   gsl_matrix_free(workspace->evec);
 193   gsl_eigen_symmv_free(workspace->wmat);
 194   free(workspace->aa);
 195   free(workspace->bb);
 196   free(workspace->dd);
 197   free(workspace->ee);
 198   free(workspace->tt);
 199   free(workspace->e2);
 200   free(workspace->zz);
 201   free(workspace);
 202 }