text:
//Gets line spectral frequencies (LSFs).
//from the linear prediction (LP) polynomial coefficients along cols or rows of X.

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include "cmp_ascend.c"

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int poly2lsf_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim);
int poly2lsf_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim);
int poly2lsf_c (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim);
int poly2lsf_z (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim);

int poly2lsf_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const float z = 0.0f, o = 1.0f;
    const int P = (dim==0) ? R-1 : C-1;
    const int job = 'E', compz = 'N';  //eigenvalues only
    const lapack_int ldh = P+1, n = P+1, ldz = 1;
    const lapack_int ilo = 1, ihi = n;  //avoids balancing
    lapack_int info;
    float *poly, *prev, *omegas, *compan, *wr, *wi, zz[1];
    int r, c;

    //Checks
    if (R<1) { fprintf(stderr,"error in poly2lsf_s: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in poly2lsf_s: ncols X must be positive\n"); return 1; }
    if (P<1) { fprintf(stderr,"error in poly2lsf_s: P (length of polynomial coeffs including a0=1) must be positive\n"); return 1; }

    //Allocate
    if (!(poly=(float *)malloc((size_t)(n+1)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(prev=(float *)malloc((size_t)(n+1)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(omegas=(float *)malloc((size_t)(2*n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wr=(float *)malloc((size_t)(n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wi=(float *)malloc((size_t)(n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(compan=(float *)malloc((size_t)(n*n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        for (c=0; c<C; c++)
        {
            if (iscolmajor) { cblas_scopy(n,&X[c*R],1,poly,1); }
            else { cblas_scopy(n,&X[c],C,poly,1); }
            cblas_sscal(n,-1.0f/poly[0],poly,1); poly[n] = z;
            for (r=1; r<=n; r++) { prev[r] = poly[n-r]; }
            for (r=1; r<=n; r++) { poly[r] -= prev[r]; }
            cblas_scopy(n*n,&z,0,compan,1); cblas_scopy(n-1,&o,0,&compan[1],n+1); cblas_scopy(n,&poly[1],1,compan,n);
            //fprintf(stderr,"compan = \n"); for (r=0; r<n; r++) { for (c2=0; c2<n; c2++) { fprintf(stderr,"%f ",compan[r+c2*n]); } fprintf(stderr,"\n"); }
            info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_s: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[r] = (wi[r]>10.0f*FLT_EPSILON) ? atan2f(wi[r],wr[r]) : 0.0f; }
            for (r=1; r<=n; r++) { poly[r] += prev[r] + prev[r]; }
            cblas_scopy(n*n,&z,0,compan,1); cblas_scopy(n-1,&o,0,&compan[1],n+1); cblas_scopy(n,&poly[1],1,compan,n);
            //fprintf(stderr,"compan = \n"); for (r=0; r<n; r++) { for (c2=0; c2<n; c2++) { fprintf(stderr,"%f ",compan[r+c2*n]); } fprintf(stderr,"\n"); }
            info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_s: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[n+r] = (wi[r]>10.0f*FLT_EPSILON) ? atan2f(wi[r],wr[r]) : 0.0f; }
            qsort(omegas,(size_t)(2*n),sizeof(float),cmp_ascend_s);
            if (iscolmajor) { cblas_scopy(n-1,&omegas[n+1],1,&Y[c*(n-1)],1); }  //cblas_scopy(2*n,omegas,1,&Y[2*c*n],1);
            else { cblas_scopy(n-1,&omegas[n+1],1,&Y[c],C); }  //cblas_scopy(2*n,omegas,1,&Y[c],C);
        }
    }
    else if (dim==1)
    {
        for (r=0; r<R; r++)
        {
            if (iscolmajor) { cblas_scopy(n,&X[r],R,poly,1); }
            else { cblas_scopy(n,&X[r*C],1,poly,1); }
            cblas_sscal(n,-1.0f/poly[0],poly,1); poly[n] = z;
            for (c=1; c<=n; c++) { prev[c] = poly[n-c]; }
            for (c=1; c<=n; c++) { poly[c] -= prev[c]; }
            cblas_scopy(n*n,&z,0,compan,1); cblas_scopy(n-1,&o,0,&compan[1],n+1); cblas_scopy(n,&poly[1],1,compan,n);
            info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_s: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[c] = (wi[c]>10.0f*FLT_EPSILON) ? atan2f(wi[c],wr[c]) : 0.0f; }
            for (c=1; c<=n; c++) { poly[c] += prev[c] + prev[c]; }
            cblas_scopy(n*n,&z,0,compan,1); cblas_scopy(n-1,&o,0,&compan[1],n+1); cblas_scopy(n,&poly[1],1,compan,n);
            info = LAPACKE_shseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_s: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[n+c] = (wi[c]>10.0f*FLT_EPSILON) ? atan2f(wi[c],wr[c]) : 0.0f; }
            qsort(omegas,(size_t)(2*n),sizeof(float),cmp_ascend_s);
            if (iscolmajor) { cblas_scopy(n-1,&omegas[n+1],1,&Y[r],R); }  //cblas_scopy(2*n,omegas,1,&Y[r],R);
            else { cblas_scopy(n-1,&omegas[n+1],1,&Y[r*(n-1)],1); }  //cblas_scopy(2*n,omegas,1,&Y[2*r*n],1);
        }
    }
    else
    {
        fprintf(stderr,"error in poly2lsf_s: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(compan); free(wr); free(wi); free(poly); free(prev); free(omegas);
    return 0;
}

int poly2lsf_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const double z = 0.0, o = 1.0;
    const int P = (dim==0) ? R-1 : C-1;
    const int job = 'E', compz = 'N';  //eigenvalues only
    const lapack_int ldh = P+1, n = P+1, ldz = 1;
    const lapack_int ilo = 1, ihi = n;  //avoids balancing
    lapack_int info;
    double *poly, *prev, *omegas, *compan, *wr, *wi, zz[1];
    int r, c;

    //Checks
    if (R<1) { fprintf(stderr,"error in poly2lsf_d: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in poly2lsf_d: ncols X must be positive\n"); return 1; }
    if (P<1) { fprintf(stderr,"error in poly2lsf_d: P (length of polynomial coeffs including a0=1) must be positive\n"); return 1; }

    //Allocate
    if (!(poly=(double *)malloc((size_t)(n+1)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(prev=(double *)malloc((size_t)(n+1)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(omegas=(double *)malloc((size_t)(2*n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wr=(double *)malloc((size_t)(n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wi=(double *)malloc((size_t)(n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(compan=(double *)malloc((size_t)(n*n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        for (c=0; c<C; c++)
        {
            if (iscolmajor) { cblas_dcopy(n,&X[c*R],1,poly,1); }
            else { cblas_dcopy(n,&X[c],C,poly,1); }
            cblas_dscal(n,-1.0/poly[0],poly,1); poly[n] = z;
            for (r=1; r<=n; r++) { prev[r] = poly[n-r]; }
            for (r=1; r<=n; r++) { poly[r] -= prev[r]; }
            cblas_dcopy(n*n,&z,0,compan,1); cblas_dcopy(n-1,&o,0,&compan[1],n+1); cblas_dcopy(n,&poly[1],1,compan,n);
            info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_d: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[r] = (wi[r]>10.0*DBL_EPSILON) ? atan2(wi[r],wr[r]) : 0.0; }
            for (r=1; r<=n; r++) { poly[r] += prev[r] + prev[r]; }
            cblas_dcopy(n*n,&z,0,compan,1); cblas_dcopy(n-1,&o,0,&compan[1],n+1); cblas_dcopy(n,&poly[1],1,compan,n);
            info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_d: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[n+r] = (wi[r]>10.0*DBL_EPSILON) ? atan2(wi[r],wr[r]) : 0.0; }
            qsort(omegas,(size_t)(2*n),sizeof(double),cmp_ascend_d);
            if (iscolmajor) { cblas_dcopy(n-1,&omegas[n+1],1,&Y[c*(n-1)],1); }  //cblas_dcopy(2*n,omegas,1,&Y[2*c*n],1);
            else { cblas_dcopy(n-1,&omegas[n+1],1,&Y[c],C); }  //cblas_dcopy(2*n,omegas,1,&Y[c],C);
        }
    }
    else if (dim==1)
    {
        for (r=0; r<R; r++)
        {
            if (iscolmajor) { cblas_dcopy(n,&X[r],R,poly,1); }
            else { cblas_dcopy(n,&X[r*C],1,poly,1); }
            cblas_dscal(n,-1.0/poly[0],poly,1); poly[n] = z;
            for (c=1; c<=n; c++) { prev[c] = poly[n-c]; }
            for (c=1; c<=n; c++) { poly[c] -= prev[c]; }
            cblas_dcopy(n*n,&z,0,compan,1); cblas_dcopy(n-1,&o,0,&compan[1],n+1); cblas_dcopy(n,&poly[1],1,compan,n);
            info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_d: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[c] = (wi[c]>10.0*DBL_EPSILON) ? atan2(wi[c],wr[c]) : 0.0; }
            for (c=1; c<=n; c++) { poly[c] += prev[c] + prev[c]; }
            cblas_dcopy(n*n,&z,0,compan,1); cblas_dcopy(n-1,&o,0,&compan[1],n+1); cblas_dcopy(n,&poly[1],1,compan,n);
            info = LAPACKE_dhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,compan,ldh,wr,wi,zz,ldz);  //eig
            if (info) { fprintf(stderr,"error in poly2lsf_d: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[n+c] = (wi[c]>10.0*DBL_EPSILON) ? atan2(wi[c],wr[c]) : 0.0; }
            qsort(omegas,(size_t)(2*n),sizeof(double),cmp_ascend_d);
            if (iscolmajor) { cblas_dcopy(n-1,&omegas[n+1],1,&Y[r],R); }  //cblas_dcopy(2*n,omegas,1,&Y[r],R);
            else { cblas_dcopy(n-1,&omegas[n+1],1,&Y[r*(n-1)],1); }
        }
    }
    else
    {
        fprintf(stderr,"error in poly2lsf_d: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(compan); free(wr); free(wi); free(poly); free(prev); free(omegas);
    return 0;
}

int poly2lsf_c (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const float z[2] =  {0.0f,0.0f}, o[2] =  {1.0f,0.0f};
    const int P = (dim==0) ? R-1 : C-1;
    const int job = 'E', compz = 'N';  //eigenvalues only
    const lapack_int ldh = P+1, n = P+1, ldz = 1;
    const lapack_int ilo = 1, ihi = n;  //avoids balancing
    lapack_int info;
    float *poly, *prev, *omegas, *compan, *wc, zz[2], sc[2];
    int r, c;

    //Checks
    if (R<1) { fprintf(stderr,"error in poly2lsf_c: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in poly2lsf_c: ncols X must be positive\n"); return 1; }
    if (P<1) { fprintf(stderr,"error in poly2lsf_c: P (length of polynomial coeffs including a0=1) must be positive\n"); return 1; }

    //Allocate
    if (!(poly=(float *)malloc((size_t)(2*n+2)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_c: problem with malloc. "); perror("malloc"); return 1; }
    if (!(prev=(float *)malloc((size_t)(2*n+2)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_c: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wc=(float *)malloc((size_t)(2*n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_c: problem with malloc. "); perror("malloc"); return 1; }
    if (!(omegas=(float *)malloc((size_t)(2*n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_c: problem with malloc. "); perror("malloc"); return 1; }
    if (!(compan=(float *)malloc((size_t)(4*n*n)*sizeof(float)))) { fprintf(stderr,"error in poly2lsf_c: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        for (c=0; c<C; c++)
        {
            if (iscolmajor) { cblas_ccopy(n,&X[2*c*R],1,poly,1); }
            else { cblas_ccopy(n,&X[2*c],C,poly,1); }
            sc[0] = 1.0f/(poly[0]*poly[0]+poly[1]*poly[1]); sc[1] = poly[1]*sc[0]; sc[0] *= -poly[0];
            cblas_cscal(n,sc,poly,1); poly[2*n] = poly[2*n+1] = 0.0f;
            for (r=1; r<=n; r++) { prev[2*r] = poly[2*(n-r)]; prev[2*r+1] = poly[2*(n-r)+1]; }
            for (r=1; r<=n; r++) { poly[2*r] -= prev[2*r]; poly[2*r+1] -= prev[2*r+1]; }
            cblas_ccopy(n*n,z,0,compan,1); cblas_ccopy(n-1,o,0,&compan[2],n+1); cblas_ccopy(n,&poly[2],1,compan,n);
            info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)wc,(lapack_complex_float *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_c: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[r] = atan2f(wc[2*r+1],wc[2*r]); }
            //for (r=0; r<n; r++) { omegas[r] = (wc[2*r+1]>10.0f*FLT_EPSILON) ? atan2f(wc[2*r+1],wc[2*r]) : 0.0f; }
            for (r=1; r<=n; r++) { poly[2*r] += prev[2*r] + prev[2*r]; poly[2*r+1] += prev[2*r+1] + prev[2*r+1]; }
            cblas_ccopy(n*n,z,0,compan,1); cblas_ccopy(n-1,o,0,&compan[2],n+1); cblas_ccopy(n,&poly[2],1,compan,n);
            info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)wc,(lapack_complex_float *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_c: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[n+r] = atan2f(wc[2*r+1],wc[2*r]); }
            //for (r=0; r<n; r++) { omegas[n+r] = (wc[2*r+1]>10.0f*FLT_EPSILON) ? atan2f(wc[2*r+1],wc[2*r]) : 0.0f; }
            qsort(omegas,(size_t)(2*n),sizeof(float),cmp_ascend_s);
            if (iscolmajor) { cblas_scopy(2*n,omegas,1,&Y[2*c*n],1); }
            else { cblas_scopy(2*n,omegas,1,&Y[c],C); }
        }
    }
    else if (dim==1)
    {
        for (r=0; r<R; r++)
        {
            if (iscolmajor) { cblas_ccopy(n,&X[2*r],R,poly,1); }
            else { cblas_ccopy(n,&X[2*r*C],1,poly,1); }
            sc[0] = 1.0f/(poly[0]*poly[0]+poly[1]*poly[1]); sc[1] = poly[1]*sc[0]; sc[0] *= -poly[0];
            cblas_cscal(n,sc,poly,1); poly[2*n] = poly[2*n+1] = 0.0f;
            for (c=1; c<=n; c++) { prev[2*c] = poly[2*(n-c)]; prev[2*c+1] = poly[2*(n-c)+1]; }
            for (c=1; c<=n; c++) { poly[2*c] -= prev[2*c]; poly[2*c+1] -= prev[2*c+1]; }
            cblas_ccopy(n*n,z,0,compan,1); cblas_ccopy(n-1,o,0,&compan[2],n+1); cblas_ccopy(n,&poly[2],1,compan,n);
            info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)wc,(lapack_complex_float *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_c: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[c] = atan2f(wc[2*c+1],wc[2*c]); }
            //for (c=0; c<n; c++) { omegas[c] = (wc[2*r+1]>10.0f*FLT_EPSILON) ? atan2f(wc[2*r+1],wc[2*r]) : 0.0f; }
            for (c=1; c<=n; c++) { poly[2*c] += prev[2*c] + prev[2*c]; poly[2*c+1] += prev[2*c+1] + prev[2*c+1]; }
            cblas_ccopy(n*n,z,0,compan,1); cblas_ccopy(n-1,o,0,&compan[2],n+1); cblas_ccopy(n,&poly[2],1,compan,n);
            info = LAPACKE_chseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_float *)compan,ldh,(lapack_complex_float *)wc,(lapack_complex_float *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_c: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[n+c] = atan2f(wc[2*c+1],wc[2*c]); }
            //for (c=0; c<n; c++) { omegas[n+c] = (wc[2*r+1]>10.0f*FLT_EPSILON) ? atan2f(wc[2*r+1],wc[2*r]) : 0.0f; }
            qsort(omegas,(size_t)(2*n),sizeof(float),cmp_ascend_s);
            if (iscolmajor) { cblas_scopy(2*n,omegas,1,&Y[r],R); }
            else { cblas_scopy(2*n,omegas,1,&Y[2*r*n],1); }
        }
    }
    else
    {
        fprintf(stderr,"error in poly2lsf_c: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(compan); free(wc); free(poly); free(prev); free(omegas);
    return 0;
}

int poly2lsf_z (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim)
{
    const double z[2] =  {0.0,0.0}, o[2] =  {1.0,0.0};
    const int P = (dim==0) ? R-1 : C-1;
    const int job = 'E', compz = 'N';  //eigenvalues only
    const lapack_int ldh = P+1, n = P+1, ldz = 1;
    const lapack_int ilo = 1, ihi = n;  //avoids balancing
    lapack_int info;
    double *poly, *prev, *omegas, *compan, *wc, zz[2], sc[2];
    int r, c;

    //Checks
    if (R<1) { fprintf(stderr,"error in poly2lsf_z: nrows X must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in poly2lsf_z: ncols X must be positive\n"); return 1; }
    if (P<1) { fprintf(stderr,"error in poly2lsf_z: P (length of polynomial coeffs including a0=1) must be positive\n"); return 1; }

    //Allocate
    if (!(poly=(double *)malloc((size_t)(2*n+2)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_z: problem with malloc. "); perror("malloc"); return 1; }
    if (!(prev=(double *)malloc((size_t)(2*n+2)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_z: problem with malloc. "); perror("malloc"); return 1; }
    if (!(wc=(double *)malloc((size_t)(2*n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_z: problem with malloc. "); perror("malloc"); return 1; }
    if (!(omegas=(double *)malloc((size_t)(2*n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_z: problem with malloc. "); perror("malloc"); return 1; }
    if (!(compan=(double *)malloc((size_t)(4*n*n)*sizeof(double)))) { fprintf(stderr,"error in poly2lsf_z: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0)
    {
        for (c=0; c<C; c++)
        {
            if (iscolmajor) { cblas_zcopy(n,&X[2*c*R],1,poly,1); }
            else { cblas_zcopy(n,&X[2*c],C,poly,1); }
            sc[0] = 1.0/(poly[0]*poly[0]+poly[1]*poly[1]); sc[1] = poly[1]*sc[0]; sc[0] *= -poly[0];
            cblas_zscal(n,sc,poly,1); poly[2*n] = poly[2*n+1] = 0.0;
            for (r=1; r<=n; r++) { prev[2*r] = poly[2*(n-r)]; prev[2*r+1] = poly[2*(n-r)+1]; }
            for (r=1; r<=n; r++) { poly[2*r] -= prev[2*r]; poly[2*r+1] -= prev[2*r+1]; }
            cblas_zcopy(n*n,z,0,compan,1); cblas_zcopy(n-1,o,0,&compan[2],n+1); cblas_zcopy(n,&poly[2],1,compan,n);
            info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)wc,(lapack_complex_double *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_z: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[r] = atan2(wc[2*r+1],wc[2*r]); }
            //for (r=0; r<n; r++) { omegas[r] = (wc[2*r+1]>10.0*DBL_EPSILON) ? atan2(wc[2*r+1],wc[2*r]) : 0.0; }
            for (r=1; r<=n; r++) { poly[2*r] += prev[2*r] + prev[2*r]; poly[2*r+1] += prev[2*r+1] + prev[2*r+1]; }
            cblas_zcopy(n*n,z,0,compan,1); cblas_zcopy(n-1,o,0,&compan[2],n+1); cblas_zcopy(n,&poly[2],1,compan,n);
            info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)wc,(lapack_complex_double *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_z: lapacke decomposition failed\n"); return 1; }
            for (r=0; r<n; r++) { omegas[n+r] = atan2(wc[2*r+1],wc[2*r]); }
            qsort(omegas,(size_t)(2*n),sizeof(double),cmp_ascend_d);
            if (iscolmajor) { cblas_dcopy(2*n,omegas,1,&Y[2*c*n],1); }
            else { cblas_dcopy(2*n,omegas,1,&Y[c],C); }
        }
    }
    else if (dim==1)
    {
        for (r=0; r<R; r++)
        {
            if (iscolmajor) { cblas_zcopy(n,&X[2*r],R,poly,1); }
            else { cblas_zcopy(n,&X[2*r*C],1,poly,1); }
            sc[0] = 1.0/(poly[0]*poly[0]+poly[1]*poly[1]); sc[1] = poly[1]*sc[0]; sc[0] *= -poly[0];
            cblas_zscal(n,sc,poly,1); poly[2*n] = poly[2*n+1] = 0.0;
            for (c=1; c<=n; c++) { prev[2*c] = poly[2*(n-c)]; prev[2*c+1] = poly[2*(n-c)+1]; }
            for (c=1; c<=n; c++) { poly[2*c] -= prev[2*c]; poly[2*c+1] -= prev[2*c+1]; }
            cblas_zcopy(n*n,z,0,compan,1); cblas_zcopy(n-1,o,0,&compan[2],n+1); cblas_zcopy(n,&poly[2],1,compan,n);
            info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)wc,(lapack_complex_double *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_z: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[c] = atan2(wc[2*c+1],wc[2*c]); }
            //for (c=0; c<n; c++) { omegas[n+c] = (wc[2*r+1]>10.0*DBL_EPSILON) ? atan2(wc[2*r+1],wc[2*r]) : 0.0; }
            for (c=1; c<=n; c++) { poly[2*c] += prev[2*c] + prev[2*c]; poly[2*c+1] += prev[2*c+1] + prev[2*c+1]; }
            cblas_zcopy(n*n,z,0,compan,1); cblas_zcopy(n-1,o,0,&compan[2],n+1); cblas_zcopy(n,&poly[2],1,compan,n);
            info = LAPACKE_zhseqr(LAPACK_COL_MAJOR,job,compz,n,ilo,ihi,(lapack_complex_double *)compan,ldh,(lapack_complex_double *)wc,(lapack_complex_double *)zz,ldz);
            if (info) { fprintf(stderr,"error in poly2lsf_z: lapacke decomposition failed\n"); return 1; }
            for (c=0; c<n; c++) { omegas[n+c] = atan2(wc[2*c+1],wc[2*c]); }
            //for (c=0; c<n; c++) { omegas[n+c] = (wc[2*r+1]>10.0*DBL_EPSILON) ? atan2(wc[2*r+1],wc[2*r]) : 0.0; }
            qsort(omegas,(size_t)(2*n),sizeof(double),cmp_ascend_d);
            if (iscolmajor) { cblas_dcopy(2*n,omegas,1,&Y[r],R); }
            else { cblas_dcopy(2*n,omegas,1,&Y[2*r*n],1); }
        }
    }
    else
    {
        fprintf(stderr,"error in poly2lsf_z: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    free(compan); free(wc); free(poly); free(prev); free(omegas);
    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include "player/instrument.h"
#include "synth/synth.h"
#include <fmidi/fmidi.h>
#include <gsl/gsl>
#include <memory>

class Midi_Synth_Instrument : public Midi_Instrument {
public:
    Midi_Synth_Instrument();
    ~Midi_Synth_Instrument();

    void flush_events() override;

    void open_midi_output(gsl::cstring_span id) override;
    void close_midi_output() override;

    bool is_synth() const override { return true; }

    void configure_audio(double audio_rate, double audio_latency);
    void generate_audio(float *output, unsigned nframes);

    void preload(const fmidi_smf_t &smf);

protected:
    void handle_send_message(const uint8_t *data, unsigned len, double ts, uint8_t flags) override;

private:
    struct Impl;
    std::unique_ptr<Impl> impl_;
};

ext: h

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
static const char help[] =
"Solves elasto-plastic torsion problem in 2D using SNESVI.  Option prefix -el_.\n"
"Equation is  - grad^2 u = f(x,y)  where  f(x,y) = 2 C.  Zero Dirichlet\n"
"boundary conditions.  Domain Omega = (0,1)x(0,1).  Constraint is upper bound\n"
"   u(x,y) <= dist((x,y), bdry Omega) = psi(x,y).\n"
"At locations where the constraint is active, the material experiences plastic\n"
"failure.  Where inactive, the equation represents the elasticity.  As with\n"
"related codes ../obstacle.c and dam.c the code reuses the residual and\n"
"Jacobian evaluation code from ch6/.\n"
"Reference: R. Kornhuber (1994) 'Monotone multigrid methods for elliptic\n"
"variational inequalities I', Numerische Mathematik, 69(2), 167-184.\n\n";

#include <petsc.h>
#include "../../ch6/poissonfunctions.h"

// z = psi(x,y) = dist((x,y), bdry Omega)  is the upper obstacle
double psi(double x, double y) {
    if (y <= x)
        if (y <= 1.0 - x)
            return y;
        else
            return 1.0 - x;
    else
        if (y <= 1.0 - x)
            return x;
        else
            return 1.0 - y;
}

double zero(double x, double y, double z, void *ctx) {
    return 0.0;
}

static double C = 2.5;
double f_fcn(double x, double y, double z, void *ctx) {
    return 2.0 * C;
}

extern PetscErrorCode FormBounds(SNES, Vec, Vec);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  DM             da, da_after;
  SNES           snes;
  Vec            u_initial, u;
  PoissonCtx     user;
  SNESConvergedReason reason;
  int            snesits;
  double         lflops,flops;
  DMDALocalInfo  info;

  PetscInitialize(&argc,&argv,NULL,help);

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"el_",
                           "elasto-plastic torsion solver options",""); CHKERRQ(ierr);
  ierr = PetscOptionsReal("-C","f(x,y)=2C is source term",
                          "elasto.c",C,&C,NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
      DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
      3,3,                       // override with -da_grid_x,_y
      PETSC_DECIDE,PETSC_DECIDE, // num of procs in each dim
      1,1,NULL,NULL,             // dof = 1 and stencil width = 1
      &da);CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,-1.0,-1.0);CHKERRQ(ierr);

  user.cx = 1.0;
  user.cy = 1.0;
  user.cz = 1.0;
  user.g_bdry = &zero;
  user.f_rhs = &f_fcn;
  user.addctx = NULL;
  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);

  ierr = SNESSetType(snes,SNESVINEWTONRSLS);CHKERRQ(ierr);
  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds);CHKERRQ(ierr);

  // reuse residual and jacobian from ch6/
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)Poisson2DFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  // initial iterate is zero
  ierr = DMCreateGlobalVector(da,&u_initial);CHKERRQ(ierr);
  ierr = VecSet(u_initial,0.0); CHKERRQ(ierr);

  /* solve; then get solution and DM after solution*/
  ierr = SNESSolve(snes,NULL,u_initial);CHKERRQ(ierr);
  ierr = VecDestroy(&u_initial); CHKERRQ(ierr);
  ierr = DMDestroy(&da); CHKERRQ(ierr);
  ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
  ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr); /* do not destroy u */

  /* performance measures */
  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
  if (reason <= 0) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "WARNING: SNES not converged ... use -snes_converged_reason to check\n"); CHKERRQ(ierr);
  }
  ierr = SNESGetIterationNumber(snes,&snesits); CHKERRQ(ierr);
  ierr = PetscGetFlops(&lflops); CHKERRQ(ierr);
  ierr = MPI_Allreduce(&lflops,&flops,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "done on %4d x %4d grid; total flops = %.3e; SNES iterations %d\n",
      info.mx,info.my,flops,snesits); CHKERRQ(ierr);

  SNESDestroy(&snes);
  return PetscFinalize();
}

// for call-back: tell SNESVI we want  -infty < u <= psi
PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu) {
  PetscErrorCode ierr;
  DM            da;
  DMDALocalInfo info;
  int           i, j;
  double        **aXu, dx, dy, x, y;
  ierr = VecSet(Xu,PETSC_NINFINITY);CHKERRQ(ierr);
  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  dx = 1.0 / (PetscReal)(info.mx-1);
  dy = 1.0 / (PetscReal)(info.my-1);
  ierr = DMDAVecGetArray(da, Xu, &aXu);CHKERRQ(ierr);
  for (j=info.ys; j<info.ys+info.ym; j++) {
    y = j * dy;
    for (i=info.xs; i<info.xs+info.xm; i++) {
      x = i * dx;
      aXu[j][i] = psi(x,y);
    }
  }
  ierr = DMDAVecRestoreArray(da, Xu, &aXu);CHKERRQ(ierr);
  return 0;
}


text:
#include <gsl/gsl_complex.h>

#define RVEC(A) int A##n, double*A##p
#define RMAT(A) int A##r, int A##c, double* A##p
#define KRVEC(A) int A##n, const double*A##p
#define KRMAT(A) int A##r, int A##c, const double* A##p

#define CVEC(A) int A##n, gsl_complex*A##p
#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
#define KCVEC(A) int A##n, const gsl_complex*A##p
#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p

#define FVEC(A) int A##n, float*A##p
#define FMAT(A) int A##r, int A##c, float* A##p
#define KFVEC(A) int A##n, const float*A##p
#define KFMAT(A) int A##r, int A##c, const float* A##p

#define QVEC(A) int A##n, gsl_complex_float*A##p
#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
#define KQVEC(A) int A##n, const gsl_complex_float*A##p
#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p

#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_multifit_nlin.h>
#include <string.h>
#include <stdio.h>

#define MACRO(B) do {B} while (0)
#define ERROR(CODE) MACRO(return CODE;)
#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
#define OK return 0;

#define MIN(A,B) ((A)<(B)?(A):(B))
#define MAX(A,B) ((A)>(B)?(A):(B))

#ifdef DBG
#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
#else
#define DEBUGMSG(M)
#endif

#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)

#ifdef DBG
#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
#else
#define DEBUGMAT(MSG,X)
#endif

#ifdef DBG
#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
#else
#define DEBUGVEC(MSG,X)
#endif

#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)

#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c)

#define V(a) (&a.vector)
#define M(a) (&a.matrix)

#define GCVEC(A) int A##n, gsl_complex*A##p
#define KGCVEC(A) int A##n, const gsl_complex*A##p

#define GQVEC(A) int A##n, gsl_complex_float*A##p
#define KGQVEC(A) int A##n, const gsl_complex_float*A##p

#define BAD_SIZE 2000
#define BAD_CODE 2001
#define MEM      2002
#define BAD_FILE 2003

void no_abort_on_error() {
    gsl_set_error_handler_off();
}

int sumF(KFVEC(x),FVEC(r)) {
    DEBUGMSG("sumF");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    float res = 0;
    for (i = 0; i < xn; i++) res += xp[i];
    rp[0] = res;
    OK
}

int sumR(KRVEC(x),RVEC(r)) {
    DEBUGMSG("sumR");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    double res = 0;
    for (i = 0; i < xn; i++) res += xp[i];
    rp[0] = res;
    OK
}

int sumQ(KQVEC(x),QVEC(r)) {
    DEBUGMSG("sumQ");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    gsl_complex_float res;
    res.dat[0] = 0;
    res.dat[1] = 0;
    for (i = 0; i < xn; i++) {
      res.dat[0] += xp[i].dat[0];
      res.dat[1] += xp[i].dat[1];
    }
    rp[0] = res;
    OK
}

int sumC(KCVEC(x),CVEC(r)) {
    DEBUGMSG("sumC");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    gsl_complex res;
    res.dat[0] = 0;
    res.dat[1] = 0;
    for (i = 0; i < xn; i++)  {
      res.dat[0] += xp[i].dat[0];
      res.dat[1] += xp[i].dat[1];
    }
    rp[0] = res;
    OK
}

int prodF(KFVEC(x),FVEC(r)) {
    DEBUGMSG("prodF");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    float res = 1;
    for (i = 0; i < xn; i++) res *= xp[i];
    rp[0] = res;
    OK
}

int prodR(KRVEC(x),RVEC(r)) {
    DEBUGMSG("prodR");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    double res = 1;
    for (i = 0; i < xn; i++) res *= xp[i];
    rp[0] = res;
    OK
}

int prodQ(KQVEC(x),QVEC(r)) {
    DEBUGMSG("prodQ");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    gsl_complex_float res;
    float temp;
    res.dat[0] = 1;
    res.dat[1] = 0;
    for (i = 0; i < xn; i++) {
      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
      res.dat[0] = temp;
    }
    rp[0] = res;
    OK
}

int prodC(KCVEC(x),CVEC(r)) {
    DEBUGMSG("prodC");
    REQUIRES(rn==1,BAD_SIZE);
    int i;
    gsl_complex res;
    double temp;
    res.dat[0] = 1;
    res.dat[1] = 0;
    for (i = 0; i < xn; i++)  {
      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
      res.dat[0] = temp;
    }
    rp[0] = res;
    OK
}

int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
    DEBUGMSG("dotF");
    REQUIRES(xn==yn,BAD_SIZE);
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("dotF");
    KFVVIEW(x);
    KFVVIEW(y);
    gsl_blas_sdot(V(x),V(y),rp);
    OK
}

int dotR(KRVEC(x), KRVEC(y), RVEC(r)) {
    DEBUGMSG("dotR");
    REQUIRES(xn==yn,BAD_SIZE);
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("dotR");
    KDVVIEW(x);
    KDVVIEW(y);
    gsl_blas_ddot(V(x),V(y),rp);
    OK
}

int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) {
    DEBUGMSG("dotQ");
    REQUIRES(xn==yn,BAD_SIZE);
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("dotQ");
    KQVVIEW(x);
    KQVVIEW(y);
    gsl_blas_cdotu(V(x),V(y),rp);
    OK
}

int dotC(KCVEC(x), KCVEC(y), CVEC(r)) {
    DEBUGMSG("dotC");
    REQUIRES(xn==yn,BAD_SIZE);
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("dotC");
    KCVVIEW(x);
    KCVVIEW(y);
    gsl_blas_zdotu(V(x),V(y),rp);
    OK
}

int toScalarR(int code, KRVEC(x), RVEC(r)) {
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("toScalarR");
    KDVVIEW(x);
    double res;
    switch(code) {
        case 0: { res = gsl_blas_dnrm2(V(x)); break; }
        case 1: { res = gsl_blas_dasum(V(x));  break; }
        case 2: { res = gsl_vector_max_index(V(x));  break; }
        case 3: { res = gsl_vector_max(V(x));  break; }
        case 4: { res = gsl_vector_min_index(V(x)); break; }
        case 5: { res = gsl_vector_min(V(x)); break; }
        default: ERROR(BAD_CODE);
    }
    rp[0] = res;
    OK
}

int toScalarF(int code, KFVEC(x), FVEC(r)) {
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("toScalarF");
    KFVVIEW(x);
    float res;
    switch(code) {
        case 0: { res = gsl_blas_snrm2(V(x)); break; }
        case 1: { res = gsl_blas_sasum(V(x));  break; }
        case 2: { res = gsl_vector_float_max_index(V(x));  break; }
        case 3: { res = gsl_vector_float_max(V(x));  break; }
        case 4: { res = gsl_vector_float_min_index(V(x)); break; }
        case 5: { res = gsl_vector_float_min(V(x)); break; }
        default: ERROR(BAD_CODE);
    }
    rp[0] = res;
    OK
}

int toScalarC(int code, KCVEC(x), RVEC(r)) {
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("toScalarC");
    KCVVIEW(x);
    double res;
    switch(code) {
        case 0: { res = gsl_blas_dznrm2(V(x)); break; }
        case 1: { res = gsl_blas_dzasum(V(x));  break; }
        default: ERROR(BAD_CODE);
    }
    rp[0] = res;
    OK
}

int toScalarQ(int code, KQVEC(x), FVEC(r)) {
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("toScalarQ");
    KQVVIEW(x);
    float res;
    switch(code) {
        case 0: { res = gsl_blas_scnrm2(V(x)); break; }
        case 1: { res = gsl_blas_scasum(V(x));  break; }
        default: ERROR(BAD_CODE);
    }
    rp[0] = res;
    OK
}

inline double sign(double x) {
    if(x>0) {
        return +1.0;
    } else if (x<0) {
        return -1.0;
    } else {
        return 0.0;
    }
}

inline float float_sign(float x) {
    if(x>0) {
        return +1.0;
    } else if (x<0) {
        return -1.0;
    } else {
        return 0.0;
    }
}

inline gsl_complex complex_abs(gsl_complex z) {
    gsl_complex r;
    r.dat[0] = gsl_complex_abs(z);
    r.dat[1] = 0;
    return r;
}

inline gsl_complex complex_signum(gsl_complex z) {
    gsl_complex r;
    double mag;
    if (z.dat[0] == 0 && z.dat[1] == 0) {
        r.dat[0] = 0;
        r.dat[1] = 0;
    } else {
        mag = gsl_complex_abs(z);
        r.dat[0] = z.dat[0]/mag;
        r.dat[1] = z.dat[1]/mag;
    }
    return r;
}

#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
int mapR(int code, KRVEC(x), RVEC(r)) {
    int k;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapR");
    switch (code) {
        OP(0,sin)
        OP(1,cos)
        OP(2,tan)
        OP(3,fabs)
        OP(4,asin)
        OP(5,acos)
        OP(6,atan) /* atan2 mediante vectorZip */
        OP(7,sinh)
        OP(8,cosh)
        OP(9,tanh)
        OP(10,gsl_asinh)
        OP(11,gsl_acosh)
        OP(12,gsl_atanh)
        OP(13,exp)
        OP(14,log)
        OP(15,sign)
        OP(16,sqrt)
        default: ERROR(BAD_CODE);
    }
}

int mapF(int code, KFVEC(x), FVEC(r)) {
    int k;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapF");
    switch (code) {
        OP(0,sin)
        OP(1,cos)
        OP(2,tan)
        OP(3,fabs)
        OP(4,asin)
        OP(5,acos)
        OP(6,atan) /* atan2 mediante vectorZip */
        OP(7,sinh)
        OP(8,cosh)
        OP(9,tanh)
        OP(10,gsl_asinh)
        OP(11,gsl_acosh)
        OP(12,gsl_atanh)
        OP(13,exp)
        OP(14,log)
        OP(15,sign)
        OP(16,sqrt)
        default: ERROR(BAD_CODE);
    }
}

int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
    int k;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapC");
    switch (code) {
        OP(0,gsl_complex_sin)
        OP(1,gsl_complex_cos)
        OP(2,gsl_complex_tan)
        OP(3,complex_abs)
        OP(4,gsl_complex_arcsin)
        OP(5,gsl_complex_arccos)
        OP(6,gsl_complex_arctan)
        OP(7,gsl_complex_sinh)
        OP(8,gsl_complex_cosh)
        OP(9,gsl_complex_tanh)
        OP(10,gsl_complex_arcsinh)
        OP(11,gsl_complex_arccosh)
        OP(12,gsl_complex_arctanh)
        OP(13,gsl_complex_exp)
        OP(14,gsl_complex_log)
        OP(15,complex_signum)
        OP(16,gsl_complex_sqrt)

        // gsl_complex_arg
        // gsl_complex_abs
        default: ERROR(BAD_CODE);
    }
}

int mapC(int code, KCVEC(x), CVEC(r)) {
    return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
}

gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
{
  gsl_complex c;
  gsl_complex r;

  gsl_complex_float float_r;

  c.dat[0] = a.dat[0];
  c.dat[1] = a.dat[1];

  r = (*cf)(c);

  float_r.dat[0] = r.dat[0];
  float_r.dat[1] = r.dat[1];

  return float_r;
}

gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex),
					gsl_complex_float a,gsl_complex_float b)
{
  gsl_complex c1;
  gsl_complex c2;
  gsl_complex r;

  gsl_complex_float float_r;

  c1.dat[0] = a.dat[0];
  c1.dat[1] = a.dat[1];

  c2.dat[0] = b.dat[0];
  c2.dat[1] = b.dat[1];

  r = (*cf)(c1,c2);

  float_r.dat[0] = r.dat[0];
  float_r.dat[1] = r.dat[1];

  return float_r;
}

#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
    int k;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapQ");
    switch (code) {
        OPC(0,gsl_complex_sin)
        OPC(1,gsl_complex_cos)
        OPC(2,gsl_complex_tan)
        OPC(3,complex_abs)
        OPC(4,gsl_complex_arcsin)
        OPC(5,gsl_complex_arccos)
        OPC(6,gsl_complex_arctan)
        OPC(7,gsl_complex_sinh)
        OPC(8,gsl_complex_cosh)
        OPC(9,gsl_complex_tanh)
        OPC(10,gsl_complex_arcsinh)
        OPC(11,gsl_complex_arccosh)
        OPC(12,gsl_complex_arctanh)
        OPC(13,gsl_complex_exp)
        OPC(14,gsl_complex_log)
        OPC(15,complex_signum)
        OPC(16,gsl_complex_sqrt)

        // gsl_complex_arg
        // gsl_complex_abs
        default: ERROR(BAD_CODE);
    }
}

int mapQ(int code, KQVEC(x), QVEC(r)) {
    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
}

int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
    int k;
    double val = *pval;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapValR");
    switch (code) {
        OPV(0,val*xp[k])
        OPV(1,val/xp[k])
        OPV(2,val+xp[k])
        OPV(3,val-xp[k])
        OPV(4,pow(val,xp[k]))
        OPV(5,pow(xp[k],val))
        default: ERROR(BAD_CODE);
    }
}

int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
    int k;
    float val = *pval;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapValF");
    switch (code) {
        OPV(0,val*xp[k])
        OPV(1,val/xp[k])
        OPV(2,val+xp[k])
        OPV(3,val-xp[k])
        OPV(4,pow(val,xp[k]))
        OPV(5,pow(xp[k],val))
        default: ERROR(BAD_CODE);
    }
}

int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
    int k;
    gsl_complex val = *pval;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapValC");
    switch (code) {
        OPV(0,gsl_complex_mul(val,xp[k]))
        OPV(1,gsl_complex_div(val,xp[k]))
        OPV(2,gsl_complex_add(val,xp[k]))
        OPV(3,gsl_complex_sub(val,xp[k]))
        OPV(4,gsl_complex_pow(val,xp[k]))
        OPV(5,gsl_complex_pow(xp[k],val))
        default: ERROR(BAD_CODE);
    }
}

int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
    return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
}

int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
    int k;
    gsl_complex_float val = *pval;
    REQUIRES(xn == rn,BAD_SIZE);
    DEBUGMSG("mapValQ");
    switch (code) {
        OPCA(0,gsl_complex_mul,val,xp[k])
	OPCA(1,gsl_complex_div,val,xp[k])
	OPCA(2,gsl_complex_add,val,xp[k])
	OPCA(3,gsl_complex_sub,val,xp[k])
	OPCA(4,gsl_complex_pow,val,xp[k])
	OPCA(5,gsl_complex_pow,xp[k],val)
        default: ERROR(BAD_CODE);
    }
}

int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
}

#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
    REQUIRES(an == bn && an == rn, BAD_SIZE);
    int k;
    switch(code) {
        OPZE(4,"zipR Pow",pow)
        OPZE(5,"zipR ATan2",atan2)
    }
    KDVVIEW(a);
    KDVVIEW(b);
    DVVIEW(r);
    gsl_vector_memcpy(V(r),V(a));
    int res;
    switch(code) {
        OPZV(0,"zipR Add",gsl_vector_add)
        OPZV(1,"zipR Sub",gsl_vector_sub)
        OPZV(2,"zipR Mul",gsl_vector_mul)
        OPZV(3,"zipR Div",gsl_vector_div)
        default: ERROR(BAD_CODE);
    }
}

int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
    REQUIRES(an == bn && an == rn, BAD_SIZE);
    int k;
    switch(code) {
        OPZE(4,"zipF Pow",pow)
        OPZE(5,"zipF ATan2",atan2)
    }
    KFVVIEW(a);
    KFVVIEW(b);
    FVVIEW(r);
    gsl_vector_float_memcpy(V(r),V(a));
    int res;
    switch(code) {
        OPZV(0,"zipF Add",gsl_vector_float_add)
        OPZV(1,"zipF Sub",gsl_vector_float_sub)
        OPZV(2,"zipF Mul",gsl_vector_float_mul)
        OPZV(3,"zipF Div",gsl_vector_float_div)
        default: ERROR(BAD_CODE);
    }
}

int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
    REQUIRES(an == bn && an == rn, BAD_SIZE);
    int k;
    switch(code) {
        OPZE(0,"zipC Add",gsl_complex_add)
        OPZE(1,"zipC Sub",gsl_complex_sub)
        OPZE(2,"zipC Mul",gsl_complex_mul)
        OPZE(3,"zipC Div",gsl_complex_div)
        OPZE(4,"zipC Pow",gsl_complex_pow)
        //OPZE(5,"zipR ATan2",atan2)
    }
    //KCVVIEW(a);
    //KCVVIEW(b);
    //CVVIEW(r);
    //gsl_vector_memcpy(V(r),V(a));
    //int res;
    switch(code) {
        default: ERROR(BAD_CODE);
    }
}

int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
    return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
}

#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
    REQUIRES(an == bn && an == rn, BAD_SIZE);
    int k;
    switch(code) {
        OPCZE(0,"zipQ Add",gsl_complex_add)
        OPCZE(1,"zipQ Sub",gsl_complex_sub)
        OPCZE(2,"zipQ Mul",gsl_complex_mul)
        OPCZE(3,"zipQ Div",gsl_complex_div)
        OPCZE(4,"zipQ Pow",gsl_complex_pow)
        //OPZE(5,"zipR ATan2",atan2)
    }
    //KCVVIEW(a);
    //KCVVIEW(b);
    //CVVIEW(r);
    //gsl_vector_memcpy(V(r),V(a));
    //int res;
    switch(code) {
        default: ERROR(BAD_CODE);
    }
}

int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
}

int fft(int code, KCVEC(X), CVEC(R)) {
    REQUIRES(Xn == Rn,BAD_SIZE);
    DEBUGMSG("fft");
    int s = Xn;
    gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
    gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
    gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
    gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
    gsl_blas_dcopy(&X.vector,&R.vector);
    if(code==0) {
        gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
    } else {
        gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
    }
    gsl_fft_complex_wavetable_free (wavetable);
    gsl_fft_complex_workspace_free (workspace);
    OK
}

int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
{
    gsl_function F;
    F.function = f;
    F.params = 0;

    if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);

    if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);

    if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);

    return 0;
}

int integrate_qng(double f(double, void*), double a, double b, double prec,
                   double *result, double*error) {
    DEBUGMSG("integrate_qng");
    gsl_function F;
    F.function = f;
    F.params = NULL;
    size_t neval;
    int res = gsl_integration_qng (&F, a,b, 0, prec, result, error, &neval);
    CHECK(res,res);
    OK
}

int integrate_qags(double f(double,void*), double a, double b, double prec, int w,
               double *result, double* error) {
    DEBUGMSG("integrate_qags");
    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
    gsl_function F;
    F.function = f;
    F.params = NULL;
    int res = gsl_integration_qags (&F, a,b, 0, prec, w,wk, result, error);
    CHECK(res,res);
    gsl_integration_workspace_free (wk);
    OK
}

int integrate_qagi(double f(double,void*), double prec, int w,
               double *result, double* error) {
    DEBUGMSG("integrate_qagi");
    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
    gsl_function F;
    F.function = f;
    F.params = NULL;
    int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error);
    CHECK(res,res);
    gsl_integration_workspace_free (wk);
    OK
}

int integrate_qagiu(double f(double,void*), double a, double prec, int w,
               double *result, double* error) {
    DEBUGMSG("integrate_qagiu");
    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
    gsl_function F;
    F.function = f;
    F.params = NULL;
    int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error);
    CHECK(res,res);
    gsl_integration_workspace_free (wk);
    OK
}

int integrate_qagil(double f(double,void*), double b, double prec, int w,
               double *result, double* error) {
    DEBUGMSG("integrate_qagil");
    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
    gsl_function F;
    F.function = f;
    F.params = NULL;
    int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error);
    CHECK(res,res);
    gsl_integration_workspace_free (wk);
    OK
}

int polySolve(KRVEC(a), CVEC(z)) {
    DEBUGMSG("polySolve");
    REQUIRES(an>1,BAD_SIZE);
    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
    int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
    CHECK(res,res);
    gsl_poly_complex_workspace_free (w);
    OK;
}

int vector_fscanf(char*filename, RVEC(a)) {
    DEBUGMSG("gsl_vector_fscanf");
    DVVIEW(a);
    FILE * f = fopen(filename,"r");
    CHECK(!f,BAD_FILE);
    int res = gsl_vector_fscanf(f,V(a));
    CHECK(res,res);
    fclose (f);
    OK
}

int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
    DEBUGMSG("gsl_vector_fprintf");
    DVVIEW(a);
    FILE * f = fopen(filename,"w");
    CHECK(!f,BAD_FILE);
    int res = gsl_vector_fprintf(f,V(a),fmt);
    CHECK(res,res);
    fclose (f);
    OK
}

int vector_fread(char*filename, RVEC(a)) {
    DEBUGMSG("gsl_vector_fread");
    DVVIEW(a);
    FILE * f = fopen(filename,"r");
    CHECK(!f,BAD_FILE);
    int res = gsl_vector_fread(f,V(a));
    CHECK(res,res);
    fclose (f);
    OK
}

int vector_fwrite(char*filename, RVEC(a)) {
    DEBUGMSG("gsl_vector_fwrite");
    DVVIEW(a);
    FILE * f = fopen(filename,"w");
    CHECK(!f,BAD_FILE);
    int res = gsl_vector_fwrite(f,V(a));
    CHECK(res,res);
    fclose (f);
    OK
}

int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
    DEBUGMSG("matrix_fprintf");
    FILE * f = fopen(filename,"w");
    CHECK(!f,BAD_FILE);
    int i,j,sr,sc;
    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
    #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
    for (i=0; i<mr; i++) {
        for (j=0; j<mc-1; j++) {
            fprintf(f,fmt,AT(m,i,j));
            fprintf(f," ");
        }
        fprintf(f,fmt,AT(m,i,j));
        fprintf(f,"\n");
    }
    fclose (f);
    OK
}

//---------------------------------------------------------------

typedef double Trawfun(int, double*);

double only_f_aux_min(const gsl_vector*x, void *pars) {
    Trawfun * f = (Trawfun*) pars;
    double* p = (double*)calloc(x->size,sizeof(double));
    int k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }
    double res = f(x->size,p);
    free(p);
    return res;
}

// this version returns info about intermediate steps
int minimize(int method, double f(int, double*), double tolsize, int maxit,
                 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
    REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
    DEBUGMSG("minimizeList (nmsimplex)");
    gsl_multimin_function my_func;
    // extract function from pars
    my_func.f = only_f_aux_min;
    my_func.n = xin;
    my_func.params = f;
    size_t iter = 0;
    int status;
    double size;
    const gsl_multimin_fminimizer_type *T;
    gsl_multimin_fminimizer *s = NULL;
    // Initial vertex size vector
    KDVVIEW(sz);
    // Starting point
    KDVVIEW(xi);
    // Minimizer nmsimplex, without derivatives
    switch(method) {
        case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
#ifdef GSL110
        case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
#else
        case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; }
#endif
        default: ERROR(BAD_CODE);
    }
    s = gsl_multimin_fminimizer_alloc (T, my_func.n);
    gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
    do {
        status = gsl_multimin_fminimizer_iterate (s);
        size = gsl_multimin_fminimizer_size (s);

        solp[iter*solc+0] = iter+1;
        solp[iter*solc+1] = s->fval;
        solp[iter*solc+2] = size;

        int k;
        for(k=0;k<xin;k++) {
            solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
        }
        iter++;
        if (status) break;
        status = gsl_multimin_test_size (size, tolsize);
    } while (status == GSL_CONTINUE && iter < maxit);
    int i,j;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        for(j=1;j<solc;j++) {
            solp[i*solc+j]=0.;
        }
    }
    gsl_multimin_fminimizer_free(s);
    OK
}

// working with the gradient

typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;

double f_aux_min(const gsl_vector*x, void *pars) {
    Tfdf * fdf = ((Tfdf*) pars);
    double* p = (double*)calloc(x->size,sizeof(double));
    int k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }
    double res = fdf->f(x->size,p);
    free(p);
    return res;
}

void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
    Tfdf * fdf = ((Tfdf*) pars);
    double* p = (double*)calloc(x->size,sizeof(double));
    double* q = (double*)calloc(g->size,sizeof(double));
    int k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }

    fdf->df(x->size,p,g->size,q);

    for(k=0;k<x->size;k++) {
        gsl_vector_set(g,k,q[k]);
    }
    free(p);
    free(q);
}

void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
    *f = f_aux_min(x,pars);
    df_aux_min(x,pars,g);
}

int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
              double initstep, double minimpar, double tolgrad, int maxit,
              KRVEC(xi), RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
    DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
    gsl_multimin_function_fdf my_func;
    // extract function from pars
    my_func.f = f_aux_min;
    my_func.df = df_aux_min;
    my_func.fdf = fdf_aux_min;
    my_func.n = xin;
    Tfdf stfdf;
    stfdf.f = f;
    stfdf.df = df;
    my_func.params = &stfdf;
    size_t iter = 0;
    int status;
    const gsl_multimin_fdfminimizer_type *T;
    gsl_multimin_fdfminimizer *s = NULL;
    // Starting point
    KDVVIEW(xi);
    // conjugate gradient fr
    switch(method) {
        case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
        case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; }
        case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; }
        case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
        case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; }
        default: ERROR(BAD_CODE);
    }
    s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
    gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
    do {
        status = gsl_multimin_fdfminimizer_iterate (s);
        solp[iter*solc+0] = iter+1;
        solp[iter*solc+1] = s->f;
        int k;
        for(k=0;k<xin;k++) {
            solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
        }
        iter++;
        if (status) break;
        status = gsl_multimin_test_gradient (s->gradient, tolgrad);
    } while (status == GSL_CONTINUE && iter < maxit);
    int i,j;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        for(j=1;j<solc;j++) {
            solp[i*solc+j]=0.;
        }
    }
    gsl_multimin_fdfminimizer_free(s);
    OK
}

//---------------------------------------------------------------

double only_f_aux_root(double x, void *pars) {
    double (*f)(double) = (double (*)(double)) pars;
    return f(x);
}

int root(int method, double f(double),
         double epsrel, int maxit,
         double xl, double xu, RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
    DEBUGMSG("root_only_f");
    gsl_function my_func;
    // extract function from pars
    my_func.function = only_f_aux_root;
    my_func.params = f;
    size_t iter = 0;
    int status;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;
    // Starting point
    switch(method) {
        case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
        case 1 : {T = gsl_root_fsolver_falsepos; break; }
        case 2 : {T = gsl_root_fsolver_brent; break; }
        default: ERROR(BAD_CODE);
    }
    s = gsl_root_fsolver_alloc (T);
    gsl_root_fsolver_set (s, &my_func, xl, xu);
    do {
           double best, current_lo, current_hi;
           status = gsl_root_fsolver_iterate (s);
           best = gsl_root_fsolver_root (s);
           current_lo = gsl_root_fsolver_x_lower (s);
           current_hi = gsl_root_fsolver_x_upper (s);
           solp[iter*solc] = iter + 1;
           solp[iter*solc+1] = best;
           solp[iter*solc+2] = current_lo;
           solp[iter*solc+3] = current_hi;
           iter++;
           if (status)   /* check if solver is stuck */
             break;

           status =
               gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
        }
        while (status == GSL_CONTINUE && iter < maxit);
    int i;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        solp[i*solc+1]=0.;
        solp[i*solc+2]=0.;
        solp[i*solc+3]=0.;
    }
    gsl_root_fsolver_free(s);
    OK
}

typedef struct {
    double (*f)(double);
    double (*jf)(double);
} uniTfjf;

double f_aux_uni(double x, void *pars) {
    uniTfjf * fjf = ((uniTfjf*) pars);
    return (fjf->f)(x);
}

double jf_aux_uni(double x, void * pars) {
    uniTfjf * fjf = ((uniTfjf*) pars);
    return (fjf->jf)(x);
}

void fjf_aux_uni(double x, void * pars, double * f, double * g) {
    *f = f_aux_uni(x,pars);
    *g = jf_aux_uni(x,pars);
}

int rootj(int method, double f(double),
          double df(double),
         double epsrel, int maxit,
         double x, RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
    DEBUGMSG("root_fjf");
    gsl_function_fdf my_func;
    // extract function from pars
    my_func.f = f_aux_uni;
    my_func.df = jf_aux_uni;
    my_func.fdf = fjf_aux_uni;
    uniTfjf stfjf;
    stfjf.f = f;
    stfjf.jf = df;
    my_func.params = &stfjf;
    size_t iter = 0;
    int status;
    const gsl_root_fdfsolver_type *T;
    gsl_root_fdfsolver *s;
    // Starting point
    switch(method) {
        case 0 : {T = gsl_root_fdfsolver_newton;; break; }
        case 1 : {T = gsl_root_fdfsolver_secant; break; }
        case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
        default: ERROR(BAD_CODE);
    }
    s = gsl_root_fdfsolver_alloc (T);

    gsl_root_fdfsolver_set (s, &my_func, x);

    do {
           double x0;
           status = gsl_root_fdfsolver_iterate (s);
           x0 = x;
           x = gsl_root_fdfsolver_root(s);
           solp[iter*solc+0] = iter+1;
           solp[iter*solc+1] = x;

           iter++;
           if (status)   /* check if solver is stuck */
             break;

           status =
               gsl_root_test_delta (x, x0, 0, epsrel);
        }
        while (status == GSL_CONTINUE && iter < maxit);

    int i;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        solp[i*solc+1]=0.;
    }
    gsl_root_fdfsolver_free(s);
    OK
}

//---------------------------------------------------------------

typedef void TrawfunV(int, double*, int, double*);

int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
    TrawfunV * f = (TrawfunV*) pars;
    double* p = (double*)calloc(x->size,sizeof(double));
    double* q = (double*)calloc(y->size,sizeof(double));
    int k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }
    f(x->size,p,y->size,q);
    for(k=0;k<y->size;k++) {
        gsl_vector_set(y,k,q[k]);
    }
    free(p);
    free(q);
    return 0; //hmmm
}

int multiroot(int method, void f(int, double*, int, double*),
         double epsabs, int maxit,
         KRVEC(xi), RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
    DEBUGMSG("root_only_f");
    gsl_multiroot_function my_func;
    // extract function from pars
    my_func.f = only_f_aux_multiroot;
    my_func.n = xin;
    my_func.params = f;
    size_t iter = 0;
    int status;
    const gsl_multiroot_fsolver_type *T;
    gsl_multiroot_fsolver *s;
    // Starting point
    KDVVIEW(xi);
    switch(method) {
        case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
        case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
        case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
        case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
        default: ERROR(BAD_CODE);
    }
    s = gsl_multiroot_fsolver_alloc (T, my_func.n);
    gsl_multiroot_fsolver_set (s, &my_func, V(xi));

    do {
           status = gsl_multiroot_fsolver_iterate (s);

           solp[iter*solc+0] = iter+1;

           int k;
           for(k=0;k<xin;k++) {
               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
           }
           for(k=xin;k<2*xin;k++) {
               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
           }

           iter++;
           if (status)   /* check if solver is stuck */
             break;

           status =
             gsl_multiroot_test_residual (s->f, epsabs);
        }
        while (status == GSL_CONTINUE && iter < maxit);

    int i,j;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        for(j=1;j<solc;j++) {
            solp[i*solc+j]=0.;
        }
    }
    gsl_multiroot_fsolver_free(s);
    OK
}

// working with the jacobian

typedef struct {int (*f)(int, double*, int, double *);
                int (*jf)(int, double*, int, int, double*);} Tfjf;

int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
    Tfjf * fjf = ((Tfjf*) pars);
    double* p = (double*)calloc(x->size,sizeof(double));
    double* q = (double*)calloc(y->size,sizeof(double));
    int k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }
    (fjf->f)(x->size,p,y->size,q);
    for(k=0;k<y->size;k++) {
        gsl_vector_set(y,k,q[k]);
    }
    free(p);
    free(q);
    return 0;
}

int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
    Tfjf * fjf = ((Tfjf*) pars);
    double* p = (double*)calloc(x->size,sizeof(double));
    double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
    int i,j,k;
    for(k=0;k<x->size;k++) {
        p[k] = gsl_vector_get(x,k);
    }

    (fjf->jf)(x->size,p,jac->size1,jac->size2,q);

    k=0;
    for(i=0;i<jac->size1;i++) {
        for(j=0;j<jac->size2;j++){
            gsl_matrix_set(jac,i,j,q[k++]);
        }
    }
    free(p);
    free(q);
    return 0;
}

int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
    f_aux(x,pars,f);
    jf_aux(x,pars,g);
    return 0;
}

int multirootj(int method, int f(int, double*, int, double*),
                      int jac(int, double*, int, int, double*),
         double epsabs, int maxit,
         KRVEC(xi), RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
    DEBUGMSG("root_fjf");
    gsl_multiroot_function_fdf my_func;
    // extract function from pars
    my_func.f = f_aux;
    my_func.df = jf_aux;
    my_func.fdf = fjf_aux;
    my_func.n = xin;
    Tfjf stfjf;
    stfjf.f = f;
    stfjf.jf = jac;
    my_func.params = &stfjf;
    size_t iter = 0;
    int status;
    const gsl_multiroot_fdfsolver_type *T;
    gsl_multiroot_fdfsolver *s;
    // Starting point
    KDVVIEW(xi);
    switch(method) {
        case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; }
        case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; }
        case 2 : {T = gsl_multiroot_fdfsolver_newton; break; }
        case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; }
        default: ERROR(BAD_CODE);
    }
    s = gsl_multiroot_fdfsolver_alloc (T, my_func.n);

    gsl_multiroot_fdfsolver_set (s, &my_func, V(xi));

    do {
           status = gsl_multiroot_fdfsolver_iterate (s);

           solp[iter*solc+0] = iter+1;

           int k;
           for(k=0;k<xin;k++) {
               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
           }
           for(k=xin;k<2*xin;k++) {
               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
           }

           iter++;
           if (status)   /* check if solver is stuck */
             break;

           status =
             gsl_multiroot_test_residual (s->f, epsabs);
        }
        while (status == GSL_CONTINUE && iter < maxit);

    int i,j;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        for(j=1;j<solc;j++) {
            solp[i*solc+j]=0.;
        }
    }
    gsl_multiroot_fdfsolver_free(s);
    OK
}

//-------------- non linear least squares fitting -------------------

int nlfit(int method, int f(int, double*, int, double*),
                      int jac(int, double*, int, int, double*),
         double epsabs, double epsrel, int maxit, int p,
         KRVEC(xi), RMAT(sol)) {
    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
    DEBUGMSG("nlfit");
    const gsl_multifit_fdfsolver_type *T;
    gsl_multifit_fdfsolver *s;
    gsl_multifit_function_fdf my_f;
    // extract function from pars
    my_f.f = f_aux;
    my_f.df = jf_aux;
    my_f.fdf = fjf_aux;
    my_f.n = p;
    my_f.p = xin;  // !!!!
    Tfjf stfjf;
    stfjf.f = f;
    stfjf.jf = jac;
    my_f.params = &stfjf;
    size_t iter = 0;
    int status;

    KDVVIEW(xi);
    //DMVIEW(cov);

    switch(method) {
        case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
        case 1 : { T = gsl_multifit_fdfsolver_lmder; break; }
        default: ERROR(BAD_CODE);
    }

    s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
    gsl_multifit_fdfsolver_set (s, &my_f, V(xi));

    do {   status = gsl_multifit_fdfsolver_iterate (s);

           solp[iter*solc+0] = iter+1;
           solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);

           int k;
           for(k=0;k<xin;k++) {
               solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
           }

           iter++;
           if (status)   /* check if solver is stuck */
             break;

           status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
        }
        while (status == GSL_CONTINUE && iter < maxit);

    int i,j;
    for (i=iter; i<solr; i++) {
        solp[i*solc+0] = iter;
        for(j=1;j<solc;j++) {
            solp[i*solc+j]=0.;
        }
    }

    //gsl_multifit_covar (s->J, 0.0, M(cov));

    gsl_multifit_fdfsolver_free (s);
    OK
}

//////////////////////////////////////////////////////

#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }

int random_vector(int seed, int code, RVEC(r)) {
    DEBUGMSG("random_vector")
    static gsl_rng * gen = NULL;
    if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
    gsl_rng_set (gen, seed);
    int k;
    switch (code) {
        RAN(0,gsl_rng_uniform)
        RAN(1,gsl_ran_ugaussian)
        default: ERROR(BAD_CODE);
    }
}
#undef RAN

//////////////////////////////////////////////////////

#include "gsl-ode.c"

//////////////////////////////////////////////////////

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
#include "scf.h"
#include "scf_decompose.h"
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_gegenbauer.h>

/* computes a!/b! =  a(a-1)(a-2)...(b+1) */
double factorial_frac(int a, int b)
{
    if(a<0 || b<0)
    {
        perror("factorial_frac: either a or b is negative!!\n");
        exit(1);
    }

    int recip = 0;
    if(a==b)
        return 1.0;
    else if(a<b)
    {
        int temp = b;
        b = a;
        a = temp;
        recip = 1;
    }

    double val = 1;

    while(a > b)
    {
        val = val * a;
        a--;
    }

    if(recip)
        val = 1.0/val;

    return val;

}

void scf_integral(double* mat_cos, double* mat_sin)
{
    long long k;
    double* rr = (double*)malloc(NumPart*sizeof(double));
    double* cos_tt = (double*)malloc(NumPart*sizeof(double));
    double* pp = (double*)malloc(NumPart*sizeof(double));

#ifdef USE_GSL_GEGENBAUER_HOST
    double* xxi = (double*)malloc(NumPart*sizeof(double));
    double this_Cna;
#endif

    /* cos_tt is calculated directly by cos(theta)=z/r because
    we never need theta by itself, but only cos(theta) in the
    legendre polynomials. Also, notice that cos_tt must be
    calculated in double precision because with single precision
    zz[k] can be greater than rr[k] by a tiny bit, which
    the plgndr routine will complain later on */

    for(k=0; k<NumPart; k++)
    {
        rr[k] = sqrt((double)xx[k]*(double)xx[k] +
                      (double)yy[k]*(double)yy[k] +
                      (double)zz[k]*(double)zz[k]);

        cos_tt[k] = (double)zz[k]/rr[k];

        pp[k] = atan2(yy[k], xx[k]);

#ifdef USE_GSL_GEGENBAUER_HOST
        xxi[k] = (rr[k]-1.0) / (rr[k]+1.0);
#endif
    }

    double N_lm;
    double K_nl, I_nl, A_nl;
    double gamm1, gamm2;

    float this_Phi_nl;
    double* P_lm_k = (double*)malloc(NumPart*sizeof(double));

    double m_Phinlk_Plm, Nlm_Anl;
    double sum_cos = 0.0;
    double sum_sin = 0.0;

    int n, l, m;
    for(l=0; l<=decompose_lmax; l++)
    {
        for(m=0; m<=l; m++)
        {
            N_lm = (2*l+1)/4.0/PI * 2.0 * factorial_frac(l-m,l+m);
            if(m==0)
                N_lm = 0.5*N_lm;

            for (k=0; k<NumPart; k++)
            {
                P_lm_k[k] = gsl_sf_legendre_Plm(l,m, cos_tt[k]);

                /* Note: This can be computed more efficiently using gsl_sf_legendre_Plm_array
                which computes a list of legendre polynomials at constant m with varying l.
                but this changes the code and I don't have time to deal with that for now.
                Numerical Recipes's plgndr is actually faster, but single precision. */
            }

            for(n=0; n<=decompose_nmax; n++)
            {
                K_nl = 0.5*n*(n+4*l+3)+(l+1)*(2*l+1);
                gamm1 = factorial_frac(n+4*l+2,n);
                gamm2 = factorial_frac(2*l+1, 2*(2*l+1));
                I_nl = -K_nl/(n+2*l+1.5) * gamm1 * gamm2 * gamm2;
                A_nl = 1.0/I_nl;

                sum_cos = 0.0;
                sum_sin = 0.0;

                for(k=0; k<NumPart; k++)
                {
#ifdef USE_GSL_GEGENBAUER_HOST
                    this_Cna = gsl_sf_gegenpoly_n(n, 2*l+1.5, xxi[k]);
                    this_Phi_nl = Phi_nl(this_Cna, l, rr[k]);
#else
                    this_Phi_nl = Phi_nl(n, l, rr[k]);
#endif

                    m_Phinlk_Plm = mm[k] * this_Phi_nl * P_lm_k[k];

                    sum_cos += m_Phinlk_Plm * cos(m*pp[k]);
                    sum_sin += m_Phinlk_Plm * sin(m*pp[k]);
                }

                Nlm_Anl = N_lm * A_nl;

                mat_cos[ind(n,l,m)] = sum_cos * Nlm_Anl;
                mat_sin[ind(n,l,m)] = sum_sin * Nlm_Anl;

            } // n loop

        } // m loop

    } // l loop

    free(rr);
    free(pp);
    free(cos_tt);
    free(P_lm_k);

}

void scf_integral_spherical(float* mat_cos)
{
    long long k;
    float* rr = (float*)malloc(NumPart*sizeof(float));
    double temp_r;

    for(k=0; k<NumPart; k++)
    {
        temp_r = sqrt(xx[k]*xx[k] + yy[k]*yy[k] + zz[k]*zz[k]);
        rr[k] = (float)temp_r;
    }

    double N00 = 0.07957747154; /* = 1/(4*pi) */
    double K_n0, I_n0, A_n0;

    float this_Phi_n0;

    double sum_cos = 0.0;

    int n;

    for(n=0; n<=decompose_nmax; n++)
    {
        K_n0 = 0.5*n*(n+3)+1;
        I_n0 = -K_n0/(n+1.5)/4.0 * (n+2) * (n+1);
        A_n0 = 1.0/I_n0;

        sum_cos = 0.0;

        for(k=0; k<NumPart; k++)
        {
            this_Phi_n0 = Phi_n0(n, rr[k]);
            sum_cos += mm[k] * this_Phi_n0;
        }

        mat_cos[n] = (float) (sum_cos * A_n0 * N00);

    } // n loop

    free(rr);

}


text:
/* multifit_nlinear/dogleg.c
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * 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.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

/*
 * This module contains an implementation of the Powell dogleg
 * algorithm for nonlinear optimization problems. This implementation
 * closely follows the following works:
 *
 * [1] H. B. Nielsen, K. Madsen, Introduction to Optimization and
 *     Data Fitting, Informatics and Mathematical Modeling,
 *     Technical University of Denmark (DTU), 2010.
 *
 * [2] J. E. Dennis and H. H. W. Mei, Two new unconstrained optimization
 *     algorithms which use function and gradient values, J. Opt. Theory and
 *     Appl., 28(4), 1979.
 */

typedef struct
{
  size_t n;                  /* number of observations */
  size_t p;                  /* number of parameters */
  gsl_vector *dx_gn;         /* Gauss-Newton step, size p */
  gsl_vector *dx_sd;         /* steepest descent step, size p */
  double norm_Dgn;           /* || D dx_gn || */
  double norm_Dsd;           /* || D dx_sd || */
  double norm_Dinvg;         /* || D^{-1} g || */
  double norm_JDinv2g;       /* || J D^{-2} g || */
  gsl_vector *workp;         /* workspace, length p */
  gsl_vector *workn;         /* workspace, length n */

  /* tunable parameters */
  gsl_multifit_nlinear_parameters params;
} dogleg_state_t;

#include "common.c"

static void * dogleg_alloc (const void * params, const size_t n, const size_t p);
static void dogleg_free(void *vstate);
static int dogleg_init(const void *vtrust_state, void *vstate);
static int dogleg_preloop(const void * vtrust_state, void * vstate);
static int dogleg_step(const void * vtrust_state, const double delta,
                       gsl_vector * dx, void * vstate);
static int dogleg_double_step(const void * vtrust_state, const double delta,
                              gsl_vector * dx, void * vstate);
static int dogleg_preduction(const void * vtrust_state, const gsl_vector * dx,
                             double * pred, void * vstate);
static int dogleg_calc_gn(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx);
static double dogleg_beta(const double t, const double delta,
                          const gsl_vector * diag, dogleg_state_t * state);

static void *
dogleg_alloc (const void * params, const size_t n, const size_t p)
{
  const gsl_multifit_nlinear_parameters *mparams = (const gsl_multifit_nlinear_parameters *) params;
  dogleg_state_t *state;

  state = calloc(1, sizeof(dogleg_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate dogleg state", GSL_ENOMEM);
    }

  state->dx_gn = gsl_vector_alloc(p);
  if (state->dx_gn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for dx_gn", GSL_ENOMEM);
    }

  state->dx_sd = gsl_vector_alloc(p);
  if (state->dx_sd == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for dx_sd", GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workp", GSL_ENOMEM);
    }

  state->workn = gsl_vector_alloc(n);
  if (state->workn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workn", GSL_ENOMEM);
    }

  state->n = n;
  state->p = p;
  state->params = *mparams;

  return state;
}

static void
dogleg_free(void *vstate)
{
  dogleg_state_t *state = (dogleg_state_t *) vstate;

  if (state->dx_gn)
    gsl_vector_free(state->dx_gn);

  if (state->dx_sd)
    gsl_vector_free(state->dx_sd);

  if (state->workp)
    gsl_vector_free(state->workp);

  if (state->workn)
    gsl_vector_free(state->workn);

  free(state);
}

/*
dogleg_init()
  Initialize dogleg solver

Inputs: vtrust_state - trust state
        vstate       - workspace

Return: success/error
*/

static int
dogleg_init(const void *vtrust_state, void *vstate)
{
  (void)vtrust_state;
  (void)vstate;

  return GSL_SUCCESS;
}

/*
dogleg_preloop()
  Initialize dogleg method prior to iteration loop.
This involves computing the steepest descent step. The
Gauss-Newton step is computed if needed in the _step()
routines.

Notes: on output,
1) state->dx_sd contains steepest descent step
2) state->norm_Dinvg contains || D^{-1} g ||
3) state->norm_JDinv2g contains || J D^{-2} g ||
*/

static int
dogleg_preloop(const void * vtrust_state, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  dogleg_state_t *state = (dogleg_state_t *) vstate;
  double u;
  double alpha; /* ||g||^2 / ||Jg||^2 */

  /* calculate the steepest descent step */

  /* compute workp = D^{-1} g and its norm */
  gsl_vector_memcpy(state->workp, trust_state->g);
  gsl_vector_div(state->workp, trust_state->diag);
  state->norm_Dinvg = gsl_blas_dnrm2(state->workp);

  /* compute workp = D^{-2} g */
  gsl_vector_div(state->workp, trust_state->diag);

  /* compute: workn = J D^{-2} g */
  gsl_blas_dgemv(CblasNoTrans, 1.0, trust_state->J, state->workp, 0.0, state->workn);
  state->norm_JDinv2g = gsl_blas_dnrm2(state->workn);

  u = state->norm_Dinvg / state->norm_JDinv2g;
  alpha = u * u;

  /* dx_sd = -alpha D^{-2} g */
  gsl_vector_memcpy(state->dx_sd, state->workp);
  gsl_vector_scale(state->dx_sd, -alpha);

  state->norm_Dsd = scaled_enorm(trust_state->diag, state->dx_sd);
  state->norm_Dgn = -1.0; /* computed later if needed */

  return GSL_SUCCESS;
}

/*
dogleg_step()
  Calculate a new step vector
*/

static int
dogleg_step(const void * vtrust_state, const double delta,
            gsl_vector * dx, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  dogleg_state_t *state = (dogleg_state_t *) vstate;

  if (state->norm_Dsd >= delta)
    {
      /* steepest descent step is outside trust region;
       * truncate steepest descent step to trust region boundary */
      gsl_vector_memcpy(dx, state->dx_sd);
      gsl_vector_scale(dx, delta / state->norm_Dsd);
    }
  else
    {
      /* compute Gauss-Newton step if needed */
      if (state->norm_Dgn < 0.0)
        {
          int status = dogleg_calc_gn(trust_state, state->dx_gn);

          if (status)
            return status;

          /* compute || D dx_gn || */
          state->norm_Dgn = scaled_enorm(trust_state->diag, state->dx_gn);
        }

      if (state->norm_Dgn <= delta)
        {
          /* Gauss-Newton step is inside trust region, use it as final step
           * since it is the global minimizer of the quadratic model function */
          gsl_vector_memcpy(dx, state->dx_gn);
        }
      else
        {
          /* Gauss-Newton step is outside trust region, but steepest
           * descent is inside; use dogleg step */

          double beta = dogleg_beta(1.0, delta, trust_state->diag, state);

          /* compute: workp = dx_gn - dx_sd */
          scaled_addition(1.0, state->dx_gn, -1.0, state->dx_sd, state->workp);

          /* dx = dx_sd + beta*(dx_gn - dx_sd) */
          scaled_addition(beta, state->workp, 1.0, state->dx_sd, dx);
        }
    }

  return GSL_SUCCESS;
}

/*
dogleg_double_step()
  Calculate a new step with double dogleg method. Based on
section 3 of [2]
*/

static int
dogleg_double_step(const void * vtrust_state, const double delta,
                   gsl_vector * dx, void * vstate)
{
  const double alpha_fac = 0.8; /* recommended value from Dennis and Mei */
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  dogleg_state_t *state = (dogleg_state_t *) vstate;

  if (state->norm_Dsd >= delta)
    {
      /* steepest descent step is outside trust region;
       * truncate steepest descent step to trust region boundary */
      gsl_vector_memcpy(dx, state->dx_sd);
      gsl_vector_scale(dx, delta / state->norm_Dsd);
    }
  else
    {
      /* compute Gauss-Newton step if needed */
      if (state->norm_Dgn < 0.0)
        {
          int status = dogleg_calc_gn(trust_state, state->dx_gn);

          if (status)
            return status;

          /* compute || D dx_gn || */
          state->norm_Dgn = scaled_enorm(trust_state->diag, state->dx_gn);
        }

      if (state->norm_Dgn <= delta)
        {
          /* Gauss-Newton step is inside trust region, use it as final step
           * since it is the global minimizer of the quadratic model function */
          gsl_vector_memcpy(dx, state->dx_gn);
        }
      else
        {
          double t, u, v, c;

          /* compute: u = ||D^{-1} g||^2 / ||J D^{-2} g||^2 */
          v = state->norm_Dinvg / state->norm_JDinv2g;
          u = v * v;

          /* compute: v = g^T dx_gn */
          gsl_blas_ddot(trust_state->g, state->dx_gn, &v);

          /* compute: c = ||D^{-1} g||^4 / (||J D^{-2} g||^2 * |g^T dx_gn|) */
          c = u * (state->norm_Dinvg / fabs(v)) * state->norm_Dinvg;

          /* compute: t = 1 - alpha_fac*(1-c) */
          t = 1.0 - alpha_fac*(1.0 - c);

          if (t * state->norm_Dgn <= delta)
            {
              /* set dx = (delta / ||D dx_gn||) dx_gn */
              gsl_vector_memcpy(dx, state->dx_gn);
              gsl_vector_scale(dx, delta / state->norm_Dgn);
            }
          else
            {
              /* Cauchy point is inside, Gauss-Newton is outside trust region;
               * use double dogleg step */

              double beta = dogleg_beta(t, delta, trust_state->diag, state);

              /* compute: workp = t*dx_gn - dx_sd */
              scaled_addition(t, state->dx_gn, -1.0, state->dx_sd, state->workp);

              /* dx = dx_sd + beta*(t*dx_gn - dx_sd) */
              scaled_addition(beta, state->workp, 1.0, state->dx_sd, dx);
            }
        }
    }

  return GSL_SUCCESS;
}

static int
dogleg_preduction(const void * vtrust_state, const gsl_vector * dx,
                  double * pred, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  dogleg_state_t *state = (dogleg_state_t *) vstate;

  *pred = quadratic_preduction(trust_state->f, trust_state->J, dx, state->workn);

  return GSL_SUCCESS;
}

/*
dogleg_calc_gn()
  Calculate Gauss-Newton step which satisfies:

J dx_gn = -f

Inputs: trust_state - trust state variables
        dx          - (output) Gauss-Newton step

Return: success/error
*/

static int
dogleg_calc_gn(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx)
{
  int status;
  const gsl_multifit_nlinear_parameters *params = trust_state->params;

  /* initialize linear least squares solver */
  status = (params->solver->init)(trust_state, trust_state->solver_state);
  if (status)
    return status;

  /* prepare the linear solver to compute Gauss-Newton step */
  status = (params->solver->presolve)(0.0, trust_state, trust_state->solver_state);
  if (status)
    return status;

  /* solve: J dx_gn = -f for Gauss-Newton step */
  status = (params->solver->solve)(trust_state->f,
                                   dx,
                                   trust_state,
                                   trust_state->solver_state);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/*
dogleg_beta()
  This function finds beta in [0,1] such that the step

dx = dx_sd + beta*(t*dx_gn - dx_sd)

has norm

||D dx|| = delta

beta is the positive root of the quadratic:

a beta^2 + b beta + c = 0

with

a = ||D(t*dx_gn - dx_sd)||^2
b = 2 dx_sd^T D^T D (t*dx_gn - dx_sd)
c = ||D dx_sd||^2 - delta^2

Inputs: t     - amount of Gauss-Newton step to use for dogleg
                (= 1 for classical dogleg, <= 1 for double dogleg)
        delta - trust region radius
        diag  - diag(D) scaling matrix
        state - workspace
*/

static double
dogleg_beta(const double t, const double delta,
            const gsl_vector * diag, dogleg_state_t * state)
{
  double beta;
  double a, b, c;

  /* compute: workp = t*dx_gn - dx_sd */
  scaled_addition(t, state->dx_gn, -1.0, state->dx_sd, state->workp);

  /* a = || D (t*dx_gn - dx_sd) ||^2 */
  a = scaled_enorm(diag, state->workp);
  a *= a;

  /* workp = D^T D (t*dx_gn - dx_sd) */
  gsl_vector_mul(state->workp, diag);
  gsl_vector_mul(state->workp, diag);

  /* b = 2 dx_sd^T D^T D (t*dx_gn - dx-sd) */
  gsl_blas_ddot(state->dx_sd, state->workp, &b);
  b *= 2.0;

  /* c = || D dx_sd ||^2 - delta^2 = (||D dx_sd|| + delta) (||D dx_sd|| - delta) */
  c = (state->norm_Dsd + delta) * (state->norm_Dsd - delta);

  if (b > 0.0)
    {
      beta = (-2.0 * c) / (b + sqrt(b*b - 4.0*a*c));
    }
  else
    {
      beta = (-b + sqrt(b*b - 4.0*a*c)) / (2.0 * a);
    }

  return beta;
}

static const gsl_multifit_nlinear_trs dogleg_type =
{
  "dogleg",
  dogleg_alloc,
  dogleg_init,
  dogleg_preloop,
  dogleg_step,
  dogleg_preduction,
  dogleg_free
};

const gsl_multifit_nlinear_trs *gsl_multifit_nlinear_trs_dogleg = &dogleg_type;

static const gsl_multifit_nlinear_trs ddogleg_type =
{
  "double-dogleg",
  dogleg_alloc,
  dogleg_init,
  dogleg_preloop,
  dogleg_double_step,
  dogleg_preduction,
  dogleg_free
};

const gsl_multifit_nlinear_trs *gsl_multifit_nlinear_trs_ddogleg = &ddogleg_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* min/gsl_min.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 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.
 */

#ifndef __GSL_MIN_H__
#define __GSL_MIN_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
  {
    const char *name;
    size_t size;
    int (*set) (void *state, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper);
    int (*iterate) (void *state, gsl_function * f, double * x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper);
  }
gsl_min_fminimizer_type;

typedef struct
  {
    const gsl_min_fminimizer_type * type;
    gsl_function * function ;
    double x_minimum ;
    double x_lower ;
    double x_upper ;
    double f_minimum, f_lower, f_upper;
    void *state;
  }
gsl_min_fminimizer;

gsl_min_fminimizer *
gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * T) ;

void gsl_min_fminimizer_free (gsl_min_fminimizer * s);

int gsl_min_fminimizer_set (gsl_min_fminimizer * s,
                            gsl_function * f, double x_minimum,
                            double x_lower, double x_upper);

int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * s,
                                        gsl_function * f,
                                        double x_minimum, double f_minimum,
                                        double x_lower, double f_lower,
                                        double x_upper, double f_upper);

int gsl_min_fminimizer_iterate (gsl_min_fminimizer * s);

const char * gsl_min_fminimizer_name (const gsl_min_fminimizer * s);

double gsl_min_fminimizer_x_minimum (const gsl_min_fminimizer * s);
double gsl_min_fminimizer_x_lower (const gsl_min_fminimizer * s);
double gsl_min_fminimizer_x_upper (const gsl_min_fminimizer * s);
double gsl_min_fminimizer_f_minimum (const gsl_min_fminimizer * s);
double gsl_min_fminimizer_f_lower (const gsl_min_fminimizer * s);
double gsl_min_fminimizer_f_upper (const gsl_min_fminimizer * s);

/* Deprecated, use x_minimum instead */
double gsl_min_fminimizer_minimum (const gsl_min_fminimizer * s);

int
gsl_min_test_interval (double x_lower, double x_upper, double epsabs, double epsrel);

GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_goldensection;
GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_brent;
GSL_VAR const gsl_min_fminimizer_type  * gsl_min_fminimizer_quad_golden;

typedef
int (*gsl_min_bracketing_function)(gsl_function *f,
                                   double *x_minimum,double * f_minimum,
                                   double *x_lower, double * f_lower,
                                   double *x_upper, double * f_upper,
                                   size_t eval_max);

int
gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum,
                     double *x_lower, double * f_lower,
                     double *x_upper, double * f_upper,
                     size_t eval_max);

__END_DECLS

#endif /* __GSL_MIN_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
//
// Created by agibsonccc on 2/6/16.
//

#ifndef NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H
#define NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H
#include <cblas.h>

/*
enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102 };
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113,
    AtlasConj=114};
enum CBLAS_UPLO  {CblasUpper=121, CblasLower=122};
enum CBLAS_DIAG  {CblasNonUnit=131, CblasUnit=132};
enum CBLAS_SIDE  {CblasLeft=141, CblasRight=142};
*/
#include <dll.h>

#ifdef __cplusplus
extern "C" {
#endif
/**
 * Converts a character
 * to its proper enum
 * for row (c) or column (f) ordering
 * default is row major
 */
CBLAS_ORDER convertOrder(int from);
/**
 * Converts a character to its proper enum
 * t -> transpose
 * n -> no transpose
 * c -> conj
 */
CBLAS_TRANSPOSE convertTranspose(int from);
/**
 * Upper or lower
 * U/u -> upper
 * L/l -> lower
 *
 * Default is upper
 */
CBLAS_UPLO convertUplo(int from);

/**
 * For diagonals:
 * u/U -> unit
 * n/N -> non unit
 *
 * Default: unit
 */
CBLAS_DIAG convertDiag(int from);
/**
 * Side of a matrix, left or right
 * l /L -> left
 * r/R -> right
 * default: left
 */
CBLAS_SIDE convertSide(int from);

#ifdef __cplusplus
}
#endif

#endif //NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H


text:
/* ieee-utils/fp-os2.c
 *
 * Copyright (C) 2001 Henry Sobotka <sobotka@axess.com>
 *
 * 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.
 */

#include <float.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  unsigned mode = 0;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      _control87(PC_24, MCW_PC);
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      _control87(PC_53, MCW_PC);
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      _control87(PC_64, MCW_PC);
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      _control87(RC_NEAR, MCW_RC);
      break ;
    case GSL_IEEE_ROUND_DOWN:
      _control87(RC_DOWN, MCW_RC);
      break ;
    case GSL_IEEE_ROUND_UP:
      _control87(RC_UP, MCW_RC);
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      _control87(RC_CHOP, MCW_RC);
      break ;
    default:
      _control87(RC_NEAR, MCW_RC);
    }

  /* Turn on all the exceptions apart from 'inexact' */

  mode = EM_INVALID | EM_DENORMAL | EM_ZERODIVIDE | EM_OVERFLOW | EM_UNDERFLOW;

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode &= ~ EM_INVALID;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    mode &= ~ EM_DENORMAL;

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode &= ~ EM_ZERODIVIDE;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode &= ~ EM_OVERFLOW;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode &= ~ EM_UNDERFLOW;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      mode |= EM_INEXACT;
    }
  else
    {
      mode &= ~ EM_INEXACT;
    }

  _control87(mode, MCW_EM);

  return GSL_SUCCESS ;
}


text:
/* randist/weibull.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The Weibull distribution has the form,

   p(x) dx = (b/a) (x/a)^(b-1) exp(-(x/a)^b) dx

 */

double
gsl_ran_weibull (const gsl_rng * r, const double a, const double b)
{
  double x = gsl_rng_uniform_pos (r);

  double z = pow (-log (x), 1 / b);

  return a * z;
}

double
gsl_ran_weibull_pdf (const double x, const double a, const double b)
{
  if (x < 0)
    {
      return 0 ;
    }
  else if (x == 0)
    {
      if (b == 1)
        return 1/a ;
      else
        return 0 ;
    }
  else if (b == 1)
    {
      return exp(-x/a)/a ;
    }
  else
    {
      double p = (b/a) * exp (-pow (x/a, b) + (b - 1) * log (x/a));
      return p;
    }
}