/* specfunc/mathieu_workspace.c

 * 

 * Copyright (C) 2003 Lowell Johnson

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

 */

/* Author:  L. Johnson */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_mathieu.h>

gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,

                                               const double qq)

{

  gsl_sf_mathieu_workspace *workspace;

  unsigned int even_order = nn/2 + 1, odd_order = (nn + 1)/2,

               extra_values;

  /* Compute the maximum number of extra terms required for 10^-18 root

     accuracy for a given value of q (contributed by Brian Gladman). */

  extra_values = (int)(2.1*pow(fabs(qq), 0.37)) + 9;

  extra_values += 20;  /* additional fudge */

  if (nn + 1 == 0)

  {

      GSL_ERROR_NULL("matrix dimension must be positive integer", GSL_EINVAL);

  }

  workspace =

         (gsl_sf_mathieu_workspace *)malloc(sizeof(gsl_sf_mathieu_workspace));

  if (workspace == NULL)

  {

      GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);

  }

  /* Extend matrices to ensure accuracy. */

  even_order += extra_values;

  odd_order += extra_values;

  workspace->size = nn;

  workspace->even_order = even_order;

  workspace->odd_order = odd_order;

  workspace->extra_values = extra_values;

  /* Allocate space for the characteristic values. */

  workspace->aa = (double *)malloc((nn+1)*sizeof(double));

  if (workspace->aa == NULL)

  {

      free(workspace);

      GSL_ERROR_NULL("Error allocating memory for characteristic a values",

                     GSL_ENOMEM);

  }

  workspace->bb = (double *)malloc((nn+1)*sizeof(double));

  if (workspace->bb == NULL)

  {

      free(workspace->aa);

      free(workspace);

      GSL_ERROR_NULL("Error allocating memory for characteristic b values",

                     GSL_ENOMEM);

  }

  /* Since even_order is always >= odd_order, dimension the arrays for

     even_order. */

  workspace->dd = (double *)malloc(even_order*sizeof(double));

  if (workspace->dd == NULL)

  {

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);

  }

  workspace->ee = (double *)malloc(even_order*sizeof(double));

  if (workspace->ee == NULL)

  {

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);

  }

  workspace->tt = (double *)malloc(3*even_order*sizeof(double));

  if (workspace->tt == NULL)

  {

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);

  }

  workspace->e2 = (double *)malloc(even_order*sizeof(double));

  if (workspace->e2 == NULL)

  {

      free(workspace->tt);

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);

  }

  workspace->zz = (double *)malloc(even_order*even_order*sizeof(double));

  if (workspace->zz == NULL)

  {

      free(workspace->e2);

      free(workspace->tt);

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for diagonal", GSL_ENOMEM);

  }

  workspace->eval = gsl_vector_alloc(even_order);

  if (workspace->eval == NULL)

    {

      free(workspace->zz);

      free(workspace->e2);

      free(workspace->tt);

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for eval", GSL_ENOMEM);

    }

  workspace->evec = gsl_matrix_alloc(even_order, even_order);

  if (workspace->evec == NULL)

    {

      gsl_vector_free (workspace->eval);

      free(workspace->zz);

      free(workspace->e2);

      free(workspace->tt);

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for evec", GSL_ENOMEM);

    }

  workspace->wmat = gsl_eigen_symmv_alloc(even_order);

  if (workspace->wmat == NULL)

    {

      gsl_matrix_free (workspace->evec);

      gsl_vector_free (workspace->eval);

      free(workspace->zz);

      free(workspace->e2);

      free(workspace->tt);

      free(workspace->ee);

      free(workspace->dd);

      free(workspace->aa);

      free(workspace->bb);

      free(workspace);

      GSL_ERROR_NULL("failed to allocate space for wmat", GSL_ENOMEM);

    }

  return workspace;

}

void gsl_sf_mathieu_free(gsl_sf_mathieu_workspace *workspace)

{

  RETURN_IF_NULL (workspace);

  gsl_vector_free(workspace->eval);

  gsl_matrix_free(workspace->evec);

  gsl_eigen_symmv_free(workspace->wmat);

  free(workspace->aa);

  free(workspace->bb);

  free(workspace->dd);

  free(workspace->ee);

  free(workspace->tt);

  free(workspace->e2);

  free(workspace->zz);

  free(workspace);

}