text:
/* blas/blas.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2009 Gerard Jungman & Brian
 * Gough
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* GSL implementation of BLAS operations for vectors and dense
 * matrices.  Note that GSL native storage is row-major.  */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_blas.h>

/* ========================================================================
 * Level 1
 * ========================================================================
 */

/* CBLAS defines vector sizes in terms of int. GSL defines sizes in
   terms of size_t, so we need to convert these into integers.  There
   is the possibility of overflow here. FIXME: Maybe this could be
   caught */

#define INT(X) ((int)(X))

int
gsl_blas_sdsdot (float alpha, const gsl_vector_float * X,
                 const gsl_vector_float * Y, float *result)
{
  if (X->size == Y->size)
    {
      *result =
        cblas_sdsdot (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
                      INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dsdot (const gsl_vector_float * X, const gsl_vector_float * Y,
                double *result)
{
  if (X->size == Y->size)
    {
      *result =
        cblas_dsdot (INT (X->size), X->data, INT (X->stride), Y->data,
                     INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_sdot (const gsl_vector_float * X, const gsl_vector_float * Y,
               float *result)
{
  if (X->size == Y->size)
    {
      *result =
        cblas_sdot (INT (X->size), X->data, INT (X->stride), Y->data,
                    INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result)
{
  if (X->size == Y->size)
    {
      *result =
        cblas_ddot (INT (X->size), X->data, INT (X->stride), Y->data,
                    INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_cdotu (const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y, gsl_complex_float * dotu)
{
  if (X->size == Y->size)
    {
      cblas_cdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
                       INT (Y->stride), GSL_COMPLEX_P (dotu));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_cdotc (const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y, gsl_complex_float * dotc)
{
  if (X->size == Y->size)
    {
      cblas_cdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
                       INT (Y->stride), GSL_COMPLEX_P (dotc));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zdotu (const gsl_vector_complex * X, const gsl_vector_complex * Y,
                gsl_complex * dotu)
{
  if (X->size == Y->size)
    {
      cblas_zdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
                       INT (Y->stride), GSL_COMPLEX_P (dotu));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zdotc (const gsl_vector_complex * X, const gsl_vector_complex * Y,
                gsl_complex * dotc)
{
  if (X->size == Y->size)
    {
      cblas_zdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
                       INT (Y->stride), GSL_COMPLEX_P (dotc));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Norms of vectors */

float
gsl_blas_snrm2 (const gsl_vector_float * X)
{
  return cblas_snrm2 (INT (X->size), X->data, INT (X->stride));
}

double
gsl_blas_dnrm2 (const gsl_vector * X)
{
  return cblas_dnrm2 (INT (X->size), X->data, INT (X->stride));
}

float
gsl_blas_scnrm2 (const gsl_vector_complex_float * X)
{
  return cblas_scnrm2 (INT (X->size), X->data, INT (X->stride));
}

double
gsl_blas_dznrm2 (const gsl_vector_complex * X)
{
  return cblas_dznrm2 (INT (X->size), X->data, INT (X->stride));
}

/* Absolute sums of vectors */

float
gsl_blas_sasum (const gsl_vector_float * X)
{
  return cblas_sasum (INT (X->size), X->data, INT (X->stride));
}

double
gsl_blas_dasum (const gsl_vector * X)
{
  return cblas_dasum (INT (X->size), X->data, INT (X->stride));
}

float
gsl_blas_scasum (const gsl_vector_complex_float * X)
{
  return cblas_scasum (INT (X->size), X->data, INT (X->stride));
}

double
gsl_blas_dzasum (const gsl_vector_complex * X)
{
  return cblas_dzasum (INT (X->size), X->data, INT (X->stride));
}

/* Maximum elements of vectors */

CBLAS_INDEX_t
gsl_blas_isamax (const gsl_vector_float * X)
{
  return cblas_isamax (INT (X->size), X->data, INT (X->stride));
}

CBLAS_INDEX_t
gsl_blas_idamax (const gsl_vector * X)
{
  return cblas_idamax (INT (X->size), X->data, INT (X->stride));
}

CBLAS_INDEX_t
gsl_blas_icamax (const gsl_vector_complex_float * X)
{
  return cblas_icamax (INT (X->size), X->data, INT (X->stride));
}

CBLAS_INDEX_t
gsl_blas_izamax (const gsl_vector_complex * X)
{
  return cblas_izamax (INT (X->size), X->data, INT (X->stride));
}

/* Swap vectors */

int
gsl_blas_sswap (gsl_vector_float * X, gsl_vector_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_sswap (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dswap (gsl_vector * X, gsl_vector * Y)
{
  if (X->size == Y->size)
    {
      cblas_dswap (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    };
}

int
gsl_blas_cswap (gsl_vector_complex_float * X, gsl_vector_complex_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_cswap (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zswap (gsl_vector_complex * X, gsl_vector_complex * Y)
{
  if (X->size == Y->size)
    {
      cblas_zswap (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Copy vectors */

int
gsl_blas_scopy (const gsl_vector_float * X, gsl_vector_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_scopy (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y)
{
  if (X->size == Y->size)
    {
      cblas_dcopy (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ccopy (const gsl_vector_complex_float * X,
                gsl_vector_complex_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_ccopy (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zcopy (const gsl_vector_complex * X, gsl_vector_complex * Y)
{
  if (X->size == Y->size)
    {
      cblas_zcopy (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Compute Y = alpha X + Y */

int
gsl_blas_saxpy (float alpha, const gsl_vector_float * X, gsl_vector_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_saxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_daxpy (double alpha, const gsl_vector * X, gsl_vector * Y)
{
  if (X->size == Y->size)
    {
      cblas_daxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_caxpy (const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                gsl_vector_complex_float * Y)
{
  if (X->size == Y->size)
    {
      cblas_caxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
                   INT (X->stride), Y->data, INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zaxpy (const gsl_complex alpha, const gsl_vector_complex * X,
                gsl_vector_complex * Y)
{
  if (X->size == Y->size)
    {
      cblas_zaxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
                   INT (X->stride), Y->data, INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Generate rotation */

int
gsl_blas_srotg (float a[], float b[], float c[], float s[])
{
  cblas_srotg (a, b, c, s);
  return GSL_SUCCESS;
}

int
gsl_blas_drotg (double a[], double b[], double c[], double s[])
{
  cblas_drotg (a, b, c, s);
  return GSL_SUCCESS;
}

/* Apply rotation to vectors */

int
gsl_blas_srot (gsl_vector_float * X, gsl_vector_float * Y, float c, float s)
{
  if (X->size == Y->size)
    {
      cblas_srot (INT (X->size), X->data, INT (X->stride), Y->data,
                  INT (Y->stride), c, s);
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const double c, const double s)
{
  if (X->size == Y->size)
    {
      cblas_drot (INT (X->size), X->data, INT (X->stride), Y->data,
                  INT (Y->stride), c, s);
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Generate modified rotation */

int
gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[])
{
  cblas_srotmg (d1, d2, b1, b2, P);
  return GSL_SUCCESS;
}

int
gsl_blas_drotmg (double d1[], double d2[], double b1[], double b2, double P[])
{
  cblas_drotmg (d1, d2, b1, b2, P);
  return GSL_SUCCESS;
}

/* Apply modified rotation */

int
gsl_blas_srotm (gsl_vector_float * X, gsl_vector_float * Y, const float P[])
{
  if (X->size == Y->size)
    {
      cblas_srotm (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride), P);
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[])
{
  if (X->size == Y->size)
    {
      cblas_drotm (INT (X->size), X->data, INT (X->stride), Y->data,
                   INT (Y->stride), P);
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* Scale vector */

void
gsl_blas_sscal (float alpha, gsl_vector_float * X)
{
  cblas_sscal (INT (X->size), alpha, X->data, INT (X->stride));
}

void
gsl_blas_dscal (double alpha, gsl_vector * X)
{
  cblas_dscal (INT (X->size), alpha, X->data, INT (X->stride));
}

void
gsl_blas_cscal (const gsl_complex_float alpha, gsl_vector_complex_float * X)
{
  cblas_cscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride));
}

void
gsl_blas_zscal (const gsl_complex alpha, gsl_vector_complex * X)
{
  cblas_zscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride));
}

void
gsl_blas_csscal (float alpha, gsl_vector_complex_float * X)
{
  cblas_csscal (INT (X->size), alpha, X->data, INT (X->stride));
}

void
gsl_blas_zdscal (double alpha, gsl_vector_complex * X)
{
  cblas_zdscal (INT (X->size), alpha, X->data, INT (X->stride));
}

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

/* GEMV */

int
gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA, float alpha,
                const gsl_matrix_float * A, const gsl_vector_float * X,
                float beta, gsl_vector_float * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
      || (TransA == CblasTrans && M == X->size && N == Y->size))
    {
      cblas_sgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
                   INT (A->tda), X->data, INT (X->stride), beta, Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,
                const gsl_vector * X, double beta, gsl_vector * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
      || (TransA == CblasTrans && M == X->size && N == Y->size))
    {
      cblas_dgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
                   INT (A->tda), X->data, INT (X->stride), beta, Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_vector_complex_float * X,
                const gsl_complex_float beta, gsl_vector_complex_float * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
      || (TransA == CblasTrans && M == X->size && N == Y->size)
      || (TransA == CblasConjTrans && M == X->size && N == Y->size))
    {
      cblas_cgemv (CblasRowMajor, TransA, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
                   INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex alpha,
                const gsl_matrix_complex * A, const gsl_vector_complex * X,
                const gsl_complex beta, gsl_vector_complex * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
      || (TransA == CblasTrans && M == X->size && N == Y->size)
      || (TransA == CblasConjTrans && M == X->size && N == Y->size))
    {
      cblas_zgemv (CblasRowMajor, TransA, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
                   INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
                   INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* HEMV */

int
gsl_blas_chemv (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_vector_complex_float * X,
                const gsl_complex_float beta, gsl_vector_complex_float * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size || N != Y->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_chemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
               INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
               Y->data, INT (Y->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_zhemv (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
                const gsl_matrix_complex * A, const gsl_vector_complex * X,
                const gsl_complex beta, gsl_vector_complex * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size || N != Y->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zhemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
               INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
               Y->data, INT (Y->stride));
  return GSL_SUCCESS;
}

/* SYMV */

int
gsl_blas_ssymv (CBLAS_UPLO_t Uplo, float alpha, const gsl_matrix_float * A,
                const gsl_vector_float * X, float beta, gsl_vector_float * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size || N != Y->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
               X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_dsymv (CBLAS_UPLO_t Uplo, double alpha, const gsl_matrix * A,
                const gsl_vector * X, double beta, gsl_vector * Y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size || N != Y->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
               X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
  return GSL_SUCCESS;
}

/* TRMV */

int
gsl_blas_strmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
                gsl_vector_float * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_strmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_dtrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dtrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_ctrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
                gsl_vector_complex_float * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ctrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_ztrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
                gsl_vector_complex * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ztrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

/* TRSV */

int
gsl_blas_strsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
                gsl_vector_float * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_strsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_dtrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dtrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_ctrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
                gsl_vector_complex_float * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ctrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

int
gsl_blas_ztrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
                gsl_vector_complex * X)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N != X->size)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ztrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}

/* GER */

int
gsl_blas_sger (float alpha, const gsl_vector_float * X,
               const gsl_vector_float * Y, gsl_matrix_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_sger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
                  INT (X->stride), Y->data, INT (Y->stride), A->data,
                  INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dger (double alpha, const gsl_vector * X, const gsl_vector * Y,
               gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_dger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
                  INT (X->stride), Y->data, INT (Y->stride), A->data,
                  INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* GERU */

int
gsl_blas_cgeru (const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_cgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zgeru (const gsl_complex alpha, const gsl_vector_complex * X,
                const gsl_vector_complex * Y, gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_zgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* GERC */

int
gsl_blas_cgerc (const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_cgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zgerc (const gsl_complex alpha, const gsl_vector_complex * X,
                const gsl_vector_complex * Y, gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_zgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* HER */

int
gsl_blas_cher (CBLAS_UPLO_t Uplo, float alpha,
               const gsl_vector_complex_float * X,
               gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher (CblasRowMajor, Uplo, INT (M), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_zher (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector_complex * X,
               gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}

/* HER2 */

int
gsl_blas_cher2 (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride), Y->data, INT (Y->stride), A->data,
               INT (A->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_zher2 (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
                const gsl_vector_complex * X, const gsl_vector_complex * Y,
                gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride), Y->data, INT (Y->stride), A->data,
               INT (A->tda));
  return GSL_SUCCESS;
}

/* SYR */

int
gsl_blas_ssyr (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
               gsl_matrix_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_dsyr (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
               gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}

/* SYR2 */

int
gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
                const gsl_vector_float * Y, gsl_matrix_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
               Y->data, INT (Y->stride), A->data, INT (A->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
                const gsl_vector * Y, gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
               Y->data, INT (Y->stride), A->data, INT (A->tda));
  return GSL_SUCCESS;
}

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/* GEMM */

int
gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                float alpha, const gsl_matrix_float * A,
                const gsl_matrix_float * B, float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_sgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                double alpha, const gsl_matrix * A, const gsl_matrix * B,
                double beta, gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_dgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_cgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_zgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* SYMM */

int
gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, float alpha,
                const gsl_matrix_float * A, const gsl_matrix_float * B,
                float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_ssymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
                   A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

}

int
gsl_blas_dsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, double alpha,
                const gsl_matrix * A, const gsl_matrix * B, double beta,
                gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_dsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
                   A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_csymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_csymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_zsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* HEMM */

int
gsl_blas_chemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_chemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

}

int
gsl_blas_zhemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_zhemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* SYRK */

int
gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
                const gsl_matrix_float * A, float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_dsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
                const gsl_matrix * A, double beta, gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}

int
gsl_blas_csyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_csyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
               GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
               GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_zsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_complex beta, gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
               GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
               GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* HERK */

int
gsl_blas_cherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
                const gsl_matrix_complex_float * A, float beta,
                gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_zherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
                const gsl_matrix_complex * A, double beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != J)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
               INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* SYR2K */

int
gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
                 const gsl_matrix_float * A, const gsl_matrix_float * B,
                 float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
                INT (A->tda), B->data, INT (B->tda), beta, C->data,
                INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
                 const gsl_matrix * A, const gsl_matrix * B, double beta,
                 gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
                INT (A->tda), B->data, INT (B->tda), beta, C->data,
                INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_csyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                 const gsl_complex_float alpha,
                 const gsl_matrix_complex_float * A,
                 const gsl_matrix_complex_float * B,
                 const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_csyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
                GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

int
gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                 const gsl_complex alpha, const gsl_matrix_complex * A,
                 const gsl_matrix_complex * B, const gsl_complex beta,
                 gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
                GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* HER2K */

int
gsl_blas_cher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                 const gsl_complex_float alpha,
                 const gsl_matrix_complex_float * A,
                 const gsl_matrix_complex_float * B, float beta,
                 gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
                GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                INT (B->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}

int
gsl_blas_zher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                 const gsl_complex alpha, const gsl_matrix_complex * A,
                 const gsl_matrix_complex * B, double beta,
                 gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
                GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                INT (B->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}

/* TRMM */

int
gsl_blas_strmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
                const gsl_matrix_float * A, gsl_matrix_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_strmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dtrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
                const gsl_matrix * A, gsl_matrix * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_dtrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ctrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                gsl_matrix_complex_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ctrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ztrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                gsl_matrix_complex * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ztrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* TRSM */

int
gsl_blas_strsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
                const gsl_matrix_float * A, gsl_matrix_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_strsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_dtrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
                const gsl_matrix * A, gsl_matrix * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_dtrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ctrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                gsl_matrix_complex_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ctrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_ztrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                gsl_matrix_complex * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ztrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/* convolvedexponential.c -- model functions for multiple exponentials + background, convolved with an instrument response */
/////////////////////////////////////////////////////////////////////////////////
//
//  Calculation of multiple exponentials convolved with an instrument response.
//  using GNU Scientific Library.  Also, calls Levenberg-Marquardt minimization routine
//  for fitting event counting histograms.  Options include using the Maximum Likelihood Estimator
//  for Poisson deviates and least squares fitting.
//
//  Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the Lawrence Livermore National Laboratory. Written by Ted Laurence (laurence2@llnl.gov)
//  LLNL-CODE-424602 All rights reserved.
//  This file is part of dlevmar_mle_der
//
//  Please also read Our Notice and GNU General Public License.
//
//  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 2 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
/////////////////////////////////////////////////////////////////////////////////

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_vector.h>
#include "levmar_mle.h"

struct data {
	double deltaT;	/*Spacing between bins*/
	size_t n;		/*Number of data elements*/
	size_t nRadix2; /*Next higher power of 2 for n, giving size of irf and fWork arrays*/
	double *irfWork;/*Instrument Response Function*/
	double *fWork;  /*Workspace for function calculation*/
	size_t m_total;	/*Number of parameters (number of exponentials = (m-1)/2 ) */
	double *a;		/*Parameter Array; fitted parameters x are subsituted into this array */
	size_t m;		/*Number of fitted parameters*/
	int *ip;		/*Mapping of fitted parameter array x to parameter array a*/
};

void complex_convolution_multiply(double *a, const double *b, int n)
{
	int i,j;
	double re,im,factor;

	factor = 1.0/n; /* to normalize inverse FFT */

	a[0]*=b[0]*factor;
	a[n/2]*=b[n/2]*factor;

	for (i=1;i<n/2;i++)
	{
		j=n-i;
		re=a[i];
		im=a[j];
		a[i]=(re*b[i]-im*b[j])*factor;
		a[j]=(im*b[i]+re*b[j])*factor;
	}
}

void convolved_exp_f(double *p, double *f, int m, int n, void *data)
{
	double deltaT=((struct data *)data)->deltaT;
	size_t nRadix2 = ((struct data *)data)->nRadix2;
	double *irfWork = ((struct data *) data)->irfWork;
	double *fWork = ((struct data *) data)->fWork;
	size_t m_total = ((struct data *)data)->m_total;
	double *a = ((struct data *) data)->a;
	int *ip = ((struct data *) data)->ip;

	int i,j;
	int numExponentials = (m_total-1)/2;
	double tau, amp, sqrtTau,sqrtAmp;

	double constantValue;
	double totalT=deltaT*n;
	double factor1;
	double factor2;

	if ( (n != ((struct data *)data)->n ) || (m != ((struct data *)data)->m) ) return;

	for (i=0; i<m; i++)
		a[ ip[i] ] = p[i];

	constantValue=a[0]*a[0]/n;

	for (i=0; i<n; i++) fWork[i] =  constantValue;
	for (i=n; i<nRadix2; i++) fWork[i] = 0.0;

	for (j=0;j<numExponentials;j++)
	{
		sqrtAmp=a[2*j+1];
		sqrtTau=a[2*j+2];
		amp=sqrtAmp*sqrtAmp;
		tau=sqrtTau*sqrtTau;
		factor1=1-exp(-deltaT/tau);
		factor2=1-exp(-totalT/tau);
		for (i=0;i<n;i++)
			fWork[i]+=amp*factor1/factor2*exp(-deltaT*i/tau);
	}

	gsl_fft_real_radix2_transform(fWork,1,nRadix2);
	complex_convolution_multiply(fWork,irfWork,nRadix2);
	gsl_fft_halfcomplex_radix2_backward(fWork,1,nRadix2);

	// This loop uses 0 offset indexing
	for (i=n;i<nRadix2;i++) fWork[i%n]+=fWork[i];
	for (i=0;i<n;i++) f[i]=fWork[i];

	return;
}

void convolved_exp_df(double *p, double *jac, int m, int n, void *data)
{
	double deltaT=((struct data *)data)->deltaT;
	size_t nRadix2 = ((struct data *)data)->nRadix2;
	double *irfWork = ((struct data *) data)->irfWork;
	double *fWork = ((struct data *) data)->fWork;
	size_t m_total = ((struct data *)data)->m_total;
	double *a = ((struct data *) data)->a;
	int *ip = ((struct data *) data)->ip;
	double *jac_row;

	int i,k;

	if ( (n != ((struct data *)data)->n ) || (m != ((struct data *)data)->m) ) return;

	for (k=0; k<m; k++)
	{
		double tau, amp, sqrtTau,sqrtAmp;
		double temp, totalT=deltaT*n;
		double factor1;
		double factor2;

		if (ip[k]==0)
		{
			double d_constant=2*a[0]/n;
			jac_row=jac;
			for (i=0;i<n;i++,jac_row+=m)
				jac_row[k]=d_constant;
			continue;
		}

		for (i=0; i<nRadix2; i++) fWork[i] = 0.0;

		if ((ip[k]-1)%2==0) /* Amplitude */
		{
			sqrtAmp=a[ip[k]];
			sqrtTau=a[ip[k]+1];
			amp=sqrtAmp*sqrtAmp;
			tau=sqrtTau*sqrtTau;
			factor1=1-exp(-deltaT/tau);
			factor2=1-exp(-totalT/tau);
			for (i=0;i<n;i++)
				fWork[i]=2*sqrtAmp*factor1/factor2*exp(-deltaT*i/tau);
		}
		else /* lifetime */
		{
			sqrtAmp=a[ip[k]-1];
			sqrtTau=a[ip[k]];
			amp=sqrtAmp*sqrtAmp;
			tau=sqrtTau*sqrtTau;
			factor1=1-exp(-deltaT/tau);
			factor2=1-exp(-totalT/tau);
			for (i=0;i<n;i++)
			{
				temp=factor1/factor2*exp(-deltaT*i/tau);
				fWork[i]= -2*amp*deltaT*exp(-(deltaT*i+deltaT)/tau)/(factor2*tau*sqrtTau)
							+2*temp*amp*deltaT*i/(tau*sqrtTau)
							+2*temp*amp*totalT*exp(-totalT/tau)/(factor2*tau*sqrtTau);
			}
		}

		gsl_fft_real_radix2_transform(fWork,1,nRadix2);
		complex_convolution_multiply(fWork,irfWork,nRadix2);
		gsl_fft_halfcomplex_radix2_backward(fWork,1,nRadix2);

		// This loop uses 0 offset indexing
		for (i=n;i<nRadix2;i++) fWork[i%n]+=fWork[i];

		jac_row=jac;
		for (i=0;i<n;i++,jac_row+=m)
			jac_row[k]=fWork[i];
	}
	return;
}

int fit_convolved_exponential(int n,			// number of measurement bins
							  double deltaT,	// time resolution for each bin
							  double *x,		// Data Measured
							  double *f,		// Space for function
							  double *irf,		// Instrument Response
							int m_total,		// Total number of parameters (including those not fitted)
							double *a,			// Parameter array (length m_total)
							int m,				// Number of fitted parameters
							int *ip,			// Array of length m with indices of fitted parameters in a
							double *chisq,		// Result for chisq
							int *nIterations,	// Number of iterations performed
							int fitType,		// 0 = MLE, 1 = Neyman weighting, 2 = Equal Weighting (sigma=1)
							double deltaChisqLimit, // Stop criterion for change in chisq
							int maxIterations)	//  maximum number of iterations
{
	struct data importedData;

	double tempN,nRadix2;
	double *irfWork, *fWork;
	double opts[LM_OPTS_SZ],info[LM_INFO_SZ];

	int i;
	int status;
	unsigned iter = 0;

	double *p;

	opts[0]=LM_INIT_MU;
	opts[1]=LM_STOP_THRESH;
	opts[2]=LM_STOP_THRESH;
	opts[3]=LM_STOP_THRESH;
	opts[4]=deltaChisqLimit;

	importedData.n=n;
	importedData.deltaT=deltaT;
	importedData.m_total=m_total;
	importedData.a=a;
	importedData.m=m;
	importedData.ip=ip;

	for (i=0; i<m_total; i++) a[i]=sqrt(a[i]);

	/* Determine power of 2 higher than n required for performing convolutions via FFT */
	tempN = pow(2, 1.0+ceil( log((double)n)/log(2.0) ) );
	nRadix2 = ( (int) tempN );
	irfWork=malloc(nRadix2*sizeof(double));
	if (irfWork==NULL) return GSL_ENOMEM;
	fWork=malloc(nRadix2*sizeof(double));
	if (fWork==NULL)
	{
		free(irfWork);
		return GSL_ENOMEM;
	}
	p=malloc(m*sizeof(double));
	if (p==NULL)
	{
		free(irfWork);
		free(fWork);
		return GSL_ENOMEM;
	}

	for (i=0; i<n; i++)
		irfWork[i]=irf[i];
	for (i=n; i<nRadix2; i++)
		irfWork[i]=0.0;
	gsl_fft_real_radix2_transform(irfWork,1,nRadix2);

	importedData.nRadix2=nRadix2;
	importedData.irfWork=irfWork;
	importedData.fWork=fWork;

	for (i=0; i<m; i++)
		p[i]=a[ip[i]];

	status = dlevmar_mle_der( &convolved_exp_f, &convolved_exp_df, p, x, m, n, maxIterations,opts,info,NULL,NULL,&importedData,fitType);

	*chisq = info[1]/(n-m);
	*nIterations=info[5];

	convolved_exp_f(p, f, m, n, &importedData);  /* Evaluate function for external use in array f */

	for (i=0; i<m; i++)
		a[ip[i]]=p[i];

	free(p);
	free(fWork);
	free(irfWork);

	for (i=0; i<m_total; i++) a[i]=a[i]*a[i];

	return info[6];
}

ext: c

max_stars_repo_name: david-hoffman/dphtools

max_issues_repo_name: david-hoffman/dphutils

max_forks_repo_name: david-hoffman/dphtools


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_3_2.c
 * \brief Source file to optimize Runge-Kutta 3 steps 2nd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_3_2.h"

#define DEBUG_RK_3_2 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 3 steps 2nd order Runge-Kutta
 * method.
 */
int
rk_tb_3_2 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_tb_3_2: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t3 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  b32 (tb) = r[3];
  b31 (tb) = (0.5L - b32 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (b31 (tb)))
    return 0;
  rk_b_3 (tb);
#if DEBUG_RK_3_2
  rk_print_tb (optimize, "rk_tb_3_2", stderr);
  fprintf (stderr, "rk_tb_3_2: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 2nd order, 3rd order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_3_2t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_tb_3_2t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t3 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  b32 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b32 (tb)))
    return 0;
  b31 (tb) = (0.5L - b32 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (b31 (tb)))
    return 0;
  rk_b_3 (tb);
#if DEBUG_RK_3_2
  rk_print_tb (optimize, "rk_tb_3_2t", stderr);
  fprintf (stderr, "rk_tb_3_2t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 1st-2nd order Runge-Kutta
 * pair.
 */
int
rk_tb_3_2p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_tb_3_2p: start\n");
#endif
  if (!rk_tb_3_2 (optimize))
    return 0;
  tb = optimize->coefficient;
  e31 (tb) = 0.L;
  rk_e_3 (tb);
#if DEBUG_RK_3_2
  rk_print_e (optimize, "rk_tb_3_2p", stderr);
  fprintf (stderr, "rk_tb_3_2p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 1st-2nd order, 1st-3rd order
 * in equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_3_2tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_tb_3_2tp: start\n");
#endif
  if (!rk_tb_3_2t (optimize))
    return 0;
  tb = optimize->coefficient;
  e31 (tb) = 0.L;
  rk_e_3 (tb);
#if DEBUG_RK_3_2
  rk_print_tb (optimize, "rk_tb_3_2tp", stderr);
  fprintf (stderr, "rk_tb_3_2pt: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 3 steps 2nd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_2 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_objective_tb_3_2: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_objective_tb_3_2: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_2: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 3 steps 2nd order, 3rd
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_2t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_objective_tb_3_2t: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_2
  fprintf (stderr, "rk_objective_tb_3_2t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_2t: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
/* specfunc/ellint.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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: G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_precision.h>
#include "gsl_sf_ellint.h"

#include "error.h"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

inline
static double locMAX3(double x, double y, double z)
{
  double xy = GSL_MAX(x, y);
  return GSL_MAX(xy, z);
}

inline
static double locMAX4(double x, double y, double z, double w)
{
  double xy  = GSL_MAX(x,  y);
  double xyz = GSL_MAX(xy, z);
  return GSL_MAX(xyz, w);
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

/* based on Carlson's algorithms:
   [B. C. Carlson Numer. Math. 33, 1 (1979)]

   see also:
   [B.C. Carlson, Special Functions of Applied Mathematics (1977)]
 */

/* According to Carlson's algorithm, the errtol parameter
   typically effects the relative error in the following way:

   relative error < 16 errtol^6 / (1 - 2 errtol)

     errtol     precision
     ------     ----------
     0.001       1.0e-17
     0.003       2.0e-14
     0.01        2.0e-11
     0.03        2.0e-8
     0.1         2.0e-5
*/

int
gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result)
{
  const double lolim = 5.0 * GSL_DBL_MIN;
  const double uplim = 0.2 * GSL_DBL_MAX;
  const gsl_prec_t goal = GSL_MODE_PREC(mode);
  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  const double prec   = gsl_prec_eps[goal];

  if(x < 0.0 || y < 0.0 || x + y < lolim) {
    DOMAIN_ERROR(result);
  }
  else if(GSL_MAX(x, y) < uplim) {
    const double c1 = 1.0 / 7.0;
    const double c2 = 9.0 / 22.0;
    double xn = x;
    double yn = y;
    double mu, sn, lamda, s;
    while(1) {
      mu = (xn + yn + yn) / 3.0;
      sn = (yn + mu) / mu - 2.0;
      if (fabs(sn) < errtol) break;
      lamda = 2.0 * sqrt(xn) * sqrt(yn) + yn;
      xn = (xn + lamda) * 0.25;
      yn = (yn + lamda) * 0.25;
    }
    s = sn * sn * (0.3 + sn * (c1 + sn * (0.375 + sn * c2)));
    result->val = (1.0 + s) / sqrt(mu);
    result->err = prec * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

int
gsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
{
  const gsl_prec_t goal = GSL_MODE_PREC(mode);
  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  const double prec   = gsl_prec_eps[goal];
  const double lolim = 2.0/pow(GSL_DBL_MAX, 2.0/3.0);
  const double uplim = pow(0.1*errtol/GSL_DBL_MIN, 2.0/3.0);

  if(GSL_MIN(x,y) < 0.0 || GSL_MIN(x+y,z) < lolim) {
    DOMAIN_ERROR(result);
  }
  else if(locMAX3(x,y,z) < uplim) {
    const double c1 = 3.0 / 14.0;
    const double c2 = 1.0 /  6.0;
    const double c3 = 9.0 / 22.0;
    const double c4 = 3.0 / 26.0;
    double xn = x;
    double yn = y;
    double zn = z;
    double sigma  = 0.0;
    double power4 = 1.0;
    double ea, eb, ec, ed, ef, s1, s2;
    double mu, xndev, yndev, zndev;
    while(1) {
      double xnroot, ynroot, znroot, lamda;
      double epslon;
      mu = (xn + yn + 3.0 * zn) * 0.2;
      xndev = (mu - xn) / mu;
      yndev = (mu - yn) / mu;
      zndev = (mu - zn) / mu;
      epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
      if (epslon < errtol) break;
      xnroot = sqrt(xn);
      ynroot = sqrt(yn);
      znroot = sqrt(zn);
      lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
      sigma  += power4 / (znroot * (zn + lamda));
      power4 *= 0.25;
      xn = (xn + lamda) * 0.25;
      yn = (yn + lamda) * 0.25;
      zn = (zn + lamda) * 0.25;
    }
    ea = xndev * yndev;
    eb = zndev * zndev;
    ec = ea - eb;
    ed = ea - 6.0 * eb;
    ef = ed + ec + ec;
    s1 = ed * (- c1 + 0.25 * c3 * ed - 1.5 * c4 * zndev * ef);
    s2 = zndev * (c2 * ef + zndev * (- c3 * ec + zndev * c4 * ea));
    result->val = 3.0 * sigma + power4 * (1.0 + s1 + s2) / (mu * sqrt(mu));
    result->err = prec * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

int
gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
{
  const double lolim = 5.0 * GSL_DBL_MIN;
  const double uplim = 0.2 * GSL_DBL_MAX;
  const gsl_prec_t goal = GSL_MODE_PREC(mode);
  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  const double prec   = gsl_prec_eps[goal];

  if(x < 0.0 || y < 0.0 || z < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x+y < lolim || x+z < lolim || y+z < lolim) {
    DOMAIN_ERROR(result);
  }
  else if(locMAX3(x,y,z) < uplim) {
    const double c1 = 1.0 / 24.0;
    const double c2 = 3.0 / 44.0;
    const double c3 = 1.0 / 14.0;
    double xn = x;
    double yn = y;
    double zn = z;
    double mu, xndev, yndev, zndev, e2, e3, s;
    while(1) {
      double epslon, lamda;
      double xnroot, ynroot, znroot;
      mu = (xn + yn + zn) / 3.0;
      xndev = 2.0 - (mu + xn) / mu;
      yndev = 2.0 - (mu + yn) / mu;
      zndev = 2.0 - (mu + zn) / mu;
      epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
      if (epslon < errtol) break;
      xnroot = sqrt(xn);
      ynroot = sqrt(yn);
      znroot = sqrt(zn);
      lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
      xn = (xn + lamda) * 0.25;
      yn = (yn + lamda) * 0.25;
      zn = (zn + lamda) * 0.25;
    }
    e2 = xndev * yndev - zndev * zndev;
    e3 = xndev * yndev * zndev;
    s = 1.0 + (c1 * e2 - 0.1 - c2 * e3) * e2 + c3 * e3;
    result->val = s / sqrt(mu);
    result->err = prec * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

int
gsl_sf_ellint_RJ_e(double x, double y, double z, double p, gsl_mode_t mode, gsl_sf_result * result)
{
  const gsl_prec_t goal = GSL_MODE_PREC(mode);
  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  const double prec   = gsl_prec_eps[goal];
  const double lolim =       pow(5.0 * GSL_DBL_MIN, 1.0/3.0);
  const double uplim = 0.3 * pow(0.2 * GSL_DBL_MAX, 1.0/3.0);

  if(x < 0.0 || y < 0.0 || z < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x + y < lolim || x + z < lolim || y + z < lolim || p < lolim) {
    DOMAIN_ERROR(result);
  }
  else if(locMAX4(x,y,z,p) < uplim) {
    const double c1 = 3.0 / 14.0;
    const double c2 = 1.0 /  3.0;
    const double c3 = 3.0 / 22.0;
    const double c4 = 3.0 / 26.0;
    double xn = x;
    double yn = y;
    double zn = z;
    double pn = p;
    double sigma = 0.0;
    double power4 = 1.0;
    double mu, xndev, yndev, zndev, pndev;
    double ea, eb, ec, e2, e3, s1, s2, s3;
    while(1) {
      double xnroot, ynroot, znroot;
      double lamda, alfa, beta;
      double epslon;
      gsl_sf_result rcresult;
      int rcstatus;
      mu = (xn + yn + zn + pn + pn) * 0.2;
      xndev = (mu - xn) / mu;
      yndev = (mu - yn) / mu;
      zndev = (mu - zn) / mu;
      pndev = (mu - pn) / mu;
      epslon = locMAX4(fabs(xndev), fabs(yndev), fabs(zndev), fabs(pndev));
      if(epslon < errtol) break;
      xnroot = sqrt(xn);
      ynroot = sqrt(yn);
      znroot = sqrt(zn);
      lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
      alfa = pn * (xnroot + ynroot + znroot) + xnroot * ynroot * znroot;
      alfa = alfa * alfa;
      beta = pn * (pn + lamda) * (pn + lamda);
      rcstatus = gsl_sf_ellint_RC_e(alfa, beta, mode, &rcresult);
      if(rcstatus != GSL_SUCCESS) {
        result->val = 0.0;
	result->err = 0.0;
        return rcstatus;
      }
      sigma  += power4 * rcresult.val;
      power4 *= 0.25;
      xn = (xn + lamda) * 0.25;
      yn = (yn + lamda) * 0.25;
      zn = (zn + lamda) * 0.25;
      pn = (pn + lamda) * 0.25;
    }

    ea = xndev * (yndev + zndev) + yndev * zndev;
    eb = xndev * yndev * zndev;
    ec = pndev * pndev;
    e2 = ea - 3.0 * ec;
    e3 = eb + 2.0 * pndev * (ea - ec);
    s1 = 1.0 + e2 * (- c1 + 0.75 * c3 * e2 - 1.5 * c4 * e3);
    s2 = eb * (0.5 * c2 + pndev * (- c3 - c3 + pndev * c4));
    s3 = pndev * ea * (c2 - pndev * c3) - c2 * pndev * ec;
    result->val = 3.0 * sigma + power4 * (s1 + s2 + s3) / (mu * sqrt(mu));
    result->err = prec * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.1)] */
int
gsl_sf_ellint_F_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
{
  double sin_phi  = sin(phi);
  double sin2_phi = sin_phi*sin_phi;
  double x = 1.0 - sin2_phi;
  double y = 1.0 - k*k*sin2_phi;
  gsl_sf_result rf;
  int status = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
  result->val = sin_phi * rf.val;
  result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin_phi*rf.err);
  return status;
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.2)] */
int
gsl_sf_ellint_E_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
{
  const double sin_phi  = sin(phi);
  const double sin2_phi = sin_phi  * sin_phi;
  const double x = 1.0 - sin2_phi;
  const double y = 1.0 - k*k*sin2_phi;
  if(x < GSL_DBL_EPSILON) {
    return gsl_sf_ellint_Ecomp_e(k, mode, result);
  }
  else {
    gsl_sf_result rf;
    gsl_sf_result rd;
    const double sin3_phi = sin2_phi * sin_phi;
    const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
    const int rdstatus = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
    result->val  = sin_phi * rf.val - k*k/3.0 * sin3_phi * rd.val;
    result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
    result->err += fabs(sin_phi*rf.err);
    result->err += k*k/3.0 * GSL_DBL_EPSILON * fabs(sin3_phi * rd.val);
    result->err += k*k/3.0 * fabs(sin3_phi*rd.err);
    return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
  }
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.3)] */
int
gsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
{
  const double sin_phi  = sin(phi);
  const double sin2_phi = sin_phi  * sin_phi;
  const double sin3_phi = sin2_phi * sin_phi;
  const double x = 1.0 - sin2_phi;
  const double y = 1.0 - k*k*sin2_phi;
  gsl_sf_result rf;
  gsl_sf_result rj;
  const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
  const int rjstatus = gsl_sf_ellint_RJ_e(x, y, 1.0, 1.0 + n*sin2_phi, mode, &rj);
  result->val  = sin_phi * rf.val - n/3.0*sin3_phi * rj.val;
  result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
  result->err += n/3.0 * fabs(sin3_phi*rj.err);
  return GSL_ERROR_SELECT_2(rfstatus, rjstatus);
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.4)] */
int
gsl_sf_ellint_D_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
{
  const double sin_phi  = sin(phi);
  const double sin2_phi = sin_phi  * sin_phi;
  const double sin3_phi = sin2_phi * sin_phi;
  const double x = 1.0 - sin2_phi;
  const double y = 1.0 - k*k*sin2_phi;
  gsl_sf_result rd;
  const int status = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
  result->val = sin3_phi/3.0 * rd.val;
  result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin3_phi/3.0 * rd.err);
  return status;
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
int
gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
{
  if(k*k >= 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
    /* [Abramowitz+Stegun, 17.3.33] */
    const double y = 1.0 - k*k;
    const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
    const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
    const double ta = a[0] + y*(a[1] + y*a[2]);
    const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
    result->val = ta + tb;
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);
  }
}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.6)] */
int
gsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
{
  if(k*k >= 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
    /* [Abramowitz+Stegun, 17.3.36] */
    const double y = 1.0 - k*k;
    const double a[] = { 0.44325141463, 0.06260601220, 0.04757383546 };
    const double b[] = { 0.24998368310, 0.09200180037, 0.04069697526 };
    const double ta = 1.0 + y*(a[0] + y*(a[1] + a[2]*y));
    const double tb = -y*log(y) * (b[0] + y*(b[1] + b[2]*y));
    result->val = ta + tb;
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result rf;
    gsl_sf_result rd;
    const double y = 1.0 - k*k;
    const int rfstatus = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, &rf);
    const int rdstatus = gsl_sf_ellint_RD_e(0.0, y, 1.0, mode, &rd);
    result->val = rf.val - k*k/3.0 * rd.val;
    result->err = rf.err + k*k/3.0 * rd.err;
    return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_Kcomp_e(k, mode, &result));
}

double gsl_sf_ellint_Ecomp(double k, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_Ecomp_e(k, mode, &result));
}

double gsl_sf_ellint_F(double phi, double k, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_F_e(phi, k, mode, &result));
}

double gsl_sf_ellint_E(double phi, double k, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_E_e(phi, k, mode, &result));
}

double gsl_sf_ellint_P(double phi, double k, double n, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_P_e(phi, k, n, mode, &result));
}

double gsl_sf_ellint_D(double phi, double k, double n, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_D_e(phi, k, n, mode, &result));
}

double gsl_sf_ellint_RC(double x, double y, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_RC_e(x, y, mode, &result));
}

double gsl_sf_ellint_RD(double x, double y, double z, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_RD_e(x, y, z, mode, &result));
}

double gsl_sf_ellint_RF(double x, double y, double z, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_RF_e(x, y, z, mode, &result));
}

double gsl_sf_ellint_RJ(double x, double y, double z, double p, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_ellint_RJ_e(x, y, z, p, mode, &result));
}