text:
#ifndef BLAS_WRAPPER_H
#define BLAS_WRAPPER_H

#include <cassert>
#include <cmath>
#include <cstring>
#include <vector>

// Useful dense kernels from BLAS, with readable, overloaded, cross-platform names and some simplified calling
//  dot    (dot-product of vectors)
//  nrm2   (2-norm of a vector)
//  asum   (1-norm of a vector)
//  amax   (index of maximum absolute value in a vector)
//  swap   (exchanging values in two vectors)
//  copy   (copying values from one vector to another)
//  axpy   (adding a scalar times a vector to another vector)
//  scal   (multiplying a vector by a scalar)
//  gemv   (multiplying a matrix times a vector, scaling, and adding result to another vector))
//  gemm   (multiplying two matrices, scaling, and adding result to another matrix)
// In addition:
//  set_zero  (zero out all entries in a vector)
//  abs_max   (return the infinity norm of a vector, i.e. the magnitude of its largest element)
// There are also version using std::vector for convenience.

// Matrices are always assumed to be in column-major format.

// You can #define one of:
//    USE_FORTRAN_BLAS (if your BLAS calls should look like dgemm_ with FORTRAN calling conventions, as in GOTO BLAS)
//    USE_AMD_BLAS     (if using the AMD Math Library)
//    USE_CBLAS        (if instead you have calls like cblas_dgemm, and have a file "cblas.h" available)
// or, if you're on the Mac, it will default to the vecLib CBLAS if none of these are specified.

namespace BLAS{

    template<class T>
    inline void set_zero(int n, T *x)
    { std::memset(x, 0, n*sizeof(T)); }

}

//============================================================================
#ifdef USE_FORTRAN_BLAS

extern "C" {
    double dsdot_(const int*, const float*, const int*, const float*, const int*);
    double sdot_(const int*, const float*, const int*, const float*, const int*);
    double ddot_(const int*, const double*, const int*, const double*, const int*);
    float snrm2_(const int*, const float*, const int*);
    double dnrm2_(const int*, const double*, const int*);
    float sasum_(const int*, const float*, const int*);
    double dasum_(const int*, const double*, const int*);
    int isamax_(const int*, const float*, const int*);
    int idamax_(const int*, const double*, const int*);
    void sswap_(const int*, float*, const int*, float*, const int*);
    void dswap_(const int*, double*, const int*, double*, const int*);
    void scopy_(const int*, const float*, const int*, float*, const int*);
    void dcopy_(const int*, const double*, const int*, double*, const int*);
    void saxpy_(const int*, const float*, const float*, const int*, float*, const int*);
    void daxpy_(const int*, const double*, const double*, const int*, double*, const int*);
    void sscal_(const int*, const float*, float*, const int*);
    void dscal_(const int*, const double*, double*, const int*);
    void sgemv_(const char*, const int*, const int*, const float*, const float*, const int*, const float*, const int*, const float*, float*, const int*);
    void dgemv_(const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*);
    void sgemm_(const char*, const char*, const int*, const int*, const int*, const float*, const float*, const int*, const float*, const int*, const float*, float*, const int*);
    void dgemm_(const char*, const char*, const int*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*);
}

namespace BLAS{

    enum Transpose {NoTrans='N', Trans='T'};
    enum UpperLower {Upper='U', Lower='L'};
    enum UnitDiag {NonUnit='N', Unit='U'};
    enum Side {Left='L', Right='R'};

    // dot products

    inline double dot(int n, const float *x, int incx, const float *y, int incy=1)
    { return dsdot_(&n, x, &incx, y, &incy); }

    inline double dot(int n, const float *x, const float *y, int incy=1)
    { const int one=1; return dsdot_(&n, x, &one, y, &incy); }

    inline float dotf(int n, const float *x, int incx, const float *y, int incy=1)
    { return (float)sdot_(&n, x, &incx, y, &incy); }

    inline float dotf(int n, const float *x, const float *y, int incy=1)
    { const int one=1; return (float)sdot_(&n, x, &one, y, &incy); }

    inline double dot(int n, const double *x, int incx, const double *y, int incy=1)
    { return ddot_(&n, x, &incx, y, &incy); }

    inline double dot(int n, const double *x, const double *y, int incy=1)
    { const int one=1; return ddot_(&n, x, &one, y, &incy); }

    // 2-norm

    inline float norm2(int n, const float *x, int incx=1)
    { return snrm2_(&n, x, &incx); }

    inline double norm2(int n, const double *x, int incx=1)
    { return dnrm2_(&n, x, &incx); }

    // 1-norm (sum of absolute values)

    inline float abs_sum(int n, const float *x, int incx=1)
    { return sasum_(&n, x, &incx); }

    inline double abs_sum(int n, const double *x, int incx=1)
    { return dasum_(&n, x, &incx); }

    // inf-norm (maximum absolute value: index of max returned)

    inline int index_abs_max(int n, const float *x, int incx=1)
    { return isamax_(&n, x, &incx)-1; }

    inline int index_abs_max(int n, const double *x, int incx=1)
    { return idamax_(&n, x, &incx)-1; }

    inline float abs_max(int n, const float *x, int incx=1)
    { return std::fabs(x[isamax_(&n, x, &incx)-1]); }

    inline double abs_max(int n, const double *x, int incx=1)
    { return std::fabs(x[idamax_(&n, x, &incx)-1]); }

    // swap (actual data exchanged, not just pointers)

    inline void swap(int n, float *x, int incx, float *y, int incy=1)
    { sswap_(&n, x, &incx, y, &incy); }

    inline void swap(int n, float *x, float *y, int incy=1)
    { const int one=1; sswap_(&n, x, &one, y, &incy); }

    inline void swap(int n, double *x, int incx, double *y, int incy=1)
    { dswap_(&n, x, &incx, y, &incy); }

    inline void swap(int n, double *x, double *y, int incy=1)
    { const int one=1; dswap_(&n, x, &one, y, &incy); }

    // copy (y=x)

    inline void copy(int n, const float *x, int incx, float *y, int incy=1)
    { scopy_(&n, x, &incx, y, &incy); }

    inline void copy(int n, const float *x, float *y, int incy=1)
    { const int one=1; scopy_(&n, x, &one, y, &incy); }

    inline void copy(int n, const double *x, int incx, double *y, int incy=1)
    { dcopy_(&n, x, &incx, y, &incy); }

    inline void copy(int n, const double *x, double *y, int incy=1)
    { const int one=1; dcopy_(&n, x, &one, y, &incy); }

    // saxpy (y=alpha*x+y)

    inline void add_scaled(int n, float alpha, const float *x, int incx, float *y, int incy=1)
    { saxpy_(&n, &alpha, x, &incx, y, &incy); }

    inline void add_scaled(int n, float alpha, const float *x, float *y, int incy=1)
    { const int one=1; saxpy_(&n, &alpha, x, &one, y, &incy); }

    inline void add_scaled(int n, double alpha, const double *x, int incx, double *y, int incy=1)
    { daxpy_(&n, &alpha, x, &incx, y, &incy); }

    inline void add_scaled(int n, double alpha, const double *x, double *y, int incy=1)
    { const int one=1; daxpy_(&n, &alpha, x, &one, y, &incy); }

    // scale (x=alpha*x)

    inline void scale(int n, float alpha, float *x, int incx=1)
    { sscal_(&n, &alpha, x, &incx); }

    inline void scale(int n, double alpha, double *x, int incx=1)
    { dscal_(&n, &alpha, x, &incx); }

    // gemv (y=alpha*A*x+beta*y, or using A^T)
    // The matrix is always m*n; the size of x and y depend on if A is transposed.

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, float alpha, const float *A, int lda,
                                       const float *x, int incx, float beta, float *y, int incy=1)
    { sgemv_((const char*)&transpose, &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }

    inline void multiply_matrix_vector(int m, int n, const float *A, const float *x, float *y, int incy=1) // y=A*x
    { const int onei=1; const float zero=0, onef=1; sgemv_("N", &m, &n, &onef, A, &m, x, &onei, &zero, y, &incy); }

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, double alpha, const double *A, int lda,
                                       const double *x, int incx, double beta, double *y, int incy=1)
    { dgemv_((const char*)&transpose, &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }

    inline void multiply_matrix_vector(int m, int n, const double *A, const double *x, double *y, int incy=1) // y=A*x
    { const int onei=1; const double zero=0, onef=1; dgemv_("N", &m, &n, &onef, A, &m, x, &onei, &zero, y, &incy); }

    // gemm (C=alpha*A*B+beta*C)

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, float alpha, const float *A, int lda,
                                       const float *B, int ldb, float beta, float *C, int ldc)
    { sgemm_((const char*)&transA, (const char*)&transB, &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const float *A, const float *B, float *C)
    { const float zero=0, one=1; sgemm_("N", "N", &m, &n, &k, &one, A, &m, B, &k, &zero, C, &m); } // C=A*B

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, double alpha, const double *A, int lda,
                                       const double *B, int ldb, double beta, double *C, int ldc)
    { dgemm_((const char*)&transA, (const char*)&transB, &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const double *A, const double *B, double *C)
    { const double zero=0, one=1; dgemm_("N", "N", &m, &n, &k, &one, A, &m, B, &k, &zero, C, &m); } // C=A*B

};

//============================================================================
#elif defined USE_AMD_BLAS

#include <acml.h>

namespace BLAS{

    enum Transpose {NoTrans='N', Trans='T'};
    enum UpperLower {Upper='U', Lower='L'};
    enum UnitDiag {NonUnit='N', Unit='U'};
    enum Side {Left='L', Right='R'};

    // dot products

    inline double dot(int n, const float *x, int incx, const float *y, int incy=1)
    { return dsdot(n, (float*)x, incx, (float*)y, incy); }

    inline double dot(int n, const float *x, const float *y, int incy=1)
    { return dsdot(n, (float*)x, 1, (float*)y, incy); }

    inline float dotf(int n, const float *x, int incx, const float *y, int incy=1)
    { return sdot(n, (float*)x, incx, (float*)y, incy); }

    inline float dotf(int n, const float *x, const float *y, int incy=1)
    { return sdot(n, (float*)x, 1, (float*)y, incy); }

    inline double dot(int n, const double *x, int incx, const double *y, int incy=1)
    { return ddot(n, (double *)x, incx, (double *)y, incy); }

    inline double dot(int n, const double *x, const double *y, int incy=1)
    { return ddot(n, (double *)x, 1, (double *)y, incy); }

    // 2-norm

    inline float norm2(int n, const float *x, int incx=1)
    { return snrm2(n, (float*)x, incx); }

    inline double norm2(int n, const double *x, int incx=1)
    { return dnrm2(n, (double*)x, incx); }

    // 1-norm (sum of absolute values)

    inline float abs_sum(int n, const float *x, int incx=1)
    { return sasum(n, (float*)x, incx); }

    inline double abs_sum(int n, const double *x, int incx=1)
    { return dasum(n, (double*)x, incx); }

    // inf-norm (maximum absolute value: index of max returned)

    inline int index_abs_max(int n, const float *x, int incx=1)
    { return isamax(n, (float*)x, incx)-1; }

    inline int index_abs_max(int n, const double *x, int incx=1)
    { return idamax(n, (double*)x, incx)-1; }

    inline float abs_max(int n, const float *x, int incx=1)
    { return std::fabs(x[isamax(n, (float*)x, incx)]-1); }

    inline double abs_max(int n, const double *x, int incx=1)
    { return std::fabs(x[idamax(n, (double*)x, incx)]-1); }

    // swap (actual data exchanged, not just pointers)

    inline void swap(int n, float *x, int incx, float *y, int incy=1)
    { sswap(n, x, incx, y, incy); }

    inline void swap(int n, float *x, float *y, int incy=1)
    { sswap(n, x, 1, y, incy); }

    inline void swap(int n, double *x, int incx, double *y, int incy=1)
    { dswap(n, x, incx, y, incy); }

    inline void swap(int n, double *x, double *y, int incy=1)
    { dswap(n, x, 1, y, incy); }

    // copy (y=x)

    inline void copy(int n, const float *x, int incx, float *y, int incy=1)
    { scopy(n, (float*)x, incx, y, incy); }

    inline void copy(int n, const float *x, float *y, int incy=1)
    { scopy(n, (float*)x, 1, y, incy); }

    inline void copy(int n, const double *x, int incx, double *y, int incy=1)
    { dcopy(n, (double *)x, incx, y, incy); }

    inline void copy(int n, const double *x, double *y, int incy=1)
    { dcopy(n, (double *)x, 1, y, incy); }

    // saxpy (y=alpha*x+y)

    inline void add_scaled(int n, float alpha, const float *x, int incx, float *y, int incy=1)
    { saxpy(n, alpha, (float*)x, incx, y, incy); }

    inline void add_scaled(int n, float alpha, const float *x, float *y, int incy=1)
    { saxpy(n, alpha, (float*)x, 1, y, incy); }

    inline void add_scaled(int n, double alpha, const double *x, int incx, double *y, int incy=1)
    { daxpy(n, alpha, (double*)x, incx, y, incy); }

    inline void add_scaled(int n, double alpha, const double *x, double *y, int incy=1)
    { daxpy(n, alpha, (double*)x, 1, y, incy); }

    // scale (x=alpha*x)

    inline void scale(int n, float alpha, float *x, int incx=1)
    { sscal(n, alpha, x, incx); }

    inline void scale(int n, double alpha, double *x, int incx=1)
    { dscal(n, alpha, x, incx); }

    // gemv (y=alpha*A*x+beta*y, or using A^T)
    // The matrix is always m*n; the size of x and y depend on if A is transposed.

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, float alpha, const float *A, int lda,
                                       const float *x, int incx, float beta, float *y, int incy=1)
    { sgemv(transpose, m, n, alpha, (float*)A, lda, (float*)x, incx, beta, y, incy); }

    inline void multiply_matrix_vector(int m, int n, const float *A, const float *x, float *y, int incy=1) // y=A*x
    { sgemv(NoTrans, m, n, 1.f, (float*)A, m, (float*)x, 1, 0.f, y, incy); }

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, double alpha, const double *A, int lda,
                                       const double *x, int incx, double beta, double *y, int incy=1)
    { dgemv(transpose, m, n, alpha, (double*)A, lda, (double*)x, incx, beta, y, incy); }

    inline void multiply_matrix_vector(int m, int n, const double *A, const double *x, double *y, int incy=1) // y=A*x
    { dgemv(NoTrans, m, n, 1., (double*)A, m, (double*)x, 1, 0., y, incy); }

    // gemm (C=alpha*A*B+beta*C)

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, float alpha, const float *A, int lda,
                                       const float *B, int ldb, float beta, float *C, int ldc)
    { sgemm(transA, transB, m, n, k, alpha, (float*)A, lda, (float*)B, ldb, beta, C, ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const float *A, const float *B, float *C)
    { sgemm(NoTrans, NoTrans, m, n, k, 1.f, (float*)A, m, (float*)B, k, 0.f, C, m); } // C=A*B

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, double alpha, const double *A, int lda,
                                       const double *B, int ldb, double beta, double *C, int ldc)
    { dgemm(transA, transB, m, n, k, alpha, (double*)A, lda, (double*)B, ldb, beta, C, ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const double *A, const double *B, double *C)
    { dgemm(NoTrans, NoTrans, m, n, k, 1., (double*)A, m, (double*)B, k, 0., C, m); } // C=A*B

};

//============================================================================
#elif defined USE_CBLAS || defined __APPLE__

#ifdef USE_CBLAS
#include <cblas.h>
#elif defined __APPLE__
#include <Accelerate/Accelerate.h>
#endif

namespace BLAS{

    enum Transpose {NoTrans=CblasNoTrans, Trans=CblasTrans};
    enum UpperLower {Upper=CblasUpper, Lower=CblasLower};
    enum UnitDiag {NonUnit=CblasNonUnit, Unit=CblasUnit};
    enum Side {Left=CblasLeft, Right=CblasRight};

    // dot products

    inline float dotf(int n, const float *x, int incx, const float *y, int incy=1)
    { return cblas_sdot(n, x, incx, y, incy); }

    inline float dotf(int n, const float *x, const float *y, int incy=1)
    { return cblas_sdot(n, x, 1, y, incy); }

    inline double dot(int n, const float *x, int incx, const float *y, int incy=1)
    { return cblas_dsdot(n, x, incx, y, incy); }

    inline double dot(int n, const float *x, const float *y, int incy=1)
    { return cblas_dsdot(n, x, 1, y, incy); }

    inline double dot(int n, const double *x, int incx, const double *y, int incy=1)
    { return cblas_ddot(n, x, incx, y, incy); }

    inline double dot(int n, const double *x, const double *y, int incy=1)
    { return cblas_ddot(n, x, 1, y, incy); }

    // 2-norm

    inline float norm2(int n, const float *x, int incx=1)
    { return cblas_snrm2(n, x, incx); }

    inline double norm2(int n, const double *x, int incx=1)
    { return cblas_dnrm2(n, x, incx); }

    // 1-norm (sum of absolute values)

    inline float abs_sum(int n, const float *x, int incx=1)
    { return cblas_sasum(n, x, incx); }

    inline double abs_sum(int n, const double *x, int incx=1)
    { return cblas_dasum(n, x, incx); }

    // inf-norm (maximum absolute value)

    inline int index_abs_max(int n, const float *x, int incx=1)
    { return cblas_isamax(n, x, incx); }

    inline int index_abs_max(int n, const double *x, int incx=1)
    { return cblas_idamax(n, x, incx); }

    inline float abs_max(int n, const float *x, int incx=1)
    { return std::fabs(x[cblas_isamax(n, x, incx)]); }

    inline double abs_max(int n, const double *x, int incx=1)
    { return std::fabs(x[cblas_idamax(n, x, incx)]); }

    // swap (actual data exchanged, not just pointers)

    inline void swap(int n, float *x, int incx, float *y, int incy=1)
    { cblas_sswap(n, x, incx, y, incy); }

    inline void swap(int n, float *x, float *y, int incy=1)
    { cblas_sswap(n, x, 1, y, incy); }

    inline void swap(int n, double *x, int incx, double *y, int incy=1)
    { cblas_dswap(n, x, incx, y, incy); }

    inline void swap(int n, double *x, double *y, int incy=1)
    { cblas_dswap(n, x, 1, y, incy); }

    // copy (y=x)

    inline void copy(int n, const float *x, int incx, float *y, int incy=1)
    { cblas_scopy(n, x, incx, y, incy); }

    inline void copy(int n, const float *x, float *y, int incy=1)
    { cblas_scopy(n, x, 1, y, incy); }

    inline void copy(int n, const double *x, int incx, double *y, int incy=1)
    { cblas_dcopy(n, x, incx, y, incy); }

    inline void copy(int n, const double *x, double *y, int incy=1)
    { cblas_dcopy(n, x, 1, y, incy); }

    // saxpy (y=alpha*x+y)

    inline void add_scaled(int n, float alpha, const float *x, int incx, float *y, int incy=1)
    { cblas_saxpy(n, alpha, x, incx, y, incy); }

    inline void add_scaled(int n, float alpha, const float *x, float *y, int incy=1)
    { cblas_saxpy(n, alpha, x, 1, y, incy); }

    inline void add_scaled(int n, double alpha, const double *x, int incx, double *y, int incy=1)
    { cblas_daxpy(n, alpha, x, incx, y, incy); }

    inline void add_scaled(int n, double alpha, const double *x, double *y, int incy=1)
    { cblas_daxpy(n, alpha, x, 1, y, incy); }

    // scale (x=alpha*x)

    inline void scale(int n, float alpha, float *x, int incx=1)
    { cblas_sscal(n, alpha, x, incx); }

    inline void scale(int n, double alpha, double *x, int incx=1)
    { cblas_dscal(n, alpha, x, incx); }

    // gemv (y=alpha*A*x+beta*y, or using A^T)
    // The matrix is always m*n; the size of x and y depend on if A is transposed.

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, float alpha, const float *A, int lda,
                                       const float *x, int incx,
                                       float beta, float *y, int incy=1)
    { cblas_sgemv(CblasColMajor, (CBLAS_TRANSPOSE)transpose, m, n, alpha, A, lda, x, incx, beta, y, incy); }

    inline void multiply_matrix_vector(int m, int n, const float *A, const float *x, float *y, int incy=1) // y=A*x
    { cblas_sgemv(CblasColMajor, CblasNoTrans, m, n, 1.f, A, m, x, 1, 0.f, y, incy); }

    inline void multiply_matrix_vector(Transpose transpose,
                                       int m, int n, double alpha, const double *A, int lda,
                                       const double *x, int incx, double beta, double *y, int incy=1)
    { cblas_dgemv(CblasColMajor, (CBLAS_TRANSPOSE)transpose, m, n, alpha, A, lda, x, incx, beta, y, incy); }

    inline void multiply_matrix_vector(int m, int n, const double *A, const double *x, double *y, int incy=1) // y=A*x
    { cblas_dgemv(CblasColMajor, CblasNoTrans, m, n, 1., A, m, x, 1, 0., y, incy); }

    // gemm (C=alpha*A*B+beta*C)

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, float alpha, const float *A, int lda,
                                       const float *B, int ldb, float beta, float *C, int ldc)
    { cblas_sgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const float *A, const float *B, float *C)
    { cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.f, A, m, B, k, 0.f, C, m); } // C=A*B

    inline void multiply_matrix_matrix(Transpose transA, Transpose transB,
                                       int m, int n, int k, double alpha, const double *A, int lda,
                                       const double *B, int ldb, double beta, double *C, int ldc)
    { cblas_dgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc); }

    inline void multiply_matrix_matrix(int m, int n, int k, const double *A, const double *B, double *C)
    { cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1., A, m, B, k, 0., C, m); } // C=A*B

    inline void rank_one_update( int m, int n, double alpha, const double* x, const double* y, double* A )
    {
        cblas_dger( CblasColMajor, m, n, alpha, x, 1, y, 1, A, m );
    }

}  // namespace BLAS

#endif

// std::vector calls =========================================================
namespace BLAS{

    template<class T>
    inline void set_zero(std::vector<T> &x)
    { set_zero((int)x.size(), &x[0]); }

    inline float dotf(const std::vector<float> &x, const std::vector<float> &y)
    { assert(x.size()==y.size()); return dotf((int)x.size(), &x[0], &y[0]); }

    inline double dot(const std::vector<float> &x, const std::vector<float> &y)
    { assert(x.size()==y.size()); return dot((int)x.size(), &x[0], &y[0]); }

    inline double dot(const std::vector<double> &x, const std::vector<double> &y)
    { assert(x.size()==y.size()); return dot((int)x.size(), &x[0], &y[0]); }

    inline float norm2(const std::vector<float> &x)
    { return norm2((int)x.size(), &x[0]); }

    inline double norm2(const std::vector<double> &x)
    { return norm2((int)x.size(), &x[0]); }

    inline float abs_sum(const std::vector<float> &x)
    { return abs_sum((int)x.size(), &x[0]); }

    inline double abs_sum(const std::vector<double> &x)
    { return abs_sum((int)x.size(), &x[0]); }

    inline int index_abs_max(const std::vector<float> &x)
    { return index_abs_max((int)x.size(), &x[0]); }

    inline int index_abs_max(const std::vector<double> &x)
    { return index_abs_max((int)x.size(), &x[0]); }

    inline float abs_max(const std::vector<float> &x)
    { return abs_max((int)x.size(), &x[0]); }

    inline double abs_max(const std::vector<double> &x)
    { return abs_max((int)x.size(), &x[0]); }

    inline void swap(std::vector<float> &x, std::vector<float> &y)
    { assert(x.size()==y.size()); swap((int)x.size(), &x[0], &y[0]); }

    inline void swap(std::vector<double> &x, std::vector<double> &y)
    { assert(x.size()==y.size()); swap((int)x.size(), &x[0], &y[0]); }

    inline void copy(const std::vector<float> &x, std::vector<float> &y)
    { assert(x.size()==y.size()); copy((int)x.size(), &x[0], &y[0]); }

    inline void copy(const std::vector<double> &x, std::vector<double> &y)
    { assert(x.size()==y.size()); copy((int)x.size(), &x[0], &y[0]); }

    inline void add_scaled(float alpha, const std::vector<float> &x, std::vector<float> &y)
    { assert(x.size()==y.size()); add_scaled((int)x.size(), alpha, &x[0], &y[0]); }

    inline void add_scaled(double alpha, const std::vector<double> &x, std::vector<double> &y)
    { assert(x.size()==y.size()); add_scaled((int)x.size(), alpha, &x[0], &y[0]); }

    inline void scale(float alpha, std::vector<float> &x)
    { scale((int)x.size(), alpha, &x[0]); }

    inline void scale(float alpha, std::vector<double> &x)
    { scale((int)x.size(), alpha, &x[0]); }

    // I'm not sure if it makes sense to include level 2 or level 3 std::vector versions,
    // since there isn't an STL matrix type...

}  // namespace BLAS

#endif


text:
/* -*- mode: C; c-basic-offset: 4 -*- */
/* ex: set shiftwidth=4 tabstop=4 expandtab: */
/*
 * Copyright (c) 2013, Georgia Tech Research Corporation
 * All rights reserved.
 *
 * Author(s): Neil T. Dantam <ntd@gatech.edu>
 * Georgia Tech Humanoid Robotics Lab
 * Under Direction of Prof. Mike Stilman <mstilman@cc.gatech.edu>
 *
 *
 * This file is provided under the following "BSD-style" License:
 *
 *
 *   Redistribution and use in source and binary forms, with or
 *   without modification, are permitted provided that the following
 *   conditions are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 *   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 *   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 *   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 *   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 *   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *   POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include <amino.h>
#include <cblas.h>
#include "reflex.h"

void x_fopen( const char *prefix, const char *suffix, size_t n, FILE *f[] ) {
    for( size_t i = 0; i < n; i ++ ) {
        char buf[512];
        sprintf(buf,"%s%lu%s", prefix, i, suffix);
        f[i] = fopen(buf, "w");
    }
}

void x_write( double t, size_t n, double e[], FILE *f[7] ) {
    for( size_t i = 0; i < n; i ++ ) {
        fprintf(f[i], "%f %f\n", t, e[i] );
        //fprintf(stdout, "%f %f\n", t, e[i] );
    }
}

void e_corrupt( double theta_max, double x_max, const double e0[7], double e1[7] )
{
    // rotation
    double aa[4];
    aa_vrand( 4, aa );
    aa_la_normalize(3,aa);
    aa[3] *= theta_max;

    // translation
    double v[3];
    aa_vrand( 3, v );
    for( size_t i = 0; i < 3; i ++ ) {
        v[i] *= x_max;
    }

    // convert
    double ec[7];
    aa_tf_axang2quat(aa, ec+AA_TF_QUTR_Q);
    AA_MEM_CPY(ec+AA_TF_QUTR_V, v, 3 );

    // mul
    aa_tf_qutr_mul( e0, ec, e1 );

}

rfx_tf_dx XX_true;
rfx_tf_dx XX_est;
rfx_tf_dx ZZ;
rfx_tf_dx UU;

double P[16*16] = {0};
double V[16*16] = {0};
double W[8*8] = {0};

double Pa[13*13] = {0};
double Wa[13*13] = {0};
double Va[13*13] = {0};

double Pb[13*13] = {0};
double Wb[7*7] = {0};
double Vb[13*13] = {0};

int main(void)
{

    double dt = .01;

    FILE *fout_z[7];
    FILE *fout_x_true[7];
    FILE *fout_x_est[7];

    FILE *fout_dx_true[6];
    FILE *fout_dx_est[6];

    // state
    memset(&XX_true,0,sizeof(XX_true));
    memset(&XX_est,0,sizeof(XX_est));
    memset(&ZZ,0,sizeof(ZZ));
    memset(&UU,0,sizeof(UU));

    aa_la_diag( 16, P, 1.0 );
    aa_la_diag( 8, W, 1.0 );
    aa_la_diag( 16, V, 1.0e-2 );

    for( size_t i = 0; i < 8; i ++ ) {
        AA_MATREF(V,16, i, 8+i ) = 1e-2;
        AA_MATREF(V,16, 8+i, i ) = 1e-2;
    }

    aa_la_diag( 13, Pb, 1 );
    aa_la_diag( 7, Wb, dt*5e-1 );

    aa_la_diag( 13, Pa, 1.0e-2 );
    aa_la_diag( 13, Wa, 1.0 );
    aa_la_diag( 13, Va, 1.0e2 );

    // files

    x_fopen("xtrue", ".dat", 7,fout_x_true);
    x_fopen("z", ".dat", 7,fout_z);
    x_fopen("xest", ".dat", 7,fout_x_est);
    x_fopen("dx_true", ".dat", 6,fout_dx_true);
    x_fopen("dx_est", ".dat", 6,fout_dx_est);

    memcpy(XX_true.tf.r.data, aa_tf_quat_ident, 4*sizeof(XX_true.tf.r.data[0]));
    memcpy(XX_est.tf.r.data, aa_tf_quat_ident, 4*sizeof(XX_true.tf.r.data[0]));

    for( double t = 0.0; t < 10.0; t += dt ) {
        fprintf(stderr, "t: %f, ", t );
        // pristine velocity
        double dx[6] = {cos(t*M_PI / 1), 0, 0,
                        0, 0, sin(t*M_PI / 1)};
        AA_MEM_CPY(XX_true.dx.data, dx, 6);

        // integrate true
        {
            double e_next[7];
            aa_tf_qutr_svel(XX_true.tf.data, XX_true.dx.data, dt, e_next );
            memcpy( XX_true.tf.data, e_next, sizeof(e_next) );
        }

        // corrupt measurement
        e_corrupt( 10*M_PI/180, 7e-2, XX_true.tf.data, ZZ.tf.data );

        // filter
        aa_tick("filter time: ");
        double S[8];
        double Sz[8];
        double dS_est[8];
        double dS_true[8];
        //double dE_est[7];
        aa_tf_qutr2duqu( XX_est.tf.data, S);
        aa_tf_qutr2duqu( ZZ.tf.data, Sz);
        aa_tf_duqu_vel2diff( S, XX_est.dx.data, dS_est );
        aa_tf_duqu_vel2diff( S, XX_true.dx.data, dS_true );

        //aa_tf_qutr_vel2diff( XX_est.tf.data, XX_est.dx.data, dE_est );

        /* rfx_lqg_duqu_predict( dt, S, dS_est, P, V ); */
        /* rfx_lqg_duqu_correct( 1, */
        /*                       S, dS_est, */
        /*                       Sz, /\*ZZ.dx.data, *\/ */
        /*                       P, W ); */
        /* aa_tf_duqu_diff2vel( S, dS_est, XX_est.dx.data ); */
        /* aa_tf_duqu2qutr( S, XX_est.tf.data); */

        rfx_lqg_qutr_process_noise( dt, 1, 1, XX_est.tf.data, Vb );

        rfx_lqg_qutr_predict( dt, XX_est.tf.data, XX_est.dx.data, Pb, Vb );
        rfx_lqg_qutr_correct( 1,
                              XX_est.tf.data, XX_est.dx.data,
                              ZZ.tf.data,
                              Pb, Wb );
        //aa_tf_qutr_diff2vel( XX_est.tf.data, dE_est, XX_est.dx.data );

        //rfx_tf_filter_update_work( dt, XX_est.tf.data, UU.tf.data, ZZ.tf.data, Pa, Va, Wa );
        aa_tock();
        //printf("--\n");
        //printf("dS_true: "); aa_dump_vec( stdout, dS_true, 8 );
        //printf("dS_est: "); aa_dump_vec(  stdout, dS_est, 8 );
        //aa_dump_mat( stdout, P, 16, 16 );

        //return 0;
        // print
        x_write( t, 7, XX_true.tf.data, fout_x_true );
        x_write( t, 6, XX_true.dx.data, fout_dx_true );

        x_write( t, 7, XX_est.tf.data, fout_x_est );
        x_write( t, 6, XX_est.dx.data, fout_dx_est );

        x_write( t, 7, ZZ.tf.data, fout_z );
    }

}

ext: c

max_stars_repo_name: golems/reflex

max_issues_repo_name: golems/reflex

max_forks_repo_name: golems/reflex


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_sf_trig.h>
#include <Accelerate/Accelerate.h>
#include "SE_trans.h"
#include "sigma.h"

#define ALPHA 0.99
#define DOMD 3.14

double* AllocVec(int dim)
{
  double* vec;
  vec = (double*)calloc(dim, sizeof(double));

  /* Check memory */
  if (vec == NULL) {
    free(vec);
    fprintf(stderr,"Vector was not allocated! \n");
    exit(1);
  }

  return vec; /* Return the first address of the allocated memory */
}

void FreeVec(double* v)
{
  free(v);
}

double f(double t)
{
  return 0.25 * ((1+t) * log1p(t) + (1-t) * log1p(-t) - 2 * M_LN2) / M_LN2;
}

double G(double x)
{
  return 0.25 * x * SE_trans_div(-1, 1, x) / M_LN2;
}

double Fapp(double* p, double t, int m)
{
  double h = sqrt(M_PI*DOMD / (ALPHA*m));

  int j;

  double sum = 0;

  for (j = -m; j <= m; j++) {
    sum += p[j+m] * gsl_sf_sinc((SE_trans_inv(-1, 1, t) - j*h) / h);
  }

  return sum;
}

double eta(double x)
{
  return 0.5*(1 + x);
}

int main()
{
  int i, j;
  int n, N;
  double t, h;
  double err, maxerr;
  double Iast;
  clock_t start, end;
  double time;
  double *J, *g, *c;

  int STEP = 1000;

  for (n = 3; n <= 180; n += 6) {
    start = clock();

    N = 2*n+1;
    J = AllocVec(N*N);
    g = AllocVec(N);
    c = AllocVec(N);

    h = sqrt(M_PI*DOMD / (ALPHA*n));

    Iast = 0;
    for (i = -n; i <= n; i++) {
      Iast += G(i*h);
    }
    Iast *= h;

    for (i = -n; i <= n; i++) {
      for (j = -n; j <= n; j++) {
        J[(i+n)*N + (j+n)] = 0.5 + (double)sigma_orig[1000 + (i - j)];
      }
      g[(i+n)] = G(i*h) - 0.5 * Iast * SE_trans_div(-1, 1, i*h);
    }

    /* c = h J g */
    cblas_dgemv(CblasRowMajor, CblasNoTrans, N, N, h, J, N, g, 1, 0.0, c, 1);

    maxerr = 0;
    for (i = - STEP + 1; i <= STEP - 1; i++) {
      t = (double)i / (double)(STEP);

      err = fabs(f(t) - Fapp(c, t, n) - Iast*eta(t));
      maxerr = fmax(err, maxerr);
    }
    end = clock();

    time = (double)(end - start) / CLOCKS_PER_SEC;
    printf("%d\t%e\t%e\n", n, err, time);

    FreeVec(J);
    FreeVec(g);
    FreeVec(c);
  }

  return EXIT_SUCCESS;
}

max_stars_repo_name: okayamat/sinc-indef

max_issues_repo_name: okayamat/sinc-indef

max_forks_repo_name: okayamat/sinc-indef


text:
#include "quac_p.h"
#include "quac.h"
#include "operators_p.h"
#include "operators.h"
#include <petsc.h>

int petsc_initialized = 0;
int nid;
int np;

/*
 * QuaC_initialize initializes petsc, gets each core's nid, and lets the
 * rest of the program know that it has been initialized.
 * Inputs:
 *       int argc, char **args - command line input, for PETSc
 */
void QuaC_initialize(int argc,char **args){

  /* Initialize Petsc */
  PetscInitialize(&argc,&args,(char*)0,NULL);
#if !defined(PETSC_USE_COMPLEX)
  SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
#endif
  /* Get core's id */
  MPI_Comm_rank(PETSC_COMM_WORLD,&nid);
  /* Get number of processors */
  MPI_Comm_size(PETSC_COMM_WORLD,&np);

  petsc_initialized = 1;
  PetscLogStageRegister("Pre-solve",&pre_solve_stage);
  PetscLogStageRegister("Solve",&solve_stage);
  PetscLogStageRegister("Post-solve",&post_solve_stage);
  /* Register events */
  PetscClassIdRegister("QuaC Class",&quac_class_id);
  PetscLogEventRegister("add_lin",quac_class_id,&add_lin_event);
  PetscLogEventRegister("add_to_ham",quac_class_id,&add_to_ham_event);
  PetscLogEventRegister("add_lin_recovery",quac_class_id,&add_lin_recovery_event);
  PetscLogEventRegister("_qc_event",quac_class_id,&_qc_event_function_event);
  PetscLogEventRegister("_qc_postevent",quac_class_id,&_qc_postevent_function_event);
  PetscLogEventRegister("_apply_gate",quac_class_id,&_apply_gate_event);

  PetscLogStagePush(pre_solve_stage);

}

/*
 * QuaC_clear clears the internal state of many of QuaC's
 * variables so that multiple systems can be run in one file.
 */

void QuaC_clear(){
  int i;
  /* Destroy Matrix */
  MatDestroy(&full_A);
  MatDestroy(&ham_A);
  MatDestroy(&full_stiff_A);
  MatDestroy(&ham_stiff_A);

  for (i=0;i<_num_time_dep;i++){
    MatDestroy(&_time_dep_list[i].mat);
  }
  //stab_added       = 0;
  _print_dense_ham = 0;
  _num_time_dep = 0;
  op_initialized = 0;
}

/*
 * QuaC_finalize finalizes petsc and destroys full_A.
 * The user is responsible for freeing all of the objects
 * using destroy_*
 */

void QuaC_finalize(){
  int i;
  /* Destroy Matrix */
  MatDestroy(&full_A);
  MatDestroy(&ham_A);
  MatDestroy(&full_stiff_A);
  MatDestroy(&ham_stiff_A);

  for (i=0;i<_num_time_dep;i++){
    MatDestroy(&_time_dep_list[i].mat);
  }
  /* Finalize Petsc */
  PetscLogStagePop();
  PetscFinalize();
  return;
}

/*
 * destroy_op frees the memory from op.
 * Inputs:
 *       operator *op - pointer to operator to be freed
 */

void destroy_op(operator *op){

  free((*op)->dag);
  free((*op)->n);
  free(*op);
}

/*
 * destroy_vec frees the memory from a vec.
 * Inputs:
 *       vec_op *op - pointer to vec_op to be freed
 */

void destroy_vec(vec_op *op){
  int num_levels,i;
  num_levels = (*op)[0]->my_levels;

  /* Free each up in the array */
  for (i=0;i<num_levels;i++){
    free((*op)[i]);
  }
  free(*op);
}

max_stars_repo_name: sgulania/QuaC

max_issues_repo_name: sgulania/QuaC

max_forks_repo_name: sgulania/QuaC


text:
/******************************************************************************
*
* IMD -- The ITAP Molecular Dynamics Program
*
* Copyright 1996-2007 Institute for Theoretical and Applied Physics,
* University of Stuttgart, D-70550 Stuttgart
*
******************************************************************************/

/******************************************************************************
*
* imd.h -- Header file for all modules of IMD
*
******************************************************************************/

/******************************************************************************
* $Revision$
* $Date$
******************************************************************************/

/* C stuff */
#ifdef MEMALIGN
#define  _XOPEN_SOURCE 600
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <time.h>
#include <unistd.h>
#include <errno.h>

//MYMOD
#include <stdbool.h>
#include "nn_interpol/nn.h"
#include "nn_interpol/delaunay.h"
//ENDOF MYMOD
//
//
#ifdef CBE
#define USE_WALLTIME
#endif

/* support for timers */
#ifndef MPI
#if defined(USE_RUSAGE) || defined(USE_WALLTIME)
#include <sys/time.h>
#include <sys/resource.h>
#else
#include <sys/times.h>
#include <sys/types.h>
#endif
#endif

/* Machine specific headers */
#if defined(MPI) || defined(NEB)
#include <mpi.h>
#ifdef MPE
#include <mpe.h>
#endif
#endif
#ifdef OMP
#include <omp.h>
#endif

/* FFT for diffraction patterns */
#ifdef DIFFPAT
#include <fftw3.h>
#endif

/* IMD version */
#include "version.h"

/* Configuration */
#include "config.h"

/* Data types */
#include "types.h"

/* Some constants */
#include "constants.h"

/* Some makros */
#include "makros.h"

/* Function Prototypes */
#include "prototypes.h"

/* Global Variables */
#include "globals.h"

/************
 * MY MOD   *
 ************/
#include <complex.h>

#ifdef COLRAD

// #include <gsl/gsl_sf_expint.h>       //FUER EXPONENTIAL INTEGRAL
#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */
#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */
#include <sunmatrix/sunmatrix_dense.h>
//#include <sundials/sundials_types.h> /* definition of type realtype */
#include "sunnonlinsol/sunnonlinsol_newton.h"
#define LAPACK  /* LAPACK SOLVER MULTO BENE*/
#ifdef LAPACK
#include <sunlinsol/sunlinsol_lapackdense.h>
#else
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
#endif
#endif //ifdef COLRAD

ext: h

max_stars_repo_name: fmqeisfeld/IMD

max_issues_repo_name: fmqeisfeld/IMD

max_forks_repo_name: fmqeisfeld/IMD


text:
/* block/gsl_block_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_BLOCK_UINT_H__
#define __GSL_BLOCK_UINT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_errno.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

struct gsl_block_uint_struct
{
  size_t size;
  unsigned int *data;
};

typedef struct gsl_block_uint_struct gsl_block_uint;

GSL_FUN gsl_block_uint *gsl_block_uint_alloc (const size_t n);
GSL_FUN gsl_block_uint *gsl_block_uint_calloc (const size_t n);
GSL_FUN void gsl_block_uint_free (gsl_block_uint * b);

GSL_FUN int gsl_block_uint_fread (FILE * stream, gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fwrite (FILE * stream, const gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fscanf (FILE * stream, gsl_block_uint * b);
GSL_FUN int gsl_block_uint_fprintf (FILE * stream, const gsl_block_uint * b, const char *format);

GSL_FUN int gsl_block_uint_raw_fread (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fwrite (FILE * stream, const unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fscanf (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
GSL_FUN int gsl_block_uint_raw_fprintf (FILE * stream, const unsigned int * b, const size_t n, const size_t stride, const char *format);

GSL_FUN size_t gsl_block_uint_size (const gsl_block_uint * b);
GSL_FUN unsigned int * gsl_block_uint_data (const gsl_block_uint * b);

__END_DECLS

#endif /* __GSL_BLOCK_UINT_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include <gsl/gsl>

namespace onnxruntime {
// Inspired by Fekir's Blog https://fekir.info/post/span-the-missing-constructor/
// Used under MIT license

// Use AsSpan for less typing on any container including initializer list to create a span
// (unnamed, untyped initializer list does not automatically convert to gsl::span).
// {1, 2, 3} as such does not have a type
// (see https://scottmeyers.blogspot.com/2014/03/if-braced-initializers-have-no-type-why.html)
//
//   Example: AsSpan({1, 2, 3}) results in gsl::span<const int>
//
// The above would deduce to std::initializer_list<int> and the result is gsl::span<const int>
//
// AsSpan<int64_t>({1, 2, 3}) produces gsl::span<const int64_t>
//
// We can also do std::array<int64_t, 3>{1, 2, 3} that can be automatically converted to span
// without memory allocation.
//
// If type conversion is not required, then for C++17 std::array template parameters are
// auto-deduced. Example: std::array{1, 2, 3}.
// We are aiming at not allocating memory dynamically.

namespace details {
template <class P>
constexpr auto AsSpanImpl(P* p, size_t s) {
  return gsl::span<P>(p, s);
}
}  // namespace details

template <class C>
constexpr auto AsSpan(C& c) {
  return details::AsSpanImpl(c.data(), c.size());
}

template <class C>
constexpr auto AsSpan(const C& c) {
  return details::AsSpanImpl(c.data(), c.size());
}

template <class C>
constexpr auto AsSpan(C&& c) {
  return details::AsSpanImpl(c.data(), c.size());
}

template <class T>
constexpr auto AsSpan(std::initializer_list<T> c) {
  return details::AsSpanImpl(c.begin(), c.size());
}

template <class T, size_t N>
constexpr auto AsSpan(T (&arr)[N]) {
  return details::AsSpanImpl(arr, N);
}

template <class T, size_t N>
constexpr auto AsSpan(const T (&arr)[N]) {
  return details::AsSpanImpl(arr, N);
}

}

max_stars_repo_name: SiriusKY/onnxruntime

max_issues_repo_name: SiriusKY/onnxruntime

max_forks_repo_name: SiriusKY/onnxruntime


text:
/* multilarge_nlinear/gsl_multilarge_nlinear.h
 *
 * Copyright (C) 2015, 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.
 */

#ifndef __GSL_MULTILARGE_NLINEAR_H__
#define __GSL_MULTILARGE_NLINEAR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_blas.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 enum
{
  GSL_MULTILARGE_NLINEAR_FWDIFF,
  GSL_MULTILARGE_NLINEAR_CTRDIFF
} gsl_multilarge_nlinear_fdtype;

/* Definition of vector-valued functions and gradient with parameters
   based on gsl_vector */

typedef struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  int (* df) (CBLAS_TRANSPOSE_t TransJ, const gsl_vector * x,
              const gsl_vector * u, void * params, gsl_vector * v,
              gsl_matrix * JTJ);
  int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
               gsl_vector * fvv);
  size_t n;        /* number of functions */
  size_t p;        /* number of independent variables */
  void * params;   /* user parameters */
  size_t nevalf;   /* number of function evaluations */
  size_t nevaldfu; /* number of Jacobian matrix-vector evaluations */
  size_t nevaldf2; /* number of Jacobian J^T J evaluations */
  size_t nevalfvv; /* number of fvv evaluations */
} gsl_multilarge_nlinear_fdf;

/* trust region subproblem method */
typedef struct
{
  const char *name;
  void * (*alloc) (const void * params, const size_t n, const size_t p);
  int (*init) (const void * vtrust_state, void * vstate);
  int (*preloop) (const void * vtrust_state, void * vstate);
  int (*step) (const void * vtrust_state, const double delta,
               gsl_vector * dx, void * vstate);
  int (*preduction) (const void * vtrust_state, const gsl_vector * dx,
                     double * pred, void * vstate);
  void (*free) (void * vstate);
} gsl_multilarge_nlinear_trs;

/* scaling matrix specification */
typedef struct
{
  const char *name;
  int (*init) (const gsl_matrix * JTJ, gsl_vector * diag);
  int (*update) (const gsl_matrix * JTJ, gsl_vector * diag);
} gsl_multilarge_nlinear_scale;

/*
 * linear least squares solvers - there are three steps to
 * solving a least squares problem using a direct method:
 *
 * 1. init: called once per iteration when a new Jacobian matrix
 *          is required; form normal equations matrix J^T J
 * 2. presolve: called each time a new LM parameter value mu is available;
 *              used for cholesky method in order to factor
 *              the (J^T J + mu D^T D) matrix
 * 3. solve: solve the least square system for a given rhs
 */
typedef struct
{
  const char *name;
  void * (*alloc) (const size_t n, const size_t p);
  int (*init) (const void * vtrust_state, void * vstate);
  int (*presolve) (const double mu, const void * vtrust_state, void * vstate);
  int (*solve) (const gsl_vector * g, gsl_vector * x,
                const void * vtrust_state, void * vstate);
  int (*rcond) (double * rcond, const gsl_matrix * JTJ, void * vstate);
  int (*covar) (const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate);
  void (*free) (void * vstate);
} gsl_multilarge_nlinear_solver;

/* tunable parameters */
typedef struct
{
  const gsl_multilarge_nlinear_trs *trs;       /* trust region subproblem method */
  const gsl_multilarge_nlinear_scale *scale;   /* scaling method */
  const gsl_multilarge_nlinear_solver *solver; /* solver method */
  gsl_multilarge_nlinear_fdtype fdtype;        /* finite difference method */
  double factor_up;                            /* factor for increasing trust radius */
  double factor_down;                          /* factor for decreasing trust radius */
  double avmax;                                /* max allowed |a|/|v| */
  double h_df;                                 /* step size for finite difference Jacobian */
  double h_fvv;                                /* step size for finite difference fvv */
  size_t max_iter;                             /* maximum iterations for trs method */
  double tol;                                  /* tolerance for solving trs */
} gsl_multilarge_nlinear_parameters;

typedef struct
{
  const char *name;
  void * (*alloc) (const gsl_multilarge_nlinear_parameters * params,
                   const size_t n, const size_t p);
  int (*init) (void * state, const gsl_vector * wts,
               gsl_multilarge_nlinear_fdf * fdf, const gsl_vector * x,
               gsl_vector * f, gsl_vector * g, gsl_matrix * JTJ);
  int (*iterate) (void * state, const gsl_vector * wts,
                  gsl_multilarge_nlinear_fdf * fdf, gsl_vector * x,
                  gsl_vector * f, gsl_vector * g, gsl_matrix * JTJ,
                  gsl_vector * dx);
  int (*rcond) (double * rcond, const gsl_matrix * JTJ, void * state);
  int (*covar) (const gsl_matrix * JTJ, gsl_matrix * covar, void * state);
  double (*avratio) (void * state);
  void (*free) (void * state);
} gsl_multilarge_nlinear_type;

/* current state passed to low-level trust region algorithms */
typedef struct
{
  const gsl_vector * x;             /* parameter values x */
  const gsl_vector * f;             /* residual vector f(x) */
  const gsl_vector * g;             /* gradient J^T f */
  const gsl_matrix * JTJ;           /* matrix J^T J */
  const gsl_vector * diag;          /* scaling matrix D */
  const gsl_vector * sqrt_wts;      /* sqrt(diag(W)) or NULL for unweighted */
  const double *mu;                 /* LM parameter */
  const gsl_multilarge_nlinear_parameters * params;
  void *solver_state;               /* workspace for direct least squares solver */
  gsl_multilarge_nlinear_fdf * fdf;
  double *avratio;                  /* |a| / |v| */
} gsl_multilarge_nlinear_trust_state;

typedef struct
{
  const gsl_multilarge_nlinear_type * type;
  gsl_multilarge_nlinear_fdf * fdf ;
  gsl_vector * x;             /* parameter values x */
  gsl_vector * f;             /* residual vector f(x) */
  gsl_vector * dx;            /* step dx */
  gsl_vector * g;             /* gradient J^T f */
  gsl_matrix * JTJ;           /* matrix J^T J */
  gsl_vector * sqrt_wts_work; /* sqrt(W) */
  gsl_vector * sqrt_wts;      /* ptr to sqrt_wts_work, or NULL if not using weights */
  size_t n;                   /* number of residuals */
  size_t p;                   /* number of parameters */
  size_t niter;               /* number of iterations performed */
  gsl_multilarge_nlinear_parameters params;
  void *state;
} gsl_multilarge_nlinear_workspace;

GSL_FUN gsl_multilarge_nlinear_workspace *
gsl_multilarge_nlinear_alloc (const gsl_multilarge_nlinear_type * T,
                              const gsl_multilarge_nlinear_parameters * params,
                              size_t n, size_t p);

GSL_FUN void gsl_multilarge_nlinear_free (gsl_multilarge_nlinear_workspace * w);

GSL_FUN gsl_multilarge_nlinear_parameters gsl_multilarge_nlinear_default_parameters(void);

GSL_FUN int
gsl_multilarge_nlinear_init (const gsl_vector * x,
                             gsl_multilarge_nlinear_fdf * fdf,
                             gsl_multilarge_nlinear_workspace * w);

GSL_FUN int gsl_multilarge_nlinear_winit (const gsl_vector * x,
                                  const gsl_vector * wts,
                                  gsl_multilarge_nlinear_fdf * fdf,
                                  gsl_multilarge_nlinear_workspace * w);

GSL_FUN int
gsl_multilarge_nlinear_iterate (gsl_multilarge_nlinear_workspace * w);

GSL_FUN double
gsl_multilarge_nlinear_avratio (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN int
gsl_multilarge_nlinear_rcond (double * rcond, const gsl_multilarge_nlinear_workspace * w);

GSL_FUN int
gsl_multilarge_nlinear_covar (gsl_matrix * covar, gsl_multilarge_nlinear_workspace * w);

GSL_FUN int
gsl_multilarge_nlinear_driver (const size_t maxiter,
                               const double xtol,
                               const double gtol,
                               const double ftol,
                               void (*callback)(const size_t iter, void *params,
                                                const gsl_multilarge_nlinear_workspace *w),
                               void *callback_params,
                               int *info,
                               gsl_multilarge_nlinear_workspace * w);

GSL_FUN const char *
gsl_multilarge_nlinear_name (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN gsl_vector *
gsl_multilarge_nlinear_position (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN gsl_vector *
gsl_multilarge_nlinear_residual (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN gsl_vector *
gsl_multilarge_nlinear_step (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN size_t
gsl_multilarge_nlinear_niter (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN const char *
gsl_multilarge_nlinear_trs_name (const gsl_multilarge_nlinear_workspace * w);

GSL_FUN int gsl_multilarge_nlinear_eval_f(gsl_multilarge_nlinear_fdf *fdf,
                                  const gsl_vector *x,
                                  const gsl_vector *swts,
                                  gsl_vector *y);

GSL_FUN int
gsl_multilarge_nlinear_eval_df(const CBLAS_TRANSPOSE_t TransJ,
                               const gsl_vector *x,
                               const gsl_vector *f,
                               const gsl_vector *u,
                               const gsl_vector *swts,
                               const double h,
                               const gsl_multilarge_nlinear_fdtype fdtype,
                               gsl_multilarge_nlinear_fdf *fdf,
                               gsl_vector *v,
                               gsl_matrix *JTJ,
                               gsl_vector *work);

GSL_FUN int
gsl_multilarge_nlinear_eval_fvv(const double h,
                                const gsl_vector *x,
                                const gsl_vector *v,
                                const gsl_vector *f,
                                const gsl_vector *swts,
                                gsl_multilarge_nlinear_fdf *fdf,
                                gsl_vector *yvv,
                                gsl_vector *work);

/* convergence.c */
GSL_FUN int
gsl_multilarge_nlinear_test (const double xtol, const double gtol,
                             const double ftol, int *info,
                             const gsl_multilarge_nlinear_workspace * w);

/* fdjac.c */
GSL_FUN int
gsl_multilarge_nlinear_df(const double h, const gsl_multilarge_nlinear_fdtype fdtype,
                          const gsl_vector *x, const gsl_vector *wts,
                          gsl_multilarge_nlinear_fdf *fdf,
                          const gsl_vector *f, gsl_matrix *J, gsl_vector *work);

/* fdfvv.c */
GSL_FUN int
gsl_multilarge_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
                             const gsl_vector *f, const gsl_matrix *J,
                             const gsl_vector *swts, gsl_multilarge_nlinear_fdf *fdf,
                             gsl_vector *fvv, gsl_vector *work);

/* top-level algorithms */
GSL_VAR const gsl_multilarge_nlinear_type * gsl_multilarge_nlinear_trust;

/* trust region subproblem methods */
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_lm;
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_lmaccel;
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_dogleg;
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_ddogleg;
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_subspace2D;
GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_cgst;

/* scaling matrix strategies */
GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_levenberg;
GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_marquardt;
GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_more;

/* linear solvers */
GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_cholesky;
GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_mcholesky;
GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_none;

__END_DECLS

#endif /* __GSL_MULTILARGE_NLINEAR_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* multiroots/hybrid.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>

#include "dogleg.c"

typedef struct
{
  size_t iter;
  size_t ncfail;
  size_t ncsuc;
  size_t nslow1;
  size_t nslow2;
  double fnorm;
  double delta;
  gsl_matrix *J;
  gsl_matrix *q;
  gsl_matrix *r;
  gsl_vector *tau;
  gsl_vector *diag;
  gsl_vector *qtf;
  gsl_vector *newton;
  gsl_vector *gradient;
  gsl_vector *x_trial;
  gsl_vector *f_trial;
  gsl_vector *df;
  gsl_vector *qtdf;
  gsl_vector *rdx;
  gsl_vector *w;
  gsl_vector *v;
}
hybrid_state_t;

static int hybrid_alloc (void *vstate, size_t n);
static int hybrid_set (void *vstate, gsl_multiroot_function * func,
                       gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static int hybrids_set (void *vstate, gsl_multiroot_function * func,
                        gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static int hybrid_set_impl (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
                gsl_vector * f, gsl_vector * dx, int scale);
static int hybrid_iterate (void *vstate, gsl_multiroot_function * func,
                           gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static void hybrid_free (void *vstate);
static int hybrid_iterate_impl (void *vstate, gsl_multiroot_function * func,
                    gsl_vector * x, gsl_vector * f, gsl_vector * dx,
                    int scale);

static int
hybrid_alloc (void *vstate, size_t n)
{
  hybrid_state_t *state = (hybrid_state_t *) vstate;
  gsl_matrix *J, *q, *r;
  gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial,
    *df, *qtdf, *rdx, *w, *v;

  J = gsl_matrix_calloc (n, n);

  if (J == 0)
    {
      GSL_ERROR ("failed to allocate space for J", GSL_ENOMEM);
    }

  state->J = J;

  q = gsl_matrix_calloc (n, n);

  if (q == 0)
    {
      gsl_matrix_free (J);

      GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM);
    }

  state->q = q;

  r = gsl_matrix_calloc (n, n);

  if (r == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);

      GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM);
    }

  state->r = r;

  tau = gsl_vector_calloc (n);

  if (tau == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);

      GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM);
    }

  state->tau = tau;

  diag = gsl_vector_calloc (n);

  if (diag == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);

      GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
    }

  state->diag = diag;

  qtf = gsl_vector_calloc (n);

  if (qtf == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);

      GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM);
    }

  state->qtf = qtf;

  newton = gsl_vector_calloc (n);

  if (newton == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);

      GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM);
    }

  state->newton = newton;

  gradient = gsl_vector_calloc (n);

  if (gradient == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);

      GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM);
    }

  state->gradient = gradient;

  x_trial = gsl_vector_calloc (n);

  if (x_trial == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);

      GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
    }

  state->x_trial = x_trial;

  f_trial = gsl_vector_calloc (n);

  if (f_trial == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);

      GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
    }

  state->f_trial = f_trial;

  df = gsl_vector_calloc (n);

  if (df == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);

      GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM);
    }

  state->df = df;

  qtdf = gsl_vector_calloc (n);

  if (qtdf == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);

      GSL_ERROR ("failed to allocate space for qtdf", GSL_ENOMEM);
    }

  state->qtdf = qtdf;

  rdx = gsl_vector_calloc (n);

  if (rdx == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (qtdf);

      GSL_ERROR ("failed to allocate space for rdx", GSL_ENOMEM);
    }

  state->rdx = rdx;

  w = gsl_vector_calloc (n);

  if (w == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (qtdf);
      gsl_vector_free (rdx);

      GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
    }

  state->w = w;

  v = gsl_vector_calloc (n);

  if (v == 0)
    {
      gsl_matrix_free (J);
      gsl_matrix_free (q);
      gsl_matrix_free (r);
      gsl_vector_free (tau);
      gsl_vector_free (diag);
      gsl_vector_free (qtf);
      gsl_vector_free (newton);
      gsl_vector_free (gradient);
      gsl_vector_free (x_trial);
      gsl_vector_free (f_trial);
      gsl_vector_free (df);
      gsl_vector_free (qtdf);
      gsl_vector_free (rdx);
      gsl_vector_free (w);

      GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM);
    }

  state->v = v;

  return GSL_SUCCESS;
}

static int
hybrid_set (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
            gsl_vector * f, gsl_vector * dx)
{
  int status = hybrid_set_impl (vstate, func, x, f, dx, 0);
  return status;
}

static int
hybrids_set (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
             gsl_vector * f, gsl_vector * dx)
{
  int status = hybrid_set_impl (vstate, func, x, f, dx, 1);
  return status;
}

static int
hybrid_set_impl (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
     gsl_vector * f, gsl_vector * dx, int scale)
{
  hybrid_state_t *state = (hybrid_state_t *) vstate;

  gsl_matrix *J = state->J;
  gsl_matrix *q = state->q;
  gsl_matrix *r = state->r;
  gsl_vector *tau = state->tau;
  gsl_vector *diag = state->diag;

  int status;

  status = GSL_MULTIROOT_FN_EVAL (func, x, f);

  if (status)
    {
      return status;
    }

  status = gsl_multiroot_fdjacobian (func, x, f, GSL_SQRT_DBL_EPSILON, J);

  if (status)
    {
      return status;
    }

  state->iter = 1;
  state->fnorm = enorm (f);
  state->ncfail = 0;
  state->ncsuc = 0;
  state->nslow1 = 0;
  state->nslow2 = 0;

  gsl_vector_set_all (dx, 0.0);

  /* Store column norms in diag */

  if (scale)
    compute_diag (J, diag);
  else
    gsl_vector_set_all (diag, 1.0);

  /* Set delta to factor |D x| or to factor if |D x| is zero */

  state->delta = compute_delta (diag, x);

  /* Factorize J into QR decomposition */

  status = gsl_linalg_QR_decomp (J, tau);

  if (status)
    {
      return status;
    }

  status = gsl_linalg_QR_unpack (J, tau, q, r);

  return status;
}

static int
hybrid_iterate (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
                gsl_vector * f, gsl_vector * dx)
{
  int status = hybrid_iterate_impl (vstate, func, x, f, dx, 0);
  return status;
}

static int
hybrids_iterate (void *vstate, gsl_multiroot_function * func, gsl_vector * x,
                 gsl_vector * f, gsl_vector * dx)
{
  int status = hybrid_iterate_impl (vstate, func, x, f, dx, 1);
  return status;
}

static int
hybrid_iterate_impl (void *vstate, gsl_multiroot_function * func,
                     gsl_vector * x,
                     gsl_vector * f, gsl_vector * dx, int scale)
{
  hybrid_state_t *state = (hybrid_state_t *) vstate;

  const double fnorm = state->fnorm;

  gsl_matrix *J = state->J;
  gsl_matrix *q = state->q;
  gsl_matrix *r = state->r;
  gsl_vector *tau = state->tau;
  gsl_vector *diag = state->diag;
  gsl_vector *qtf = state->qtf;
  gsl_vector *x_trial = state->x_trial;
  gsl_vector *f_trial = state->f_trial;
  gsl_vector *df = state->df;
  gsl_vector *qtdf = state->qtdf;
  gsl_vector *rdx = state->rdx;
  gsl_vector *w = state->w;
  gsl_vector *v = state->v;

  double prered, actred;
  double pnorm, fnorm1, fnorm1p;
  double ratio;
  double p1 = 0.1, p5 = 0.5, p001 = 0.001, p0001 = 0.0001;

  /* Compute qtf = Q^T f */

  compute_qtf (q, f, qtf);

  /* Compute dogleg step */

  dogleg (r, qtf, diag, state->delta, state->newton, state->gradient, dx);

  /* Take a trial step */

  compute_trial_step (x, dx, state->x_trial);

  pnorm = scaled_enorm (diag, dx);

  if (state->iter == 1)
    {
      if (pnorm < state->delta)
        {
          state->delta = pnorm;
        }
    }

  /* Evaluate function at x + p */

  {
    int status = GSL_MULTIROOT_FN_EVAL (func, x_trial, f_trial);

    if (status != GSL_SUCCESS)
      {
        return GSL_EBADFUNC;
      }
  }

  /* Set df = f_trial - f */

  compute_df (f_trial, f, df);

  /* Compute the scaled actual reduction */

  fnorm1 = enorm (f_trial);

  actred = compute_actual_reduction (fnorm, fnorm1);

  /* Compute rdx = R dx */

  compute_rdx (r, dx, rdx);

  /* Compute the scaled predicted reduction phi1p = |Q^T f + R dx| */

  fnorm1p = enorm_sum (qtf, rdx);

  prered = compute_predicted_reduction (fnorm, fnorm1p);

  /* Compute the ratio of the actual to predicted reduction */

  if (prered > 0)
    {
      ratio = actred / prered;
    }
  else
    {
      ratio = 0;
    }

  /* Update the step bound */

  if (ratio < p1)
    {
      state->ncsuc = 0;
      state->ncfail++;
      state->delta *= p5;
    }
  else
    {
      state->ncfail = 0;
      state->ncsuc++;

      if (ratio >= p5 || state->ncsuc > 1)
        state->delta = GSL_MAX (state->delta, pnorm / p5);
      if (fabs (ratio - 1) <= p1)
        state->delta = pnorm / p5;
    }

  /* Test for successful iteration */

  if (ratio >= p0001)
    {
      gsl_vector_memcpy (x, x_trial);
      gsl_vector_memcpy (f, f_trial);
      state->fnorm = fnorm1;
      state->iter++;
    }

  /* Determine the progress of the iteration */

  state->nslow1++;
  if (actred >= p001)
    state->nslow1 = 0;

  if (actred >= p1)
    state->nslow2 = 0;

  if (state->ncfail == 2)
    {
      gsl_multiroot_fdjacobian (func, x, f, GSL_SQRT_DBL_EPSILON, J);

      state->nslow2++;

      if (state->iter == 1)
        {
          if (scale)
            compute_diag (J, diag);
          state->delta = compute_delta (diag, x);
        }
      else
        {
          if (scale)
            update_diag (J, diag);
        }

      /* Factorize J into QR decomposition */

      gsl_linalg_QR_decomp (J, tau);
      gsl_linalg_QR_unpack (J, tau, q, r);

      return GSL_SUCCESS;
    }

  /* Compute qtdf = Q^T df, w = (Q^T df - R dx)/|dx|,  v = D^2 dx/|dx| */

  compute_qtf (q, df, qtdf);

  compute_wv (qtdf, rdx, dx, diag, pnorm, w, v);

  /* Rank-1 update of the jacobian Q'R' = Q(R + w v^T) */

  gsl_linalg_QR_update (q, r, w, v);

  /* No progress as measured by jacobian evaluations */

  if (state->nslow2 == 5)
    {
      return GSL_ENOPROGJ;
    }

  /* No progress as measured by function evaluations */

  if (state->nslow1 == 10)
    {
      return GSL_ENOPROG;
    }

  return GSL_SUCCESS;
}

static void
hybrid_free (void *vstate)
{
  hybrid_state_t *state = (hybrid_state_t *) vstate;

  gsl_vector_free (state->v);
  gsl_vector_free (state->w);
  gsl_vector_free (state->rdx);
  gsl_vector_free (state->qtdf);
  gsl_vector_free (state->df);
  gsl_vector_free (state->f_trial);
  gsl_vector_free (state->x_trial);
  gsl_vector_free (state->gradient);
  gsl_vector_free (state->newton);
  gsl_vector_free (state->qtf);
  gsl_vector_free (state->diag);
  gsl_vector_free (state->tau);
  gsl_matrix_free (state->r);
  gsl_matrix_free (state->q);
  gsl_matrix_free (state->J);
}

static const gsl_multiroot_fsolver_type hybrid_type = {
  "hybrid",                     /* name */
  sizeof (hybrid_state_t),
  &hybrid_alloc,
  &hybrid_set,
  &hybrid_iterate,
  &hybrid_free
};

static const gsl_multiroot_fsolver_type hybrids_type = {
  "hybrids",                    /* name */
  sizeof (hybrid_state_t),
  &hybrid_alloc,
  &hybrids_set,
  &hybrids_iterate,
  &hybrid_free
};

const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrid = &hybrid_type;
const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrids =
  &hybrids_type;

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet