text:
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "extra_compile_macros.h"

#if USE_GSL == 1
#include <gsl/gsl_sf_gamma.h>
#endif

double nan_density(double t, double *pars, double *q, int n_dim) { return NAN; }
double nan_value(double t, double *pars, double *q, int n_dim) { return NAN; }
void nan_gradient(double t, double *pars, double *q, int n_dim, double *grad) {}
void nan_hessian(double t, double *pars, double *q, int n_dim, double *hess) {}

double null_density(double t, double *pars, double *q, int n_dim) { return 0; }
double null_value(double t, double *pars, double *q, int n_dim) { return 0; }
void null_gradient(double t, double *pars, double *q, int n_dim, double *grad){}
void null_hessian(double t, double *pars, double *q, int n_dim, double *hess) {}

/* Note: many Hessians generated with sympy in
    gala-notebooks/Make-all-Hessians.ipynb
*/

/* ---------------------------------------------------------------------------
    Henon-Heiles potential
*/
double henon_heiles_value(double t, double *pars, double *q, int n_dim) {
    /*  no parameters... */
    return 0.5 * (q[0]*q[0] + q[1]*q[1] + 2*q[0]*q[0]*q[1] - 2/3.*q[1]*q[1]*q[1]);
}

void henon_heiles_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  no parameters... */
    grad[0] = grad[0] + q[0] + 2*q[0]*q[1];
    grad[1] = grad[1] + q[1] + q[0]*q[0] - q[1]*q[1];
}

void henon_heiles_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  no parameters... */
    double x = q[0];
    double y = q[1];

    double tmp_0 = 2.0 * y;
    double tmp_1 = 2.0 * x;

    hess[0] = hess[0] + tmp_0 + 1.0;
    hess[1] = hess[1] + tmp_1;
    hess[2] = hess[2] + tmp_1;
    hess[3] = hess[3] + 1.0 - tmp_0;
}

/* ---------------------------------------------------------------------------
    Kepler potential
*/
double kepler_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
    */
    double R;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    return -pars[0] * pars[1] / R;
}

void kepler_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
    */
    double R, fac;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    fac = pars[0] * pars[1] / (R*R*R);

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2];
}

double kepler_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
    */
    double r2;
    r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];

    if (r2 == 0.) {
        return INFINITY;
    } else {
        return 0.;
    }
}

void kepler_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
    */
    double G = pars[0];
    double m = pars[1];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = G*m;
    double tmp_5 = tmp_4/pow(tmp_3, 3.0/2.0);
    double tmp_6 = 3*tmp_4/pow(tmp_3, 5.0/2.0);
    double tmp_7 = tmp_6*x;
    double tmp_8 = -tmp_7*y;
    double tmp_9 = -tmp_7*z;
    double tmp_10 = -tmp_6*y*z;

    hess[0] = hess[0] + -tmp_0*tmp_6 + tmp_5;
    hess[1] = hess[1] + tmp_8;
    hess[2] = hess[2] + tmp_9;
    hess[3] = hess[3] + tmp_8;
    hess[4] = hess[4] + -tmp_1*tmp_6 + tmp_5;
    hess[5] = hess[5] + tmp_10;
    hess[6] = hess[6] + tmp_9;
    hess[7] = hess[7] + tmp_10;
    hess[8] = hess[8] + -tmp_2*tmp_6 + tmp_5;
}

/* ---------------------------------------------------------------------------
    Isochrone potential
*/
double isochrone_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (core scale)
    */
    double R2;
    R2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
    return -pars[0] * pars[1] / (sqrt(R2 + pars[2]*pars[2]) + pars[2]);
}

void isochrone_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (core scale)
    */
    double sqrt_r2_b2, fac, denom;
    sqrt_r2_b2 = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*pars[2]);
    denom = sqrt_r2_b2 * (sqrt_r2_b2 + pars[2])*(sqrt_r2_b2 + pars[2]);
    fac = pars[0] * pars[1] / denom;

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2];
}

double isochrone_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (core scale)
    */
    double r2, a, b;
    b = pars[2];
    r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
    a = sqrt(b*b + r2);

    return pars[1] * (3*(b+a)*a*a - r2*(b+3*a)) / (4*M_PI*pow(b+a,3)*a*a*a);
}

void isochrone_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (length scale)
    */
    double G = pars[0];
    double m = pars[1];
    double b = pars[2];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = pow(b, 2) + tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = sqrt(tmp_3);
    double tmp_5 = b + tmp_4;
    double tmp_6 = G*m;
    double tmp_7 = tmp_6/pow(tmp_5, 2);
    double tmp_8 = tmp_7/tmp_4;
    double tmp_9 = 2*tmp_6/(tmp_3*pow(tmp_5, 3));
    double tmp_10 = tmp_7/pow(tmp_3, 3.0/2.0);
    double tmp_11 = tmp_9*x;
    double tmp_12 = tmp_10*x;
    double tmp_13 = -tmp_11*y - tmp_12*y;
    double tmp_14 = -tmp_11*z - tmp_12*z;
    double tmp_15 = y*z;
    double tmp_16 = -tmp_10*tmp_15 - tmp_15*tmp_9;

    hess[0] = hess[0] + -tmp_0*tmp_10 - tmp_0*tmp_9 + tmp_8;
    hess[1] = hess[1] + tmp_13;
    hess[2] = hess[2] + tmp_14;
    hess[3] = hess[3] + tmp_13;
    hess[4] = hess[4] + -tmp_1*tmp_10 - tmp_1*tmp_9 + tmp_8;
    hess[5] = hess[5] + tmp_16;
    hess[6] = hess[6] + tmp_14;
    hess[7] = hess[7] + tmp_16;
    hess[8] = hess[8] + -tmp_10*tmp_2 - tmp_2*tmp_9 + tmp_8;

}

/* ---------------------------------------------------------------------------
    Hernquist sphere
*/
double hernquist_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double R;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    return -pars[0] * pars[1] / (R + pars[2]);
}

void hernquist_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double R, fac;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    fac = pars[0] * pars[1] / ((R + pars[2]) * (R + pars[2]) * R);

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2];
}

double hernquist_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double r, rho0;
    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    rho0 = pars[1]/(2*M_PI*pars[2]*pars[2]*pars[2]);
    return rho0 / ((r/pars[2]) * pow(1+r/pars[2],3));
}

void hernquist_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double G = pars[0];
    double m = pars[1];
    double c = pars[2];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = sqrt(tmp_3);
    double tmp_5 = c + tmp_4;
    double tmp_6 = G*m;
    double tmp_7 = tmp_6/pow(tmp_5, 2);
    double tmp_8 = tmp_7/tmp_4;
    double tmp_9 = 2*tmp_6/(tmp_3*pow(tmp_5, 3));
    double tmp_10 = tmp_7/pow(tmp_3, 3.0/2.0);
    double tmp_11 = tmp_9*x;
    double tmp_12 = tmp_10*x;
    double tmp_13 = -tmp_11*y - tmp_12*y;
    double tmp_14 = -tmp_11*z - tmp_12*z;
    double tmp_15 = y*z;
    double tmp_16 = -tmp_10*tmp_15 - tmp_15*tmp_9;

    hess[0] = hess[0] + -tmp_0*tmp_10 - tmp_0*tmp_9 + tmp_8;
    hess[1] = hess[1] + tmp_13;
    hess[2] = hess[2] + tmp_14;
    hess[3] = hess[3] + tmp_13;
    hess[4] = hess[4] + -tmp_1*tmp_10 - tmp_1*tmp_9 + tmp_8;
    hess[5] = hess[5] + tmp_16;
    hess[6] = hess[6] + tmp_14;
    hess[7] = hess[7] + tmp_16;
    hess[8] = hess[8] + -tmp_10*tmp_2 - tmp_2*tmp_9 + tmp_8;
}

/* ---------------------------------------------------------------------------
    Plummer sphere
*/
double plummer_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (length scale)
    */
    double R2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
    return -pars[0]*pars[1] / sqrt(R2 + pars[2]*pars[2]);
}

void plummer_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (length scale)
    */
    double R2b, fac;
    R2b = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*pars[2];
    fac = pars[0] * pars[1] / sqrt(R2b) / R2b;

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2];
}

double plummer_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (length scale)
    */
    double r2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
    return 3*pars[1] / (4*M_PI*pars[2]*pars[2]*pars[2]) * pow(1 + r2/(pars[2]*pars[2]), -2.5);
}

void plummer_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - b (length scale)
    */
    double G = pars[0];
    double m = pars[1];
    double b = pars[2];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = pow(b, 2) + tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = G*m;
    double tmp_5 = tmp_4/pow(tmp_3, 3.0/2.0);
    double tmp_6 = 3*tmp_4/pow(tmp_3, 5.0/2.0);
    double tmp_7 = tmp_6*x;
    double tmp_8 = -tmp_7*y;
    double tmp_9 = -tmp_7*z;
    double tmp_10 = -tmp_6*y*z;

    hess[0] = hess[0] + -tmp_0*tmp_6 + tmp_5;
    hess[1] = hess[1] + tmp_8;
    hess[2] = hess[2] + tmp_9;
    hess[3] = hess[3] + tmp_8;
    hess[4] = hess[4] + -tmp_1*tmp_6 + tmp_5;
    hess[5] = hess[5] + tmp_10;
    hess[6] = hess[6] + tmp_9;
    hess[7] = hess[7] + tmp_10;
    hess[8] = hess[8] + -tmp_2*tmp_6 + tmp_5;
}

/* ---------------------------------------------------------------------------
    Jaffe sphere
*/
double jaffe_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double R;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    return -pars[0] * pars[1] / pars[2] * log(1 + pars[2]/R);
}

void jaffe_gradient(double t, double *pars, double *q, int n_dim, double *grad){
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double R, fac;
    R = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    fac = pars[0] * pars[1] / pars[2] * (pars[2] / (R * (pars[2] + R)));

    grad[0] = grad[0] + fac*q[0]/R;
    grad[1] = grad[1] + fac*q[1]/R;
    grad[2] = grad[2] + fac*q[2]/R;
}

double jaffe_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double r, rho0;
    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    rho0 = pars[1] / (4*M_PI*pars[2]*pars[2]*pars[2]);
    return rho0 / (pow(r/pars[2],2) * pow(1+r/pars[2],2));
}

void jaffe_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - c (length scale)
    */
    double G = pars[0];
    double m = pars[1];
    double c = pars[2];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = 1.0/tmp_3;
    double tmp_5 = sqrt(tmp_3);
    double tmp_6 = c + tmp_5;
    double tmp_7 = pow(tmp_6, -2);
    double tmp_8 = tmp_7*x;
    double tmp_9 = 1.0/tmp_5;
    double tmp_10 = 1.0/tmp_6;
    double tmp_11 = tmp_10*tmp_9;
    double tmp_12 = G*m/c;
    double tmp_13 = tmp_12*(tmp_11*x - tmp_8);
    double tmp_14 = tmp_13*tmp_4;
    double tmp_15 = pow(tmp_3, -3.0/2.0);
    double tmp_16 = tmp_13*tmp_15*tmp_6;
    double tmp_17 = tmp_10*tmp_15;
    double tmp_18 = tmp_4*tmp_7;
    double tmp_19 = 2*tmp_9/pow(tmp_6, 3);
    double tmp_20 = tmp_11 - tmp_7;
    double tmp_21 = tmp_12*tmp_6;
    double tmp_22 = tmp_21*tmp_9;
    double tmp_23 = tmp_19*x;
    double tmp_24 = tmp_4*tmp_8;
    double tmp_25 = tmp_17*x;
    double tmp_26 = tmp_14*y - tmp_16*y + tmp_22*(tmp_23*y - tmp_24*y - tmp_25*y);
    double tmp_27 = tmp_4*z;
    double tmp_28 = tmp_13*tmp_27 - tmp_16*z + tmp_22*(tmp_23*z - tmp_24*z - tmp_25*z);
    double tmp_29 = tmp_7*y;
    double tmp_30 = tmp_11*y - tmp_29;
    double tmp_31 = tmp_12*tmp_30;
    double tmp_32 = tmp_15*tmp_21;
    double tmp_33 = tmp_30*tmp_32;
    double tmp_34 = y*z;
    double tmp_35 = tmp_22*(-tmp_17*tmp_34 + tmp_19*tmp_34 - tmp_27*tmp_29) + tmp_27*tmp_31 - tmp_33*z;
    double tmp_36 = tmp_11*z - tmp_7*z;

    hess[0] = hess[0] + tmp_14*x - tmp_16*x + tmp_22*(-tmp_0*tmp_17 - tmp_0*tmp_18 + tmp_0*tmp_19 + tmp_20);
    hess[1] = hess[1] + tmp_26;
    hess[2] = hess[2] + tmp_28;
    hess[3] = hess[3] + tmp_26;
    hess[4] = hess[4] + tmp_22*(-tmp_1*tmp_17 - tmp_1*tmp_18 + tmp_1*tmp_19 + tmp_20) + tmp_31*tmp_4*y - tmp_33*y;
    hess[5] = hess[5] + tmp_35;
    hess[6] = hess[6] + tmp_28;
    hess[7] = hess[7] + tmp_35;
    hess[8] = hess[8] + tmp_12*tmp_27*tmp_36 + tmp_22*(-tmp_17*tmp_2 - tmp_18*tmp_2 + tmp_19*tmp_2 + tmp_20) - tmp_32*tmp_36*z;
}

/* ---------------------------------------------------------------------------
    Power-law potential with exponential cutoff
*/
#if USE_GSL == 1

double safe_gamma_inc(double a, double x) {
    int N, m, n;
    double A = 1.;
    double B = 0.;
    double tmp;

    if (a > 0) {
        return gsl_sf_gamma_inc_P(a, x) * gsl_sf_gamma(a);;
    } else {
        N = (int) ceil(-a);

        for (n=0; n < N; n++) {
            A = A * (a + n);

            tmp = 1.;
            for (m=N-1; m > n; m--) {
                tmp = tmp * (a + m);
            }
            B = B + pow(x, a+n) * exp(-x) * tmp;
        }
        return (B + gsl_sf_gamma_inc_P(a + N, x) * gsl_sf_gamma(a + N)) / A;
    }
}

double powerlawcutoff_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            0 - G (Gravitational constant)
            1 - m (total mass)
            2 - a (power-law index)
            3 - c (cutoff radius)
    */
    double G = pars[0];
    double m = pars[1];
    double alpha = pars[2];
    double r_c = pars[3];
    double x = q[0];
    double y = q[1];
    double z = q[2];
    double r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);

    if (r == 0.) {
        return -INFINITY;
    } else {
        double tmp_0 = (1.0/2.0)*alpha;
        double tmp_1 = -tmp_0;
        double tmp_2 = tmp_1 + 1.5;
        double tmp_3 = pow(x, 2) + pow(y, 2) + pow(z, 2);
        double tmp_4 = tmp_3/pow(r_c, 2);
        double tmp_5 = G*m;
        double tmp_6 = tmp_5*safe_gamma_inc(tmp_2, tmp_4)/(sqrt(tmp_3)*tgamma(tmp_1 + 2.5));
        return tmp_0*tmp_6 - 3.0/2.0*tmp_6 + tmp_5*safe_gamma_inc(tmp_1 + 1, tmp_4)/(r_c*tgamma(tmp_2));
    }
}

double powerlawcutoff_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            0 - G (Gravitational constant)
            1 - m (total mass)
            2 - a (power-law index)
            3 - c (cutoff radius)
    */
    double r, A;
    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    A = pars[1] / (2*M_PI) * pow(pars[3], pars[2] - 3) / gsl_sf_gamma(0.5 * (3 - pars[2]));
    return A * pow(r, -pars[2]) * exp(-r*r / (pars[3]*pars[3]));
}

void powerlawcutoff_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            0 - G (Gravitational constant)
            1 - m (total mass)
            2 - a (power-law index)
            3 - c (cutoff radius)
    */
    double r, dPhi_dr;
    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    dPhi_dr = (pars[0] * pars[1] / (r*r) *
        gsl_sf_gamma_inc_P(0.5 * (3-pars[2]), r*r/(pars[3]*pars[3]))); // / gsl_sf_gamma(0.5 * (3-pars[2])));

    grad[0] = grad[0] + dPhi_dr * q[0]/r;
    grad[1] = grad[1] + dPhi_dr * q[1]/r;
    grad[2] = grad[2] + dPhi_dr * q[2]/r;
}

void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - alpha (exponent)
            - r_c (cutoff radius)
    */
    double G = pars[0];
    double m = pars[1];
    double alpha = pars[2];
    double r_c = pars[3];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = (1.0/2.0)*alpha;
    double tmp_5 = -tmp_4;
    double tmp_6 = tmp_5 + 1.5;
    double tmp_7 = pow(r_c, -2);
    double tmp_8 = tmp_3*tmp_7;
    double tmp_9 = G*m;
    double tmp_10 = tmp_9/tgamma(tmp_5 + 2.5);
    double tmp_11 = tmp_10*safe_gamma_inc(tmp_6, tmp_8);
    double tmp_12 = tmp_11/pow(tmp_3, 5.0/2.0);
    double tmp_13 = (9.0/2.0)*tmp_12;
    double tmp_14 = exp(-tmp_8);
    double tmp_15 = tmp_0*tmp_14;
    double tmp_16 = pow(tmp_8, -tmp_4)*tmp_9/tgamma(tmp_6);
    double tmp_17 = 4*tmp_16/pow(r_c, 5);
    double tmp_18 = alpha*tmp_0;
    double tmp_19 = (3.0/2.0)*tmp_12;
    double tmp_20 = 6*tmp_15;
    double tmp_21 = pow(r_c, -4);
    double tmp_22 = tmp_5 + 0.5;
    double tmp_23 = pow(tmp_8, tmp_22);
    double tmp_24 = tmp_10*tmp_23/sqrt(tmp_3);
    double tmp_25 = tmp_21*tmp_24;
    double tmp_26 = pow(tmp_3, -3.0/2.0);
    double tmp_27 = tmp_10*tmp_23*tmp_26*tmp_7;
    double tmp_28 = tmp_20*tmp_27;
    double tmp_29 = 2*tmp_14;
    double tmp_30 = tmp_18*tmp_29;
    double tmp_31 = tmp_16*tmp_29/pow(r_c, 3);
    double tmp_32 = tmp_31/tmp_3;
    double tmp_33 = tmp_27*tmp_30;
    double tmp_34 = tmp_11*tmp_26;
    double tmp_35 = tmp_14*tmp_24;
    double tmp_36 = tmp_35*tmp_7;
    double tmp_37 = alpha*tmp_36 + tmp_31 - tmp_34*tmp_4 + (3.0/2.0)*tmp_34 - 3*tmp_36;
    double tmp_38 = tmp_13*x;
    double tmp_39 = alpha*tmp_19;
    double tmp_40 = x*y;
    double tmp_41 = tmp_14*tmp_17;
    double tmp_42 = alpha*tmp_40;
    double tmp_43 = tmp_21*tmp_35;
    double tmp_44 = 6*tmp_40;
    double tmp_45 = tmp_14*tmp_27;
    double tmp_46 = tmp_44*tmp_45;
    double tmp_47 = tmp_29*tmp_42;
    double tmp_48 = tmp_27*tmp_47;
    double tmp_49 = -tmp_22*tmp_46 + tmp_22*tmp_48 - tmp_25*tmp_47 - tmp_32*tmp_42 - tmp_38*y + tmp_39*tmp_40 - tmp_40*tmp_41 + tmp_43*tmp_44 + tmp_46 - tmp_48;
    double tmp_50 = x*z;
    double tmp_51 = alpha*tmp_50;
    double tmp_52 = 6*tmp_50;
    double tmp_53 = tmp_45*tmp_52;
    double tmp_54 = tmp_29*tmp_51;
    double tmp_55 = tmp_27*tmp_54;
    double tmp_56 = -tmp_22*tmp_53 + tmp_22*tmp_55 - tmp_25*tmp_54 - tmp_32*tmp_51 - tmp_38*z + tmp_39*tmp_50 - tmp_41*tmp_50 + tmp_43*tmp_52 + tmp_53 - tmp_55;
    double tmp_57 = 6*tmp_1;
    double tmp_58 = tmp_45*tmp_57;
    double tmp_59 = alpha*tmp_1;
    double tmp_60 = tmp_29*tmp_59;
    double tmp_61 = tmp_27*tmp_60;
    double tmp_62 = y*z;
    double tmp_63 = alpha*tmp_62;
    double tmp_64 = 6*tmp_62;
    double tmp_65 = tmp_45*tmp_64;
    double tmp_66 = tmp_29*tmp_63;
    double tmp_67 = tmp_27*tmp_66;
    double tmp_68 = -tmp_13*tmp_62 - tmp_22*tmp_65 + tmp_22*tmp_67 - tmp_25*tmp_66 - tmp_32*tmp_63 + tmp_39*tmp_62 - tmp_41*tmp_62 + tmp_43*tmp_64 + tmp_65 - tmp_67;
    double tmp_69 = 6*tmp_2;
    double tmp_70 = tmp_45*tmp_69;
    double tmp_71 = alpha*tmp_2;
    double tmp_72 = tmp_29*tmp_71;
    double tmp_73 = tmp_27*tmp_72;

    hess[0] = hess[0] + -tmp_0*tmp_13 - tmp_15*tmp_17 + tmp_18*tmp_19 - tmp_18*tmp_32 + tmp_20*tmp_25 - tmp_22*tmp_28 + tmp_22*tmp_33 - tmp_25*tmp_30 + tmp_28 - tmp_33 + tmp_37;
    hess[1] = hess[1] + tmp_49;
    hess[2] = hess[2] + tmp_56;
    hess[3] = hess[3] + tmp_49;
    hess[4] = hess[4] + -tmp_1*tmp_13 + tmp_1*tmp_39 - tmp_1*tmp_41 - tmp_22*tmp_58 + tmp_22*tmp_61 - tmp_25*tmp_60 - tmp_32*tmp_59 + tmp_37 + tmp_43*tmp_57 + tmp_58 - tmp_61;
    hess[5] = hess[5] + tmp_68;
    hess[6] = hess[6] + tmp_56;
    hess[7] = hess[7] + tmp_68;
    hess[8] = hess[8] + -tmp_13*tmp_2 + tmp_2*tmp_39 - tmp_2*tmp_41 - tmp_22*tmp_70 + tmp_22*tmp_73 - tmp_25*tmp_72 - tmp_32*tmp_71 + tmp_37 + tmp_43*tmp_69 + tmp_70 - tmp_73;
}

#endif

/* ---------------------------------------------------------------------------
    Stone-Ostriker potential from Stone & Ostriker (2015)
*/
double stone_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - M (total mass)
            - r_c (core radius)
            - r_h (halo radius)
    */
    double r, u_c, u_h, fac;

    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    u_c = r / pars[2];
    u_h = r / pars[3];

    fac = 2*pars[0]*pars[1] / M_PI / (pars[3] - pars[2]);

    return -fac * (atan(u_h)/u_h - atan(u_c)/u_c +
                   0.5*log((r*r + pars[3]*pars[3])/(r*r + pars[2]*pars[2])));
}

void stone_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - M (total mass)
            - r_c (core radius)
            - r_h (halo radius)
    */
    double r, u_c, u_h, fac, dphi_dr;

    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    u_c = r / pars[2];
    u_h = r / pars[3];

    fac = 2*pars[0]*pars[1] / (M_PI*r*r) / (pars[2] - pars[3]);  // order flipped from value
    dphi_dr = fac * (pars[2]*atan(u_c) - pars[3]*atan(u_h));

    grad[0] = grad[0] + dphi_dr*q[0]/r;
    grad[1] = grad[1] + dphi_dr*q[1]/r;
    grad[2] = grad[2] + dphi_dr*q[2]/r;
}

double stone_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - M (total mass)
            - r_c (core radius)
            - r_h (halo radius)
    */
    double r, u_c, u_t, rho;
    rho = pars[1] * (pars[2] + pars[3]) / (2*M_PI*M_PI*pars[2]*pars[2]*pars[3]*pars[3]);

    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    u_c = r / pars[2];
    u_t = r / pars[3];

    return rho / ((1 + u_c*u_c)*(1 + u_t*u_t));
}

void stone_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_c (core radius)
            - r_h (halo radius)
    */
    double G = pars[0];
    double m = pars[1];
    double r_c = pars[2];
    double r_h = pars[3];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(r_h, 2);
    double tmp_1 = 1.0/tmp_0;
    double tmp_2 = pow(x, 2);
    double tmp_3 = pow(y, 2);
    double tmp_4 = pow(z, 2);
    double tmp_5 = tmp_2 + tmp_3 + tmp_4;
    double tmp_6 = tmp_1*tmp_5 + 1;
    double tmp_7 = 1.0/tmp_6;
    double tmp_8 = 3/pow(tmp_5, 2);
    double tmp_9 = tmp_2*tmp_8;
    double tmp_10 = pow(r_c, 2);
    double tmp_11 = 1.0/tmp_10;
    double tmp_12 = tmp_11*tmp_5 + 1;
    double tmp_13 = 1.0/tmp_12;
    double tmp_14 = tmp_10 + tmp_5;
    double tmp_15 = pow(tmp_14, -2);
    double tmp_16 = 8*tmp_15;
    double tmp_17 = tmp_0 + tmp_5;
    double tmp_18 = 8*tmp_17/pow(tmp_14, 3);
    double tmp_19 = 2/tmp_14;
    double tmp_20 = 2*tmp_15*tmp_17;
    double tmp_21 = tmp_19 - tmp_20;
    double tmp_22 = 1.0/tmp_17;
    double tmp_23 = 0.5*tmp_14*tmp_22;
    double tmp_24 = tmp_19*x - tmp_20*x;
    double tmp_25 = 1.0*tmp_22;
    double tmp_26 = tmp_24*tmp_25;
    double tmp_27 = sqrt(tmp_5);
    double tmp_28 = r_c*atan(tmp_27/r_c);
    double tmp_29 = 3/pow(tmp_5, 5.0/2.0);
    double tmp_30 = tmp_2*tmp_29;
    double tmp_31 = 1.0/tmp_5;
    double tmp_32 = 2*tmp_31;
    double tmp_33 = tmp_2*tmp_32;
    double tmp_34 = tmp_1/pow(tmp_6, 2);
    double tmp_35 = tmp_11/pow(tmp_12, 2);
    double tmp_36 = r_h*atan(tmp_27/r_h);
    double tmp_37 = 1.0*tmp_14/pow(tmp_17, 2);
    double tmp_38 = tmp_24*tmp_37;
    double tmp_39 = pow(tmp_5, -3.0/2.0);
    double tmp_40 = -tmp_13*tmp_31 + tmp_28*tmp_39 + tmp_31*tmp_7 - tmp_36*tmp_39;
    double tmp_41 = 2*G*m/(-3.1415926535897931*r_c + 3.1415926535897931*r_h);
    double tmp_42 = x*y;
    double tmp_43 = tmp_42*tmp_8;
    double tmp_44 = tmp_29*tmp_42;
    double tmp_45 = tmp_32*tmp_42;
    double tmp_46 = tmp_16*x;
    double tmp_47 = -tmp_41*(tmp_13*tmp_43 + tmp_23*(tmp_18*tmp_42 - tmp_46*y) + tmp_26*y - tmp_28*tmp_44 - tmp_34*tmp_45 + tmp_35*tmp_45 + tmp_36*tmp_44 - tmp_38*y - tmp_43*tmp_7);
    double tmp_48 = x*z;
    double tmp_49 = tmp_48*tmp_8;
    double tmp_50 = tmp_29*tmp_48;
    double tmp_51 = tmp_32*tmp_48;
    double tmp_52 = -tmp_41*(tmp_13*tmp_49 + tmp_23*(tmp_18*tmp_48 - tmp_46*z) + tmp_26*z - tmp_28*tmp_50 - tmp_34*tmp_51 + tmp_35*tmp_51 + tmp_36*tmp_50 - tmp_38*z - tmp_49*tmp_7);
    double tmp_53 = tmp_3*tmp_8;
    double tmp_54 = tmp_19*y - tmp_20*y;
    double tmp_55 = tmp_25*tmp_54;
    double tmp_56 = tmp_29*tmp_3;
    double tmp_57 = tmp_3*tmp_32;
    double tmp_58 = tmp_37*tmp_54;
    double tmp_59 = y*z;
    double tmp_60 = tmp_59*tmp_8;
    double tmp_61 = tmp_29*tmp_59;
    double tmp_62 = tmp_32*tmp_59;
    double tmp_63 = -tmp_41*(tmp_13*tmp_60 + tmp_23*(-tmp_16*tmp_59 + tmp_18*tmp_59) - tmp_28*tmp_61 - tmp_34*tmp_62 + tmp_35*tmp_62 + tmp_36*tmp_61 + tmp_55*z - tmp_58*z - tmp_60*tmp_7);
    double tmp_64 = tmp_4*tmp_8;
    double tmp_65 = z*(tmp_19*z - tmp_20*z);
    double tmp_66 = tmp_29*tmp_4;
    double tmp_67 = tmp_32*tmp_4;

    hess[0] = hess[0] + -tmp_41*(tmp_13*tmp_9 + tmp_23*(-tmp_16*tmp_2 + tmp_18*tmp_2 + tmp_21) + tmp_26*x - tmp_28*tmp_30 + tmp_30*tmp_36 - tmp_33*tmp_34 + tmp_33*tmp_35 - tmp_38*x + tmp_40 - tmp_7*tmp_9);
    hess[1] = hess[1] + tmp_47;
    hess[2] = hess[2] + tmp_52;
    hess[3] = hess[3] + tmp_47;
    hess[4] = hess[4] + -tmp_41*(tmp_13*tmp_53 + tmp_23*(-tmp_16*tmp_3 + tmp_18*tmp_3 + tmp_21) - tmp_28*tmp_56 - tmp_34*tmp_57 + tmp_35*tmp_57 + tmp_36*tmp_56 + tmp_40 - tmp_53*tmp_7 + tmp_55*y - tmp_58*y);
    hess[5] = hess[5] + tmp_63;
    hess[6] = hess[6] + tmp_52;
    hess[7] = hess[7] + tmp_63;
    hess[8] = hess[8] + -tmp_41*(tmp_13*tmp_64 + tmp_23*(-tmp_16*tmp_4 + tmp_18*tmp_4 + tmp_21) + tmp_25*tmp_65 - tmp_28*tmp_66 - tmp_34*tmp_67 + tmp_35*tmp_67 + tmp_36*tmp_66 - tmp_37*tmp_65 + tmp_40 - tmp_64*tmp_7);
}

/* ---------------------------------------------------------------------------
    Spherical NFW
*/
double sphericalnfw_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
    */
    double u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = -pars[0] * pars[1] / pars[2];
    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / pars[2];
    return v_h2 * log(1 + u) / u;
}

void sphericalnfw_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
    */
    double fac, u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = pars[0] * pars[1] / pars[2];

    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / pars[2];
    fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2];
}

double sphericalnfw_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
    */
    // double v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    double v_h2 = pars[0] * pars[1] / pars[2];
    double r, rho0;
    r = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);

    rho0 = v_h2 / (4*M_PI*pars[0]*pars[2]*pars[2]);
    return rho0 / ((r/pars[2]) * pow(1+r/pars[2],2));
}

void sphericalnfw_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
    */
    double G = pars[0];
    double m = pars[1];
    double r_s = pars[2];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = tmp_0 + tmp_1 + tmp_2;
    double tmp_4 = pow(tmp_3, 7);
    double tmp_5 = 3*tmp_0;
    double tmp_6 = sqrt(tmp_3);
    double tmp_7 = r_s + tmp_6;
    double tmp_8 = pow(tmp_3, 13.0/2.0)*tmp_7;
    double tmp_9 = pow(tmp_7, 2);
    double tmp_10 = 1.0/r_s;
    double tmp_11 = tmp_9*log(tmp_10*tmp_7);
    double tmp_12 = tmp_11*pow(tmp_3, 6);
    double tmp_13 = tmp_11*tmp_4 - pow(tmp_3, 15.0/2.0)*tmp_7;
    double tmp_14 = G*m;
    double tmp_15 = tmp_14/tmp_9;
    double tmp_16 = tmp_15/pow(tmp_3, 17.0/2.0);
    double tmp_17 = x*y;
    double tmp_18 = 4*tmp_15/pow(tmp_3, 3.0/2.0);
    double tmp_19 = 3*tmp_17;
    double tmp_20 = r_s*tmp_15/pow(tmp_3, 2);
    double tmp_21 = tmp_14*log(tmp_10*tmp_6 + 1)/pow(tmp_3, 5.0/2.0);
    double tmp_22 = tmp_17*tmp_18 + tmp_19*tmp_20 - tmp_19*tmp_21;
    double tmp_23 = x*z;
    double tmp_24 = 3*tmp_20;
    double tmp_25 = 3*tmp_21;
    double tmp_26 = tmp_18*tmp_23 + tmp_23*tmp_24 - tmp_23*tmp_25;
    double tmp_27 = 3*tmp_8;
    double tmp_28 = 3*tmp_12;
    double tmp_29 = y*z;
    double tmp_30 = tmp_18*tmp_29 + tmp_24*tmp_29 - tmp_25*tmp_29;

    hess[0] = hess[0] + tmp_16*(tmp_0*tmp_4 - tmp_12*tmp_5 + tmp_13 + tmp_5*tmp_8);
    hess[1] = hess[1] + tmp_22;
    hess[2] = hess[2] + tmp_26;
    hess[3] = hess[3] + tmp_22;
    hess[4] = hess[4] + tmp_16*(tmp_1*tmp_27 - tmp_1*tmp_28 + tmp_1*tmp_4 + tmp_13);
    hess[5] = hess[5] + tmp_30;
    hess[6] = hess[6] + tmp_26;
    hess[7] = hess[7] + tmp_30;
    hess[8] = hess[8] + tmp_16*(tmp_13 + tmp_2*tmp_27 - tmp_2*tmp_28 + tmp_2*tmp_4);
}

/* ---------------------------------------------------------------------------
    Flattened NFW
*/
double flattenednfw_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (scale mass)
            - r_s (scale radius)
            - qz (flattening)
    */
    double u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = -pars[0] * pars[1] / pars[2];
    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]/(pars[3]*pars[3])) / pars[2];
    return v_h2 * log(1 + u) / u;
}

void flattenednfw_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (scale mass)
            - r_s (scale radius)
            - qz (flattening)
    */
    double fac, u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = pars[0] * pars[1] / pars[2];
    u = sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]/(pars[3]*pars[3])) / pars[2];

    fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2]/(pars[3]*pars[3]);
}

void flattenednfw_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
            - c (flattening)
    */
    double G = pars[0];
    double m = pars[1];
    double r_s = pars[2];
    double c = pars[3];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(z, 2);
    double tmp_2 = pow(c, 2);
    double tmp_3 = pow(y, 2);
    double tmp_4 = tmp_1 + tmp_2*(tmp_0 + tmp_3);
    double tmp_5 = pow(tmp_4, 4);
    double tmp_6 = 3*tmp_0;
    double tmp_7 = tmp_4/tmp_2;
    double tmp_8 = sqrt(tmp_7);
    double tmp_9 = r_s + tmp_8;
    double tmp_10 = pow(c, 8);
    double tmp_11 = tmp_10*tmp_9;
    double tmp_12 = tmp_11*pow(tmp_7, 7.0/2.0);
    double tmp_13 = pow(tmp_4, 3);
    double tmp_14 = pow(tmp_9, 2);
    double tmp_15 = tmp_14*log(tmp_9/r_s);
    double tmp_16 = tmp_15*tmp_2;
    double tmp_17 = -tmp_11*pow(tmp_7, 9.0/2.0) + tmp_15*tmp_5;
    double tmp_18 = G*m/tmp_14;
    double tmp_19 = tmp_18/pow(tmp_7, 11.0/2.0);
    double tmp_20 = tmp_19/tmp_10;
    double tmp_21 = pow(c, 4);
    double tmp_22 = pow(tmp_4, 2);
    double tmp_23 = 3*tmp_9;
    double tmp_24 = pow(tmp_7, 3.0/2.0);
    double tmp_25 = tmp_18*x;
    double tmp_26 = tmp_21*tmp_25*y*(-3*tmp_15*tmp_21*tmp_24 + tmp_21*pow(tmp_7, 5.0/2.0) + tmp_22*tmp_23)/tmp_5;
    double tmp_27 = 3*tmp_16;
    double tmp_28 = tmp_2*z*(tmp_2*tmp_24 + tmp_23*tmp_4 - tmp_27*tmp_8)/tmp_13;
    double tmp_29 = tmp_25*tmp_28;
    double tmp_30 = tmp_18*tmp_28*y;

    hess[0] = hess[0] + tmp_20*(tmp_0*tmp_5 + tmp_12*tmp_6 - tmp_13*tmp_16*tmp_6 + tmp_17);
    hess[1] = hess[1] + tmp_26;
    hess[2] = hess[2] + tmp_29;
    hess[3] = hess[3] + tmp_26;
    hess[4] = hess[4] + tmp_20*(3*tmp_12*tmp_3 - tmp_13*tmp_27*tmp_3 + tmp_17 + tmp_3*tmp_5);
    hess[5] = hess[5] + tmp_30;
    hess[6] = hess[6] + tmp_29;
    hess[7] = hess[7] + tmp_30;
    hess[8] = hess[8] + tmp_13*tmp_19*(tmp_1*tmp_2*tmp_23*tmp_8 - tmp_1*tmp_27 + tmp_1*tmp_4 + tmp_16*tmp_4 - tmp_22*tmp_9/tmp_8)/pow(c, 12);

}

/* ---------------------------------------------------------------------------
    Triaxial NFW - triaxiality in potential!
*/
double triaxialnfw_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (scale mass)
            - r_s (scale radius)
            - a (major axis)
            - b (intermediate axis)
            - c (minor axis)
    */
    double u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = -pars[0] * pars[1] / pars[2];
    u = sqrt(q[0]*q[0]/(pars[3]*pars[3])
           + q[1]*q[1]/(pars[4]*pars[4])
           + q[2]*q[2]/(pars[5]*pars[5])) / pars[2];
    return v_h2 * log(1 + u) / u;
}

void triaxialnfw_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - v_c (circular velocity at the scale radius)
            - r_s (scale radius)
            - a (major axis)
            - b (intermediate axis)
            - c (minor axis)
    */
    double fac, u, v_h2;
    // v_h2 = pars[1]*pars[1] / (log(2.) - 0.5);
    v_h2 = pars[0] * pars[1] / pars[2];
    u = sqrt(q[0]*q[0]/(pars[3]*pars[3])
           + q[1]*q[1]/(pars[4]*pars[4])
           + q[2]*q[2]/(pars[5]*pars[5])) / pars[2];

    fac = v_h2 / (u*u*u) / (pars[2]*pars[2]) * (log(1+u) - u/(1+u));

    grad[0] = grad[0] + fac*q[0]/(pars[3]*pars[3]);
    grad[1] = grad[1] + fac*q[1]/(pars[4]*pars[4]);
    grad[2] = grad[2] + fac*q[2]/(pars[5]*pars[5]);
}

void triaxialnfw_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - r_s (scale radius)
            - a (major axis)
            - b (intermediate axis)
            - c (minor axis)
    */
    double G = pars[0];
    double m = pars[1];
    double r_s = pars[2];
    double a = pars[3];
    double b = pars[4];
    double c = pars[5];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(a, -2);
    double tmp_1 = G*m;
    double tmp_2 = tmp_0*tmp_1;
    double tmp_3 = pow(x, 2);
    double tmp_4 = pow(b, -2);
    double tmp_5 = pow(y, 2);
    double tmp_6 = pow(c, -2);
    double tmp_7 = pow(z, 2);
    double tmp_8 = tmp_0*tmp_3 + tmp_4*tmp_5 + tmp_6*tmp_7;
    double tmp_9 = pow(tmp_8, -3.0/2.0);
    double tmp_10 = 1.0/r_s;
    double tmp_11 = tmp_10*sqrt(tmp_8) + 1;
    double tmp_12 = log(tmp_11);
    double tmp_13 = tmp_12*tmp_9;
    double tmp_14 = tmp_3/pow(a, 4);
    double tmp_15 = 3*tmp_1;
    double tmp_16 = tmp_12/pow(tmp_8, 5.0/2.0);
    double tmp_17 = tmp_15*tmp_16;
    double tmp_18 = tmp_10/tmp_11;
    double tmp_19 = tmp_18/tmp_8;
    double tmp_20 = tmp_9/(pow(r_s, 2)*pow(tmp_11, 2));
    double tmp_21 = tmp_1*tmp_20;
    double tmp_22 = tmp_18/pow(tmp_8, 2);
    double tmp_23 = tmp_15*tmp_22;
    double tmp_24 = tmp_4*y;
    double tmp_25 = tmp_2*x;
    double tmp_26 = 3*tmp_25;
    double tmp_27 = tmp_16*tmp_26;
    double tmp_28 = tmp_20*tmp_25;
    double tmp_29 = tmp_22*tmp_26;
    double tmp_30 = -tmp_24*tmp_27 + tmp_24*tmp_28 + tmp_24*tmp_29;
    double tmp_31 = tmp_6*z;
    double tmp_32 = -tmp_27*tmp_31 + tmp_28*tmp_31 + tmp_29*tmp_31;
    double tmp_33 = tmp_1*tmp_13;
    double tmp_34 = tmp_5/pow(b, 4);
    double tmp_35 = tmp_1*tmp_19;
    double tmp_36 = tmp_24*tmp_31;
    double tmp_37 = -tmp_17*tmp_36 + tmp_21*tmp_36 + tmp_23*tmp_36;
    double tmp_38 = tmp_7/pow(c, 4);

    hess[0] = hess[0] + tmp_13*tmp_2 - tmp_14*tmp_17 + tmp_14*tmp_21 + tmp_14*tmp_23 - tmp_19*tmp_2;
    hess[1] = hess[1] + tmp_30;
    hess[2] = hess[2] + tmp_32;
    hess[3] = hess[3] + tmp_30;
    hess[4] = hess[4] + -tmp_17*tmp_34 + tmp_21*tmp_34 + tmp_23*tmp_34 + tmp_33*tmp_4 - tmp_35*tmp_4;
    hess[5] = hess[5] + tmp_37;
    hess[6] = hess[6] + tmp_32;
    hess[7] = hess[7] + tmp_37;
    hess[8] = hess[8] + -tmp_17*tmp_38 + tmp_21*tmp_38 + tmp_23*tmp_38 + tmp_33*tmp_6 - tmp_35*tmp_6;
}

/* ---------------------------------------------------------------------------
    Satoh potential
*/
double satoh_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */
    double S2;
    S2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
    return -pars[0] * pars[1] / sqrt(S2);
}

void satoh_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */

    double S2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + pars[2]*(pars[2] + 2*sqrt(q[2]*q[2] + pars[3]*pars[3]));
    double dPhi_dS = pars[0] * pars[1] / S2;

    grad[0] = grad[0] + dPhi_dS*q[0]/sqrt(S2);
    grad[1] = grad[1] + dPhi_dS*q[1]/sqrt(S2);
    grad[2] = grad[2] + dPhi_dS/sqrt(S2) * q[2]*(1 + pars[2] / sqrt(q[2]*q[2] + pars[3]*pars[3]));
}

double satoh_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */
    double z2b2 = q[2]*q[2] + pars[3]*pars[3];
    double xyz2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2];
    double S2 = xyz2 + pars[2]*(pars[2] + 2*sqrt(z2b2));
    double A = pars[1] * pars[2] * pars[3]*pars[3] / (4*M_PI*S2*sqrt(S2)*z2b2);
    return A * (1/sqrt(z2b2) + 3/pars[2]*(1 - xyz2/S2));
}

void satoh_hessian(double t, double *pars, double *q, int n_dim, double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a ()
            - b ()
    */
    double G = pars[0];
    double m = pars[1];
    double a = pars[2];
    double b = pars[3];

    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = pow(b, 2) + tmp_2;
    double tmp_4 = sqrt(tmp_3);
    double tmp_5 = a*(a + 2*tmp_4) + tmp_0 + tmp_1 + tmp_2;
    double tmp_6 = G*m;
    double tmp_7 = tmp_6/pow(tmp_5, 3.0/2.0);
    double tmp_8 = tmp_6/pow(tmp_5, 5.0/2.0);
    double tmp_9 = 3*tmp_8;
    double tmp_10 = -tmp_9*x*y;
    double tmp_11 = 3*z;
    double tmp_12 = a/tmp_4;
    double tmp_13 = tmp_8*(-tmp_11*tmp_12 - tmp_11);
    double tmp_14 = tmp_13*x;
    double tmp_15 = tmp_13*y;

    hess[0] = hess[0] + -tmp_0*tmp_9 + tmp_7;
    hess[1] = hess[1] + tmp_10;
    hess[2] = hess[2] + tmp_14;
    hess[3] = hess[3] + tmp_10;
    hess[4] = hess[4] + -tmp_1*tmp_9 + tmp_7;
    hess[5] = hess[5] + tmp_15;
    hess[6] = hess[6] + tmp_14;
    hess[7] = hess[7] + tmp_15;
    hess[8] = hess[8] + -tmp_13*(-tmp_12*z - z) - tmp_7*(a*tmp_2/pow(tmp_3, 3.0/2.0) - tmp_12 - 1);
}

/* ---------------------------------------------------------------------------
    Kuzmin potential
*/
double kuzmin_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
    */
    double S2 = q[0]*q[0] + q[1]*q[1] + pow(pars[2] + fabs(q[2]), 2);
    return -pars[0] * pars[1] / sqrt(S2);
}

void kuzmin_gradient(double t, double *pars, double *q, int n_dim,
                     double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
    */

    double S2 = q[0]*q[0] + q[1]*q[1] + pow(pars[2] + fabs(q[2]), 2);
    double fac = pars[0] * pars[1] * pow(S2, -1.5);
    double zsign;

    if (q[2] > 0) {
        zsign = 1.;
    } else if (q[2] < 0) {
        zsign = -1.;
    } else {
        zsign = 0.;
    }

    grad[0] = grad[0] + fac * q[0];
    grad[1] = grad[1] + fac * q[1];
    grad[2] = grad[2] + fac * zsign * (pars[2] + fabs(q[2]));
}

double kuzmin_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
    */
    if (q[2] != 0.) {
        return 0.;
    } else {
        return pars[1] * pars[2] / (2 * M_PI) *
            pow(q[0]*q[0] + q[1]*q[1] + pars[2]*pars[2], -1.5);
   }

}

/* ---------------------------------------------------------------------------
    Miyamoto-Nagai flattened potential
*/
double miyamotonagai_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */
    double zd;
    zd = (pars[2] + sqrt(q[2]*q[2] + pars[3]*pars[3]));
    return -pars[0] * pars[1] / sqrt(q[0]*q[0] + q[1]*q[1] + zd*zd);
}

void miyamotonagai_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */
    double sqrtz, zd, fac;

    sqrtz = sqrt(q[2]*q[2] + pars[3]*pars[3]);
    zd = pars[2] + sqrtz;
    fac = pars[0]*pars[1] * pow(q[0]*q[0] + q[1]*q[1] + zd*zd, -1.5);

    grad[0] = grad[0] + fac*q[0];
    grad[1] = grad[1] + fac*q[1];
    grad[2] = grad[2] + fac*q[2] * (1. + pars[2] / sqrtz);
}

double miyamotonagai_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */

    double M, a, b;
    M = pars[1];
    a = pars[2];
    b = pars[3];

    double R2 = q[0]*q[0] + q[1]*q[1];
    double sqrt_zb = sqrt(q[2]*q[2] + b*b);
    double numer = (b*b*M / (4*M_PI)) * (a*R2 + (a + 3*sqrt_zb)*(a + sqrt_zb)*(a + sqrt_zb));
    double denom = pow(R2 + (a + sqrt_zb)*(a + sqrt_zb), 2.5) * sqrt_zb*sqrt_zb*sqrt_zb;

    return numer/denom;
}

void miyamotonagai_hessian(double t, double *pars, double *q, int n_dim,
                           double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - m (mass scale)
            - a (length scale 1) TODO
            - b (length scale 2) TODO
    */
    double G = pars[0];
    double m = pars[1];
    double a = pars[2];
    double b = pars[3];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(x, 2);
    double tmp_1 = pow(y, 2);
    double tmp_2 = pow(z, 2);
    double tmp_3 = pow(b, 2) + tmp_2;
    double tmp_4 = sqrt(tmp_3);
    double tmp_5 = a + tmp_4;
    double tmp_6 = pow(tmp_5, 2);
    double tmp_7 = tmp_0 + tmp_1 + tmp_6;
    double tmp_8 = G*m;
    double tmp_9 = tmp_8/pow(tmp_7, 3.0/2.0);
    double tmp_10 = 3*tmp_8/pow(tmp_7, 5.0/2.0);
    double tmp_11 = tmp_10*x;
    double tmp_12 = -tmp_11*y;
    double tmp_13 = tmp_5/tmp_4;
    double tmp_14 = tmp_13*z;
    double tmp_15 = -tmp_11*tmp_14;
    double tmp_16 = -tmp_10*tmp_14*y;
    double tmp_17 = 1.0/tmp_3;
    double tmp_18 = tmp_2*tmp_9;

    hess[0] = hess[0] + -tmp_0*tmp_10 + tmp_9;
    hess[1] = hess[1] + tmp_12;
    hess[2] = hess[2] + tmp_15;
    hess[3] = hess[3] + tmp_12;
    hess[4] = hess[4] + -tmp_1*tmp_10 + tmp_9;
    hess[5] = hess[5] + tmp_16;
    hess[6] = hess[6] + tmp_15;
    hess[7] = hess[7] + tmp_16;
    hess[8] = hess[8] + -tmp_10*tmp_17*tmp_2*tmp_6 + tmp_13*tmp_9 + tmp_17*tmp_18 - tmp_18*tmp_5/pow(tmp_3, 3.0/2.0);
}

/* ---------------------------------------------------------------------------
    Lee-Suto triaxial NFW from Lee & Suto (2003)
*/
double leesuto_value(double t, double *pars, double *q, int n_dim) {
    /*  pars: (alpha = 1)
            0 - G
            1 - v_c
            2 - r_s
            3 - a
            4 - b
            5 - c
    */
    double x, y, z, _r, u, phi0;
    double e_b2 = 1-pow(pars[4]/pars[3],2);
    double e_c2 = 1-pow(pars[5]/pars[3],2);
    double F1,F2,F3,costh2,sinth2,sinph2;

    phi0 = pars[1]*pars[1] / (log(2.) - 0.5 + (log(2.)-0.75)*e_b2 + (log(2.)-0.75)*e_c2);

    x = q[0];
    y = q[1];
    z = q[2];

    _r = sqrt(x*x + y*y + z*z);
    u = _r / pars[2];

    F1 = -log(1+u)/u;
    F2 = -1/3. + (2*u*u - 3*u + 6)/(6*u*u) + (1/u - pow(u,-3.))*log(1+u);
    F3 = (u*u - 3*u - 6)/(2*u*u*(1+u)) + 3*pow(u,-3)*log(1+u);
    costh2 = z*z / (_r*_r);
    sinth2 = 1 - costh2;
    sinph2 = y*y / (x*x + y*y);
    //return phi0 * ((e_b2/2 + e_c2/2)*((1/u - 1/(u*u*u))*log(u + 1) - 1 + (2*u*u - 3*u + 6)/(6*u*u)) + (e_b2*y*y/(2*_r*_r) + e_c2*z*z/(2*_r*_r))*((u*u - 3*u - 6)/(2*u*u*(u + 1)) + 3*log(u + 1)/(u*u*u)) - log(u + 1)/u);
    return phi0 * (F1 + (e_b2+e_c2)/2.*F2 + (e_b2*sinth2*sinph2 + e_c2*costh2)/2. * F3);
}

void leesuto_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars: (alpha = 1)
            0 - G
            1 - v_c
            2 - r_s
            3 - a
            4 - b
            5 - c
    */
    double x, y, z, _r, _r2, _r4, ax, ay, az;
    double v_h2, x0, x2, x22;
    double x20, x21, x7, x1;
    double x10, x13, x15, x16, x17;
    double e_b2 = 1-pow(pars[4]/pars[3],2);
    double e_c2 = 1-pow(pars[5]/pars[3],2);

    v_h2 = pars[1]*pars[1] / (log(2.) - 0.5 + (log(2.)-0.75)*e_b2 + (log(2.)-0.75)*e_c2);

    x = q[0];
    y = q[1];
    z = q[2];

    _r2 = x*x + y*y + z*z;
    _r = sqrt(_r2);
    _r4 = _r2*_r2;

    x0 = _r + pars[2];
    x1 = x0*x0;
    x2 = v_h2/(12.*_r4*_r2*_r*x1);
    x10 = log(x0/pars[2]);

    x13 = _r*3.*pars[2];
    x15 = x13 - _r2;
    x16 = x15 + 6.*(pars[2]*pars[2]);
    x17 = 6.*pars[2]*x0*(_r*x16 - x0*x10*6.*(pars[2]*pars[2]));
    x20 = x0*_r2;
    x21 = 2.*_r*x0;
    x7 = e_b2*y*y + e_c2*z*z;
    x22 = -12.*_r4*_r*pars[2]*x0 + 12.*_r4*pars[2]*x1*x10 + 3.*pars[2]*x7*(x16*_r2 - 18.*x1*x10*(pars[2]*pars[2]) + x20*(2.*_r - 3.*pars[2]) + x21*(x15 + 9.*(pars[2]*pars[2]))) - x20*(e_b2 + e_c2)*(-6.*_r*pars[2]*(_r2 - (pars[2]*pars[2])) + 6.*pars[2]*x0*x10*(_r2 - 3.*(pars[2]*pars[2])) + x20*(-4.*_r + 3.*pars[2]) + x21*(-x13 + 2.*_r2 + 6.*(pars[2]*pars[2])));

    ax = x2*x*(x17*x7 + x22);
    ay = x2*y*(x17*(x7 - _r2*e_b2) + x22);
    az = x2*z*(x17*(x7 - _r2*e_c2) + x22);

    grad[0] = grad[0] + ax;
    grad[1] = grad[1] + ay;
    grad[2] = grad[2] + az;
}

double leesuto_density(double t, double *pars, double *q, int n_dim) {
    /*  pars: (alpha = 1)
            0 - G
            1 - v_c
            2 - r_s
            3 - a
            4 - b
            5 - c
    */
    double x, y, z, u, v_h2;
    double b_a2, c_a2;
    b_a2 = pars[4]*pars[4] / (pars[3]*pars[3]);
    c_a2 = pars[5]*pars[5] / (pars[3]*pars[3]);
    double e_b2 = 1-b_a2;
    double e_c2 = 1-c_a2;
    v_h2 = pars[1]*pars[1] / (log(2.) - 0.5 + (log(2.)-0.75)*e_b2 + (log(2.)-0.75)*e_c2);

    x = q[0];
    y = q[1];
    z = q[2];

    u = sqrt(x*x + y*y/b_a2 + z*z/c_a2) / pars[2];
    return v_h2 / (u * (1+u)*(1+u)) / (4.*M_PI*pars[2]*pars[2]*pars[0]);
}

/* ---------------------------------------------------------------------------
    Logarithmic (triaxial)
*/
double logarithmic_value(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - v_c (velocity scale)
            - r_h (length scale)
            - q1
            - q2
            - q3
    */
    double x, y, z;

    x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
    y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
    z = q[2];

    return 0.5*pars[1]*pars[1] * log(pars[2]*pars[2] + // scale radius
                                     x*x/(pars[3]*pars[3]) +
                                     y*y/(pars[4]*pars[4]) +
                                     z*z/(pars[5]*pars[5]));
}

double logarithmic_density(double t, double *pars, double *q, int n_dim) {
    /*  pars:
            - G (Gravitational constant)
            - v_c (velocity scale)
            - r_h (length scale)
            - q1
            - q2
            - q3
    */
    double tmp_0 = pow(pars[3], 2);
    double tmp_1 = pow(pars[4], 2);
    double tmp_2 = tmp_0*tmp_1;
    double tmp_3 = tmp_2*pow(q[2], 2);
    double tmp_4 = pow(pars[5], 2);
    double tmp_5 = tmp_0*tmp_4;
    double tmp_6 = tmp_5*pow(q[1], 2);
    double tmp_7 = tmp_1*tmp_4;
    double tmp_8 = tmp_7*pow(q[0], 2);
    double tmp_9 = pow(pars[2], 2)*tmp_2*tmp_4;
    double tmp_10 = tmp_6 + tmp_8 + tmp_9;
    double tmp_11 = tmp_3 + tmp_9;
    return pow(pars[1], 2)*(tmp_2*(tmp_10 - tmp_3) + tmp_5*(tmp_11 - tmp_6 + tmp_8) + tmp_7*(tmp_11 + tmp_6 - tmp_8))/pow(tmp_10 + tmp_3, 2) / (4*M_PI*pars[0]);
}

void logarithmic_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  pars:
            - G (Gravitational constant)
            - v_c (velocity scale)
            - r_h (length scale)
            - q1
            - q2
            - q3
    */
    double x, y, z, ax, ay, az, fac;

    x = q[0]*cos(pars[6]) + q[1]*sin(pars[6]);
    y = -q[0]*sin(pars[6]) + q[1]*cos(pars[6]);
    z = q[2];

    fac = pars[1]*pars[1] / (pars[2]*pars[2] + x*x/(pars[3]*pars[3]) + y*y/(pars[4]*pars[4]) + z*z/(pars[5]*pars[5]));
    ax = fac*x/(pars[3]*pars[3]);
    ay = fac*y/(pars[4]*pars[4]);
    az = fac*z/(pars[5]*pars[5]);

    grad[0] = grad[0] + (ax*cos(pars[6]) - ay*sin(pars[6]));
    grad[1] = grad[1] + (ax*sin(pars[6]) + ay*cos(pars[6]));
    grad[2] = grad[2] + az;
}

void logarithmic_hessian(double t, double *pars, double *q, int n_dim,
                         double *hess) {
    /*  pars:
            - G (Gravitational constant)
            - v_c (velocity scale)
            - r_h (length scale)
            - q1
            - q2
            - q3
    */
    double v_c = pars[1];
    double r_h = pars[2];
    double q1 = pars[3];
    double q2 = pars[4];
    double q3 = pars[5];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = pow(q1, -2);
    double tmp_1 = pow(v_c, 2);
    double tmp_2 = tmp_0*tmp_1;
    double tmp_3 = pow(x, 2);
    double tmp_4 = pow(q2, -2);
    double tmp_5 = pow(y, 2);
    double tmp_6 = pow(q3, -2);
    double tmp_7 = pow(z, 2);
    double tmp_8 = pow(r_h, 2) + tmp_0*tmp_3 + tmp_4*tmp_5 + tmp_6*tmp_7;
    double tmp_9 = 1.0/tmp_8;
    double tmp_10 = 2.0/pow(tmp_8, 2);
    double tmp_11 = tmp_1*tmp_10;
    double tmp_12 = tmp_4*y;
    double tmp_13 = tmp_10*tmp_2*x;
    double tmp_14 = tmp_12*tmp_13;
    double tmp_15 = tmp_6*z;
    double tmp_16 = tmp_13*tmp_15;
    double tmp_17 = tmp_1*tmp_9;
    double tmp_18 = tmp_11*tmp_12*tmp_15;

    // minus signs because I initially borked the sympy definition
    hess[0] = hess[0] - (-tmp_2*tmp_9 + tmp_11*tmp_3/pow(q1, 4));
    hess[1] = hess[1] - (tmp_14);
    hess[2] = hess[2] - (tmp_16);
    hess[3] = hess[3] - (tmp_14);
    hess[4] = hess[4] - (-tmp_17*tmp_4 + tmp_11*tmp_5/pow(q2, 4));
    hess[5] = hess[5] - (tmp_18);
    hess[6] = hess[6] - (tmp_16);
    hess[7] = hess[7] - (tmp_18);
    hess[8] = hess[8] - (-tmp_17*tmp_6 + tmp_11*tmp_7/pow(q3, 4));
}

/* ---------------------------------------------------------------------------
    Logarithmic (triaxial)
*/
double longmuralibar_value(double t, double *pars, double *q, int n_dim) {
    /*  http://adsabs.harvard.edu/abs/1992ApJ...397...44L

        pars:
        - G (Gravitational constant)
        - m (mass scale)
        - a
        - b
        - c
        - alpha
    */
    double x, y, z;
    double a, b, c;
    double Tm, Tp;

    x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
    y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
    z = q[2];

    a = pars[2];
    b = pars[3];
    c = pars[4];

    Tm = sqrt((a-x)*(a-x) + y*y + pow(b + sqrt(c*c + z*z),2));
    Tp = sqrt((a+x)*(a+x) + y*y + pow(b + sqrt(c*c + z*z),2));

    return pars[0]*pars[1]/(2*a) * log((x - a + Tm) / (x + a + Tp));
}

void longmuralibar_gradient(double t, double *pars, double *q, int n_dim, double *grad) {
    /*  http://adsabs.harvard.edu/abs/1992ApJ...397...44L

        pars:
        - G (Gravitational constant)
        - m (mass scale)
        - a
        - b
        - c
        - alpha
    */
    double x, y, z;
    double a, b, c;
    double Tm, Tp, fac1, fac2, fac3, bcz;
    double gx, gy, gz;

    x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
    y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
    z = q[2];

    a = pars[2];
    b = pars[3];
    c = pars[4];

    bcz = b + sqrt(c*c + z*z);
    Tm = sqrt((a-x)*(a-x) + y*y + bcz*bcz);
    Tp = sqrt((a+x)*(a+x) + y*y + bcz*bcz);

    fac1 = pars[0]*pars[1] / (2*Tm*Tp);
    fac2 = 1 / (y*y + bcz*bcz);
    fac3 = Tp + Tm - (4*x*x)/(Tp+Tm);

    gx = 4 * fac1 * x / (Tp + Tm);
    gy = fac1 * y * fac2 * fac3;
    gz = fac1 * z * fac2 * fac3 * bcz / sqrt(z*z + c*c);

    grad[0] = grad[0] + (gx*cos(pars[5]) - gy*sin(pars[5]));
    grad[1] = grad[1] + (gx*sin(pars[5]) + gy*cos(pars[5]));
    grad[2] = grad[2] + gz;
}

double longmuralibar_density(double t, double *pars, double *q, int n_dim) {
    /*
        Generated by sympy...

        pars:
        - G (Gravitational constant)
        - m (mass scale)
        - a
        - b
        - c
        - alpha
    */
    double a = pars[2];
    double b = pars[3];
    double c = pars[4];

    double x = q[0]*cos(pars[5]) + q[1]*sin(pars[5]);
    double y = -q[0]*sin(pars[5]) + q[1]*cos(pars[5]);
    double z = q[2];

    double tmp0 = a - x;
    double tmp1 = pow(tmp0, 2);
    double tmp2 = pow(y, 2);
    double tmp3 = pow(z, 2);
    double tmp4 = pow(c, 2) + tmp3;
    double tmp5 = sqrt(tmp4);
    double tmp6 = b + tmp5;
    double tmp7 = pow(tmp6, 2);
    double tmp8 = tmp2 + tmp7;
    double tmp9 = tmp1 + tmp8;
    double tmp10 = sqrt(tmp9);
    double tmp11 = -a + tmp10 + x;
    double tmp12 = 1.0/tmp11;
    double tmp13 = 1.0/tmp10;
    double tmp14 = pow(tmp9, -1.5);
    double tmp15 = 1.0/tmp4;
    double tmp16 = tmp13*tmp3;
    double tmp17 = tmp6/tmp5;
    double tmp18 = pow(tmp4, -1.5);
    double tmp19 = tmp15*tmp3*tmp7;
    double tmp20 = 2*tmp2;
    double tmp21 = a + x;
    double tmp22 = pow(tmp21, 2);
    double tmp23 = tmp22 + tmp8;
    double tmp24 = sqrt(tmp23);
    double tmp25 = 1.0/tmp24;
    double tmp26 = tmp21 + tmp24;
    double tmp27 = 1.0/tmp26;
    double tmp28 = tmp25*tmp27;
    double tmp29 = tmp11*tmp28;
    double tmp30 = tmp11*tmp27/pow(tmp23, 1.5);
    double tmp31 = 1.0/tmp23;
    double tmp32 = pow(tmp26, -2);
    double tmp33 = tmp11*tmp31*tmp32;
    double tmp34 = tmp21*tmp25 + 1;
    double tmp35 = tmp27*tmp34;
    double tmp36 = tmp13*tmp15*tmp3*tmp7;
    double tmp37 = -tmp13 + tmp29;
    double tmp38 = tmp2*tmp37;
    double tmp39 = tmp0*tmp13;
    double tmp40 = tmp11*tmp27*tmp34 + tmp39 - 1;
    return pars[1]/8.*tmp12*(2*tmp11*tmp32*pow(tmp34, 2) +
        tmp12*tmp13*tmp38 + tmp12*tmp36*tmp37 + tmp12*tmp40*(-tmp39 + 1) +
        tmp13*tmp17 - tmp13*tmp20*tmp25*tmp27 + tmp13*(-tmp1/tmp9 + 1) + tmp13 -
        tmp14*tmp19 - tmp14*tmp2 + tmp15*tmp16 - tmp15*tmp28*tmp3*tmp37*tmp7 -
        tmp15*tmp29*tmp3 + 2*tmp15*tmp3*tmp33*tmp7 - tmp16*tmp18*tmp6 -
        tmp17*tmp29 + tmp18*tmp29*tmp3*tmp6 + tmp19*tmp30 + tmp2*tmp30 +
        tmp20*tmp33 - 2*tmp25*tmp27*tmp36 - tmp28*tmp38 - tmp29*(-tmp22*tmp31 +
        1) - tmp29 - tmp35*tmp40 - tmp35*(-2*tmp0*tmp13 + 2))/(M_PI*a);
}

void longmuralibar_hessian(double t, double *pars, double *q, int n_dim,
                         double *hess) {
    /* Generated by sympy...

        pars:
        - G (Gravitational constant)
        - m (mass scale)
        - a
        - b
        - c
        - alpha
    */
    double G = pars[0];
    double m = pars[1];
    double a = pars[2];
    double b = pars[3];
    double c = pars[4];
    double alpha = pars[5];
    double x = q[0];
    double y = q[1];
    double z = q[2];

    double tmp_0 = cos(alpha);
    double tmp_1 = tmp_0*x;
    double tmp_2 = sin(alpha);
    double tmp_3 = tmp_2*y;
    double tmp_4 = tmp_1 + tmp_3;
    double tmp_5 = a + tmp_4;
    double tmp_6 = tmp_0*tmp_5;
    double tmp_7 = tmp_0*y - tmp_2*x;
    double tmp_8 = tmp_2*tmp_7;
    double tmp_9 = -tmp_8;
    double tmp_10 = tmp_6 + tmp_9;
    double tmp_11 = pow(z, 2);
    double tmp_12 = pow(c, 2) + tmp_11;
    double tmp_13 = sqrt(tmp_12);
    double tmp_14 = b + tmp_13;
    double tmp_15 = pow(tmp_14, 2);
    double tmp_16 = tmp_15 + pow(tmp_7, 2);
    double tmp_17 = tmp_16 + pow(tmp_5, 2);
    double tmp_18 = sqrt(tmp_17);
    double tmp_19 = 1.0/tmp_18;
    double tmp_20 = tmp_10*tmp_19;
    double tmp_21 = tmp_18 + tmp_5;
    double tmp_22 = 1.0/tmp_21;
    double tmp_23 = a - tmp_1 - tmp_3;
    double tmp_24 = tmp_0*tmp_23;
    double tmp_25 = -tmp_24 + tmp_9;
    double tmp_26 = tmp_16 + pow(tmp_23, 2);
    double tmp_27 = sqrt(tmp_26);
    double tmp_28 = 1.0/tmp_27;
    double tmp_29 = tmp_25*tmp_28;
    double tmp_30 = tmp_0 + tmp_29;
    double tmp_31 = -tmp_0;
    double tmp_32 = -tmp_20 + tmp_31;
    double tmp_33 = pow(tmp_21, -2);
    double tmp_34 = -a + tmp_27 + tmp_4;
    double tmp_35 = tmp_33*tmp_34;
    double tmp_36 = tmp_22*tmp_30 + tmp_32*tmp_35;
    double tmp_37 = (1.0/2.0)*G*m/a;
    double tmp_38 = tmp_37/tmp_34;
    double tmp_39 = tmp_36*tmp_38;
    double tmp_40 = tmp_21*tmp_37/pow(tmp_34, 2);
    double tmp_41 = tmp_36*tmp_40;
    double tmp_42 = tmp_32*tmp_33;
    double tmp_43 = pow(tmp_0, 2) + pow(tmp_2, 2);
    double tmp_44 = tmp_28*tmp_43;
    double tmp_45 = pow(tmp_26, -3.0/2.0);
    double tmp_46 = tmp_25*tmp_45;
    double tmp_47 = -tmp_19*tmp_43;
    double tmp_48 = pow(tmp_17, -3.0/2.0);
    double tmp_49 = tmp_10*tmp_48;
    double tmp_50 = tmp_34/pow(tmp_21, 3);
    double tmp_51 = tmp_32*tmp_50;
    double tmp_52 = tmp_21*tmp_38;
    double tmp_53 = tmp_0*tmp_7;
    double tmp_54 = tmp_2*tmp_5;
    double tmp_55 = tmp_53 + tmp_54;
    double tmp_56 = tmp_19*tmp_55;
    double tmp_57 = tmp_2 + tmp_56;
    double tmp_58 = -tmp_2;
    double tmp_59 = tmp_2*tmp_23;
    double tmp_60 = tmp_53 - tmp_59;
    double tmp_61 = tmp_28*tmp_60;
    double tmp_62 = tmp_58 - tmp_61;
    double tmp_63 = -tmp_53;
    double tmp_64 = tmp_59 + tmp_63;
    double tmp_65 = tmp_22*tmp_46;
    double tmp_66 = -tmp_56 + tmp_58;
    double tmp_67 = tmp_33*tmp_66;
    double tmp_68 = tmp_2 + tmp_61;
    double tmp_69 = -tmp_54 + tmp_63;
    double tmp_70 = tmp_35*tmp_49;
    double tmp_71 = -2*tmp_2 - 2*tmp_56;
    double tmp_72 = tmp_39*tmp_57 + tmp_41*tmp_62 + tmp_52*(tmp_30*tmp_67 + tmp_42*tmp_68 + tmp_51*tmp_71 + tmp_64*tmp_65 - tmp_69*tmp_70);
    double tmp_73 = 1.0/tmp_13;
    double tmp_74 = tmp_14*tmp_73;
    double tmp_75 = tmp_74*z;
    double tmp_76 = tmp_19*tmp_75;
    double tmp_77 = tmp_28*tmp_75;
    double tmp_78 = tmp_19*tmp_33;
    double tmp_79 = tmp_75*tmp_78;
    double tmp_80 = 2*tmp_76;
    double tmp_81 = tmp_39*tmp_76 - tmp_41*tmp_77 + tmp_52*(-tmp_30*tmp_79 + tmp_42*tmp_77 - tmp_51*tmp_80 - tmp_65*tmp_75 + tmp_70*tmp_75);
    double tmp_82 = tmp_22*tmp_68 + tmp_35*tmp_66;
    double tmp_83 = tmp_38*tmp_82;
    double tmp_84 = tmp_40*tmp_82;
    double tmp_85 = tmp_45*tmp_60;
    double tmp_86 = tmp_50*tmp_66;
    double tmp_87 = tmp_48*tmp_55;
    double tmp_88 = tmp_52*(-tmp_22*tmp_75*tmp_85 + tmp_35*tmp_75*tmp_87 + tmp_67*tmp_77 - tmp_68*tmp_79 - tmp_80*tmp_86) + tmp_76*tmp_83 - tmp_77*tmp_84;
    double tmp_89 = tmp_22*tmp_28;
    double tmp_90 = tmp_14*tmp_89;
    double tmp_91 = tmp_73*tmp_90;
    double tmp_92 = tmp_19*tmp_35;
    double tmp_93 = tmp_74*tmp_92;
    double tmp_94 = tmp_91*z - tmp_93*z;
    double tmp_95 = tmp_11/tmp_12;
    double tmp_96 = tmp_11/pow(tmp_12, 3.0/2.0);
    double tmp_97 = tmp_15*tmp_95;
    double tmp_98 = 2*tmp_97;

    hess[0] = hess[0] + tmp_39*(tmp_0 + tmp_20) + tmp_41*(-tmp_29 + tmp_31) + tmp_52*(tmp_22*(tmp_44 + tmp_46*(tmp_24 + tmp_8)) + 2*tmp_30*tmp_42 + tmp_35*(tmp_47 - tmp_49*(-tmp_6 + tmp_8)) + tmp_51*(-2*tmp_0 - 2*tmp_20));
    hess[1] = hess[1] + tmp_72;
    hess[2] = hess[2] + tmp_81;
    hess[3] = hess[3] + tmp_72;
    hess[4] = hess[4] + tmp_52*(tmp_22*(tmp_44 + tmp_64*tmp_85) + tmp_35*(tmp_47 - tmp_69*tmp_87) + 2*tmp_67*tmp_68 + tmp_71*tmp_86) + tmp_57*tmp_83 + tmp_62*tmp_84;
    hess[5] = hess[5] + tmp_88;
    hess[6] = hess[6] + tmp_81;
    hess[7] = hess[7] + tmp_88;
    hess[8] = hess[8] + tmp_38*tmp_76*tmp_94 - tmp_40*tmp_77*tmp_94 + tmp_52*(tmp_14*tmp_92*tmp_96 - tmp_22*tmp_45*tmp_97 - tmp_28*tmp_78*tmp_98 + tmp_35*tmp_48*tmp_97 + tmp_89*tmp_95 - tmp_90*tmp_96 + tmp_91 - tmp_92*tmp_95 - tmp_93 + tmp_50*tmp_98/tmp_17);
}

max_stars_repo_name: akeemlh/gala

max_issues_repo_name: akeemlh/gala

max_forks_repo_name: akeemlh/gala


text:
#ifndef INCLUDED_BLAS_L3_SYRK_H
#define INCLUDED_BLAS_L3_SYRK_H

#include <complex>
#include <cblas.h>

namespace blas{

    void syrk( CBLAS_ORDER order, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans,
        std::size_t n, std::size_t k, float alpha, const float* A,
        std::size_t lda, float beta, float* C, std::size_t ldc ){
        return cblas_ssyrk( order, uplo, trans, n, k, alpha, A, lda, beta, C,
            ldc );
    }
    void syrk( CBLAS_ORDER order, CBLAS_UPLO uplo,
        CBLAS_TRANSPOSE trans, std::size_t n, std::size_t k, double alpha,
        const double* A, std::size_t lda, double beta, double* C,
        std::size_t ldc ){
        return cblas_dsyrk( order, uplo, trans, n, k, alpha, A, lda, beta, C,
            ldc );
    }
    void syrk( CBLAS_ORDER order, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans,
        std::size_t n, std::size_t k, const std::complex<float>* alpha,
        const std::complex<float>* A, std::size_t lda,
        const std::complex<float>* beta, std::complex<float>* C,
        std::size_t ldc ){
        return cblas_csyrk( order, uplo, trans, n, k, (float*) alpha,
            (float*) A, lda, (float*) beta, (float*) C, ldc );
    }
    void syrk( CBLAS_ORDER order, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans,
        std::size_t n, std::size_t k, const std::complex<double>* alpha,
        const std::complex<double>* A, std::size_t lda,
        const std::complex<double>* beta, std::complex<double>* C,
        std::size_t ldc ){
        return cblas_zsyrk( order, uplo, trans, n, k, (double*) alpha,
            (double*) A, lda, (double*) beta, (double*) C, ldc );
    }

} //namespace blas

#endif

ext: h

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


text:
/********************************\
* Matrix Multiply Test Benchmark *
*                                *
* by                             *
*   Elliott Forney               *
*   3.15.2010                    *
\********************************/

/*
 *  Libraries
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <unistd.h>
#include <getopt.h>

#include <cblas.h>

#include "matrix.h"
#include "benchmark.h"
#include "errcheck.h"

/*
 *  Macros
 */

// print command line usage
#define print_usage() fprintf(stdout, "Usage: %s [-m a_rows] [-n a_cols] [-p b_cols] [-s]\n", arg[0])

// acceptable maximum relative error tolerance
#define TOLERANCE 0.01

#define DEBUG 0

/*
 *  Global variables
 */

// flag for simple, one line output
bool simple_out = false;

/*
 *  Function prototypes
 */

// parse command line arguments
void parse_args(int narg, char **arg,
                unsigned *m, unsigned *n,
                unsigned *p);

/*
 *  Function bodies
 */

// setup network
int main(int narg, char **arg)
{
  // default matrix dimensions
  unsigned m = 6400;
  unsigned n = 0;
  unsigned p = 0;

  // parse command line arguments
  parse_args(narg, arg, &m, &n, &p);

  // if n or p == 0 then set to m
  if (n == 0)
    n = m;
  if (p == 0)
    p = m;

  if (DEBUG > 0)
    printf("m: %d\nn: %d\np: %d\n", m, n, p);

  // create matrix to hold result
  matrix result;
  matrix_init(&result, m, p);
  matrix_load_runif(result, 500, 600);

  // load test matrix a
  matrix a;
  matrix_init(&a, m, p);
//  matrix_load_testa(a, n);

  // load test matrix b
  matrix b;
  matrix_init(&b, m, n);
//  matrix_load_testb(b);
  matrix_load_runif(b, 0, 1);

  // load test matrix c
  matrix c;
  matrix_init(&c, n, p);
  matrix_load_testc(c);
  matrix_load_runif(c, 0, 1);

  // figure a with BLAS
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
              b.r, c.c, b.c, 1.0f,
              b.cpu_data, b.cstride, c.cpu_data, c.cstride,
              0.0f, a.cpu_data, a.cstride);

  // run warm up
  matrix_mult(result, b, c);

  // wait for all kernels to finish
  matrix_wait();

  // create a new benchmark timer
  benchmark ben;
  benchmark_init(&ben);

  // start timer
  benchmark_start_timer(&ben);

  // run multiplication
  matrix_mult(result, b, c);

  // wait for all kernels to finish
  matrix_wait();

  if (DEBUG > 5)
  {
    printf("a:\n");
    matrix_print_padded(a);
    printf("b:\n");
    matrix_print_padded(b);
    printf("c:\n");
    matrix_print_padded(c);
    printf("result:\n");
    matrix_print_padded(result);
  }

  // stop timer
  benchmark_stop_timer(&ben);

  // figure giga floating point operatins per second
  benchmark_add_flop(&ben, matmult_nflop(m,n,p));
  benchmark_add_byte(&ben, (m*n+n*p+m*p)*sizeof(float));
  double gbytes = benchmark_check_gbytes(ben);
  double gflops = benchmark_check_gflops(ben);
  double time   = benchmark_check_timer(ben);

  //
  float rmse = matrix_rmse(a, result);
  bool passed = true;
  if (rmse > TOLERANCE)
    passed = false;

  // if simple output requested
  if (simple_out)
    // simply print time, gflops and error
    printf("%f %f\n", time*1e3f, gflops);

  // if full output requested
  else
  {
    if (passed)
      printf("Test Passed!\n=======\n");
    else
      printf("Test Failed!\n=======\n");

    // print time and gflops
    printf("Time: %f  GByteS: %f  GFlopS: %f  RMSE: %f\n",
           time, gbytes, gflops, rmse);
  }

  // clean up
  matrix_dest(&a);
  matrix_dest(&b);
  matrix_dest(&c);
  matrix_dest(&result);
  benchmark_dest(&ben);

  // Come back soon now ya'hear!
  return 0;
}

// parse command line arguments
void parse_args(int narg, char **arg,
                unsigned *m, unsigned *n,
                unsigned *p)
{
  int opt; // getopt output

  // for each argument
  while ((opt = getopt(narg, arg, "sm:n:p:")) != -1)
  {
    if (opt == 's')
      simple_out = true;

    else if (opt == 'm')
      *m = (unsigned)atoi(arg[optind-1]);

    else if (opt == 'n')
      *n = (unsigned)atoi(arg[optind-1]);

    else if (opt == 'p')
      *p = (unsigned)atoi(arg[optind-1]);

    // print usage and quit on unknown option
    else
    {
      print_usage();
      exit(1);
    }
  }

  // assume last non dash arg is m
  if (optind < narg)
    *m = (unsigned)atoi(arg[optind]);
}

ext: c

max_stars_repo_name: idfah/badger

max_issues_repo_name: idfah/badger

max_forks_repo_name: idfah/badger


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#if 0
/*

Performs y = ( G^T G )^{-1} ( G^T K G ) ( G^T G )^{-1} x

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>

#include "common-driver-utils.h"

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/pcimpl.h"
     #include "petsc/private/kspimpl.h"
  #else
     #include "petsc-private/pcimpl.h"
     #include "petsc-private/kspimpl.h"
  #endif
#else
  #include "private/pcimpl.h"
  #include "private/kspimpl.h"
#endif

#include "pc_ScaledGtKG.h"

#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>

#define PCTYPE_SCGtKG "scgtkg"

/* private data */

typedef struct {
	Mat         B, Bt, F, C;   /* Bt \in [M x N], F \in [M x M] */
	Vec         X1,X2, Y1,Y2;  /* the scaling vectors */
	Vec         s,t,X;         /* s \in [M], t \in [N], X \in [M] */
	KSP         ksp_BBt;       /* The user MUST provide this operator */
	PetscTruth  BBt_has_cnst_nullspace;
	PetscTruth  monitor_activated;
} _PC_SC_GtKG;

typedef _PC_SC_GtKG* PC_SC_GtKG;

/* private prototypes */
PetscErrorCode BSSCR_PCApply_ScGtKG( PC pc, Vec x, Vec y );
PetscErrorCode BSSCR_PCApplyTranspose_ScGtKG( PC pc, Vec x, Vec y );
PetscErrorCode BSSCR_PCSetUp_GtKG( PC pc );

PetscErrorCode BSSCR_BSSCR_pc_error_ScGtKG( PC pc, const char func_name[] )
{
	const PCType type;
	PCGetType( pc, &type );

	if( strcmp(type,PCTYPE_SCGtKG)!=0 ) {
		printf("Error(%s): PC type (%s) should be scgtkg. \n",func_name, type );
		PetscFinalize();
		exit(0);
	}
	PetscFunctionReturn(0);
}

void BSSCR_BSSCR_get_number_nonzeros_AIJ_ScGtKG( Mat A, PetscInt *nnz )
{
	MatInfo info;

	MatGetInfo( A, MAT_GLOBAL_SUM, &info );
	*nnz = info.nz_used;

}

PetscErrorCode BSSCR_PCScGtKGBBtContainsConstantNullSpace( PC pc, PetscTruth *has_cnst_nullsp )
{
	PC_SC_GtKG    ctx = (PC_SC_GtKG)pc->data;
	PetscInt N;
	PetscScalar sum;
	PetscReal nrm;
	Vec l,r;
	Mat BBt,A;

	Stg_KSPGetOperators( ctx->ksp_BBt, &BBt, PETSC_NULL, PETSC_NULL );
	A = BBt;

	MatGetVecs( A, &r, &l ); // l = A r

	VecGetSize(r,&N);
	sum  = 1.0/N;
	VecSet(r,sum);

	/* scale */
	VecPointwiseMult( r, r, ctx->Y2 );
	/* {l} = [A] {r} */
	MatMult( A,r, l );
	/* scale */
	VecPointwiseMult( l, l, ctx->X2 );

	VecNorm(l,NORM_2,&nrm);
	if (nrm < 1.e-7) {
		*has_cnst_nullsp = PETSC_TRUE;
	}
	else {
		*has_cnst_nullsp = PETSC_FALSE;
	}

	Stg_VecDestroy(&l);
	Stg_VecDestroy(&r);

	PetscFunctionReturn(0);
}

/*
I should not modify setup called!!
This is handled via petsc.
*/
PetscErrorCode BSSCR_PCSetUp_ScGtKG( PC pc )
{
	PC_SC_GtKG    ctx = (PC_SC_GtKG)pc->data;

	if( ctx->F == PETSC_NULL ) {   Stg_SETERRQ( PETSC_ERR_SUP, "gtkg: F not set" );   }
	if( ctx->Bt == PETSC_NULL ) {  Stg_SETERRQ( PETSC_ERR_SUP, "gtkg: Bt not set" );  }

	if(ctx->ksp_BBt==PETSC_NULL) {
		BSSCR_PCScGtKGUseStandardBBtOperator( pc ) ;
	}

	BSSCR_PCScGtKGBBtContainsConstantNullSpace( pc, &ctx->BBt_has_cnst_nullspace );
	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		PetscPrintf( PETSC_COMM_WORLD, "\t* Detected prescence of constant nullspace in BBt-C\n" );
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCDestroy_ScGtKG( PC pc )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	if( ctx == PETSC_NULL ) {	PetscFunctionReturn(0); }

	if( ctx->ksp_BBt != PETSC_NULL ) {	Stg_KSPDestroy(&ctx->ksp_BBt );		}

	if( ctx->s != PETSC_NULL ) {		Stg_VecDestroy(&ctx->s );		}
	if( ctx->X != PETSC_NULL ) {		Stg_VecDestroy(&ctx->X );		}
	if( ctx->t != PETSC_NULL ) {		Stg_VecDestroy(&ctx->t );		}

	if( ctx->X1 != PETSC_NULL ) {		Stg_VecDestroy(&ctx->X1 );		}
	if( ctx->X2 != PETSC_NULL ) {		Stg_VecDestroy(&ctx->X2 );		}
	if( ctx->Y1 != PETSC_NULL ) {		Stg_VecDestroy(&ctx->Y1 );		}
	if( ctx->Y2 != PETSC_NULL ) {		Stg_VecDestroy(&ctx->Y2 );		}

	PetscFree( ctx );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCView_ScGtKG( PC pc, PetscViewer viewer )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	PetscViewerASCIIPushTab(viewer); //1

	PetscViewerASCIIPrintf( viewer, "gtkg-ksp \n" );
	PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
	PetscViewerASCIIPushTab(viewer);
		KSPView( ctx->ksp_BBt, viewer );
	PetscViewerASCIIPopTab(viewer);
	PetscViewerASCIIPrintf(viewer,"---------------------------------\n");

	PetscViewerASCIIPopTab(viewer); //1

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCRBSSCR_Lp_monitor_ScGtKG( KSP ksp, PetscInt index )
{
	PetscInt max_it;
	PetscReal rnorm;
	KSPConvergedReason reason;

	KSPGetIterationNumber( ksp, &max_it );
	KSPGetResidualNorm( ksp, &rnorm );
	KSPGetConvergedReason( ksp, &reason );
	if (reason >= 0) {
		PetscPrintf(((PetscObject)ksp)->comm,"\t<Lp(%d)>: Linear solve converged. its.=%.4d ; |r|=%5.5e ; Reason=%s\n",
				index, max_it, rnorm, KSPConvergedReasons[reason] );
	} else {
		PetscPrintf(((PetscObject)ksp)->comm,"\t<Lp(%d)>: Linear solve did not converge. its.=%.4d ; |r|=%5.5e ; Residual reduction=%2.2e ; Reason=%s\n",
				index, max_it, rnorm, (ksp->rnorm0/rnorm), KSPConvergedReasons[reason]);
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScBFBTSubKSPMonitor( KSP ksp, PetscInt index, PetscLogDouble time )
{
	PetscInt max_it;
	PetscReal rnorm;
	KSPConvergedReason reason;

	KSPGetIterationNumber( ksp, &max_it );
	KSPGetResidualNorm( ksp, &rnorm );
	KSPGetConvergedReason( ksp, &reason );
	PetscPrintf(((PetscObject)ksp)->comm,"    PCScBFBTSubKSP (%d): %D Residual norm; r0 %12.12e, r %12.12e: Reason %s: Time %5.5e \n",
			index, max_it, ksp->rnorm0, rnorm, KSPConvergedReasons[reason], time );

	PetscFunctionReturn(0);
}

/*
Performs y <- S*^{-1} x
  S*^{-1} = ( B' Bt' )^{-1} B' F' Bt' ( B' Bt' )^{-1}
where
  F'  = X1 F Y1
  Bt' = X1 Bt Y2
  B'  = X2 B Y1

Thus, S*^{-1} = [ X2 B Y1 X1 Bt Y2 ]^{-1} . [ X2 B Y1 . X1 F Y1 . X1 Bt Y2 ] . [ X2 B Y1 X1 Bt Y2 ]^{-1}
              = Y2^{-1} ksp_BBt X2^{-1} . [ B' F' Bt' ] . Y2^{-1} ksp_BBt X2^{-1}
              = Y2^{-1} ksp_BBt . [ B Y1 . X1 F Y1 . X1 Bt ] . ksp_BBt X2^{-1}
*/
PetscErrorCode BSSCR_PCApply_ScGtKG( PC pc, Vec x, Vec y )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	KSP ksp;
	Mat F, Bt;
	Vec s,t,X;
	PetscLogDouble t0,t1;

	ksp = ctx->ksp_BBt;
	F = ctx->F;
	Bt = ctx->Bt;
	s = ctx->s;
	t = ctx->t;
	X = ctx->X;

	/* Apply scaled Poisson operator */
	/* scale x */
	/* ========================================================
	     NOTE:
	       I THINK TO OMIT THESE AS WE WANT TO UNSCALE THE
	        PRECONDITIONER AS S IN THIS CASE IS NOT SCALED
	======================================================== */
//	VecPointwiseDivide( x, x, ctx->X2 ); /* x <- x/X2 */  /* NEED TO BE SURE */

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
	  BSSCR_VecRemoveConstNullspace( x, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolve( ksp, x, t ); /* t <- GtG_inv x */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,1,(t1-t0));
	}

	/* Apply Bt */
	MatMult( Bt, t, s ); /* s <- G t */
	VecPointwiseMult( s, s, ctx->X1 ); /* s <- s * X1 */

	/* Apply F */
	VecPointwiseMult( s, s, ctx->Y1 ); /* s <- s * Y1 */
	MatMult( F, s, X ); /* X <- K s */
	VecPointwiseMult( X, X, ctx->X1 ); /* X <- X * X1 */

	/* Apply B */
	VecPointwiseMult( X, X, ctx->Y1 ); /* X <- X * Y1 */
	MatMultTranspose( Bt, X, t ); /* t <- Gt X */

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		BSSCR_VecRemoveConstNullspace( t, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolve( ksp, t, y ); /* y <- GtG_inv t */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,2,(t1-t0));
	}

	VecPointwiseMult( y, y, ctx->Y2 ); /* y <- y/Y2 */

	/* undo modification made to x on entry */
//	VecPointwiseMult( x, x, ctx->X2 ); /* x <- x/X2 */ /* NEED TO BE SURE */

	PetscFunctionReturn(0);
}

/*
Performs y <- S*^{-1} x
  S*^{-1} = ( B' Bt' )^{-1} B' F' Bt' - C' ( B' Bt' )^{-1}
where
  F'  = X1 F Y1
  Bt' = X1 Bt Y2
  B'  = X2 B Y1
  C'  = X2 C Y2

Thus, S*^{-1} = [ X2 B Y1 X1 Bt Y2 ]^{-1} . [ X2 B Y1 . X1 F Y1 . X1 Bt Y2 - X2 C Y2 ] . [ X2 B Y1 X1 Bt Y2 ]^{-1}
              = Y2^{-1} ksp_BBt X2^{-1} . [ B' F' Bt' - C' ] . Y2^{-1} ksp_BBt X2^{-1}
              = Y2^{-1} ksp_BBt . [ B Y1 . X1 F Y1 . X1 Bt - C ] . ksp_BBt X2^{-1}
*/
PetscErrorCode BSSCR_BSSCR_PCApply_ScGtKG_C( PC pc, Vec x, Vec y )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	KSP ksp;
	Mat F, Bt,C;
	Vec s,t,X;
	PetscLogDouble t0,t1;

	ksp = ctx->ksp_BBt;
	F = ctx->F;
	Bt = ctx->Bt;
	C = ctx->C;
	s = ctx->s;
	t = ctx->t;
	X = ctx->X;

	/* Apply scaled Poisson operator */
	/* scale x */
	/* ========================================================
	     NOTE:
	       I THINK TO OMIT THESE AS WE WANT TO UNSCALE THE
	        PRECONDITIONER AS S IN THIS CASE IS NOT SCALED
	======================================================== */
//	VecPointwiseDivide( x, x, ctx->X2 ); /* x <- x/X2 */  /* NEED TO BE SURE */

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		BSSCR_VecRemoveConstNullspace( x, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolve( ksp, x, t ); /* t <- GtG_inv x */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,1,(t1-t0));
	}

	/* Apply Bt */
	MatMult( Bt, t, s ); /* s <- G t */
	VecPointwiseMult( s, s, ctx->X1 ); /* s <- s * X1 */

	/* Apply F */
	VecPointwiseMult( s, s, ctx->Y1 ); /* s <- s * Y1 */
	MatMult( F, s, X ); /* X <- K s */
	VecPointwiseMult( X, X, ctx->X1 ); /* X <- X * X1 */

	/* Apply B */
	VecPointwiseMult( X, X, ctx->Y1 ); /* X <- X * Y1 */
	MatMultTranspose( Bt, X, s ); /* s <- Gt X */

	/* s <- s - C t */
	VecScale( s, -1.0 );
	MatMultAdd( C, t, s, s );
	VecScale( s, -1.0 );

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		BSSCR_VecRemoveConstNullspace( s, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolve( ksp, s, y ); /* y <- GtG_inv s */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,2,(t1-t0));
	}

	VecPointwiseMult( y, y, ctx->Y2 ); /* y <- y/Y2 */

	/* undo modification made to x on entry */
//	VecPointwiseMult( x, x, ctx->X2 ); /* x <- x/X2 */ /* NEED TO BE SURE */

	PetscFunctionReturn(0);
}

/* Need to check this one if correct */
/*
S^{-1} = ( G^T G )^{-1} G^T K G ( G^T G )^{-1}
       = A C A
S^{-T} = A^T (A C)^T
       = A^T C^T A^T, but A = G^T G which is symmetric
       = A C^T A
       = A G^T ( G^T K )^T A
       = A G^T K^T G A

*/
PetscErrorCode BSSCR_PCApplyTranspose_ScGtKG( PC pc, Vec x, Vec y )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	KSP ksp;
	Mat F, Bt;
	Vec s,t,X;
	PetscLogDouble t0,t1;

	ksp = ctx->ksp_BBt;
	F = ctx->F;
	Bt = ctx->Bt;
	s = ctx->s;
	t = ctx->t;
	X = ctx->X;

	/* Apply scaled Poisson operator */
	/* scale x */
	/* ========================================================
	     NOTE:
	       I THINK TO OMIT THESE AS WE WANT TO UNSCALE THE
	        PRECONDITIONER AS S IN THIS CASE IS NOT SCALED
	======================================================== */
//	VecPointwiseDivide( x, x, ctx->X2 ); /* x <- x/X2 */  /* NEED TO BE SURE */

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		BSSCR_VecRemoveConstNullspace( x, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolveTranspose( ksp, x, t ); /* t <- GtG_inv x */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,1,(t1-t0));
	}

	/* Apply Bt */
	MatMult( Bt, t, s ); /* s <- G t */
	VecPointwiseMult( s, s, ctx->X1 ); /* s <- s * X1 */

	/* Apply F */
	VecPointwiseMult( s, s, ctx->Y1 ); /* s <- s * Y1 */
	MatMultTranspose( F, s, X ); /* X <- K s */
	VecPointwiseMult( X, X, ctx->X1 ); /* X <- X * X1 */

	/* Apply B */
	VecPointwiseMult( X, X, ctx->Y1 ); /* X <- X * Y1 */
	MatMultTranspose( Bt, X, t ); /* t <- Gt X */

	if( ctx->BBt_has_cnst_nullspace == PETSC_TRUE ) {
		BSSCR_VecRemoveConstNullspace( t, PETSC_NULL );
	}
	PetscGetTime(&t0);
	KSPSolveTranspose( ksp, t, y ); /* y <- GtG_inv t */
	PetscGetTime(&t1);
	if (ctx->monitor_activated) {
		BSSCR_PCScBFBTSubKSPMonitor(ksp,2,(t1-t0));
	}

	VecPointwiseMult( y, y, ctx->Y2 ); /* y <- y/Y2 */

	/* undo modification made to x on entry */
//	VecPointwiseMult( x, x, ctx->X2 ); /* x <- x/X2 */ /* NEED TO BE SURE */

	PetscFunctionReturn(0);
}

/*
Performs y <- S^{-1} x
S^{-1} = ( G^T Di G )^{-1} G^T Di K Di G ( G^T Di G )^{-1}
where Di = diag(M)^{-1}
*/

/*
Only the options related to GtKG should be set here.
*/
PetscErrorCode BSSCR_PCSetFromOptions_ScGtKG( PC pc )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	PetscTruth ivalue, flg;

	if(ctx->ksp_BBt!=PETSC_NULL) {
		PetscOptionsGetTruth( PETSC_NULL, "-pc_gtkg_monitor", &ivalue, &flg );
		BSSCR_PCScGtKGSetSubKSPMonitor( pc, ivalue );

	}

	PetscFunctionReturn(0);
}

/* ---- Exposed functions ---- */

PetscErrorCode BSSCR_PCCreate_ScGtKG( PC pc )
{
	PC_SC_GtKG      pc_data;
	PetscErrorCode  ierr;

	/* create memory for ctx */
	ierr = Stg_PetscNew( _PC_SC_GtKG,&pc_data);CHKERRQ(ierr);

	/* init ctx */
	pc_data->F   = PETSC_NULL;
	pc_data->Bt  = PETSC_NULL;
	pc_data->B   = PETSC_NULL;
	pc_data->BBt_has_cnst_nullspace = PETSC_FALSE;
	pc_data->ksp_BBt   = PETSC_NULL;
	pc_data->monitor_activated = PETSC_FALSE;

	pc_data->X1 = PETSC_NULL;
	pc_data->X2 = PETSC_NULL;
	pc_data->Y1 = PETSC_NULL;
	pc_data->Y2 = PETSC_NULL;

	pc_data->s = PETSC_NULL;
	pc_data->t = PETSC_NULL;
	pc_data->X = PETSC_NULL;

	/* set ctx onto pc */
	pc->data  = (void*)pc_data;

	ierr = PetscLogObjectMemory(pc,sizeof(_PC_SC_GtKG));CHKERRQ(ierr);

	/* define operations */
	pc->ops->setup   = BSSCR_PCSetUp_ScGtKG;
	pc->ops->view    = BSSCR_PCView_ScGtKG;
	pc->ops->destroy = BSSCR_PCDestroy_ScGtKG;
	pc->ops->setfromoptions = BSSCR_PCSetFromOptions_ScGtKG;

	pc->ops->apply          = BSSCR_PCApply_ScGtKG;
	pc->ops->applytranspose = BSSCR_PCApplyTranspose_ScGtKG;

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGGetScalings( PC pc, Vec *X1, Vec *X2, Vec *Y1, Vec *Y2 )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	if( X1 != PETSC_NULL ) {  *X1 = ctx->X1;  }
	if( X2 != PETSC_NULL ) {  *X2 = ctx->X2;  }

	if( Y1 != PETSC_NULL ) {  *Y1 = ctx->Y1;  }
	if( Y2 != PETSC_NULL ) {  *Y2 = ctx->Y2;  }

	PetscFunctionReturn(0);
}

/*
F & Bt must different to PETSC_NULL
B may be PETSC_NULL
C can be PETSC_NULL
*/
PetscErrorCode BSSCR_PCScGtKGSetOperators( PC pc, Mat F, Mat Bt, Mat B, Mat C )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	ctx->F  = F;
	ctx->Bt = Bt;
	ctx->B  = B;
	ctx->C  = C;

	if( C != PETSC_NULL ) {
		pc->ops->apply          = BSSCR_BSSCR_PCApply_ScGtKG_C;
		pc->ops->applytranspose = PETSC_NULL;
	}

	/* Create vectors */
	if( ctx->s == PETSC_NULL ) {  MatGetVecs( ctx->F, &ctx->s, PETSC_NULL );  }
	if( ctx->X == PETSC_NULL ) {  MatGetVecs( ctx->F, PETSC_NULL, &ctx->X );  }
	if( ctx->t == PETSC_NULL ) {  MatGetVecs( ctx->Bt, &ctx->t, PETSC_NULL ); }

	if( ctx->F == PETSC_NULL ) {   Stg_SETERRQ( PETSC_ERR_SUP, "gtkg: F not set" );   }
	if( ctx->Bt == PETSC_NULL ) {  Stg_SETERRQ( PETSC_ERR_SUP, "gtkg: Bt not set" );  }

	if( ctx->X1 == PETSC_NULL ) {  MatGetVecs( ctx->F, &ctx->X1, PETSC_NULL );  }
	if( ctx->Y1 == PETSC_NULL ) {  MatGetVecs( ctx->F, &ctx->Y1, PETSC_NULL );  }
	if( ctx->X2 == PETSC_NULL ) {  MatGetVecs( ctx->Bt, &ctx->X2, PETSC_NULL ); }
	if( ctx->Y2 == PETSC_NULL ) {  MatGetVecs( ctx->Bt, &ctx->Y2, PETSC_NULL ); }

	VecSet( ctx->X1, 1.0 );
	VecSet( ctx->Y1, 1.0 );
	VecSet( ctx->X2, 1.0 );
	VecSet( ctx->Y2, 1.0 );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGAttachNullSpace( PC pc )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	MatNullSpace nsp;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	/* Attach a null space */
	MatNullSpaceCreate( PETSC_COMM_WORLD, PETSC_TRUE, PETSC_NULL, PETSC_NULL, &nsp );
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR <6) )
	KSPSetNullSpace( ctx->ksp_BBt, nsp );
#else
    Mat A;
    KSPGetOperators(ctx->ksp_BBt,&A,NULL);//Note: DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.
    MatSetNullSpace( A, nsp);
#endif
	/*
	NOTE: This does NOT destroy the memory for nsp, it just decrements the nsp->refct, so that
	the next time MatNullSpaceDestroy() is called, the memory will be released. The next time this
	is called will be by KSPDestroy();
	*/
	MatNullSpaceDestroy( nsp );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGGetKSP( PC pc, KSP *ksp )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	if( ksp != PETSC_NULL ) {
		(*ksp) = ctx->ksp_BBt;
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGSetKSP( PC pc, KSP ksp )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	if( ctx->ksp_BBt != PETSC_NULL ) {
		Stg_KSPDestroy(&ctx->ksp_BBt);
	}

	PetscObjectReference( (PetscObject)ksp );
	ctx->ksp_BBt = ksp;

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGSetSubKSPMonitor( PC pc, PetscTruth flg )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	ctx->monitor_activated = flg;

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCScGtKGUseStandardScaling( PC pc )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	Mat K,G,D,C;
	Vec rG;
	PetscScalar rg2, rg, ra;
	PetscInt N;
	Vec rA, rC;
	Vec L1,L2, R1,R2;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	L1 = ctx->X1;
	L2 = ctx->X2;

	R1 = ctx->Y1;
	R2 = ctx->Y2;

	rA = L1;
	rC = L2;

	K = ctx->F;
	G = ctx->Bt;
	D = ctx->B;
	C = ctx->C;

	VecDuplicate( rA, &rG );

	/* Get magnitude of K */
	MatGetRowMax( K, rA, PETSC_NULL );

	VecSqrt( rA );
	VecReciprocal( rA );

	VecDot( rA,rA, &ra );
	VecGetSize( rA, &N );
	ra = PetscSqrtScalar( ra/N );

	/* Get magnitude of G */
	MatGetRowMax( G, rG, PETSC_NULL );

	VecDot( rG, rG, &rg2 );
	VecGetSize( rG, &N );
	rg = PetscSqrtScalar(rg2/N);
	//      printf("rg = %f \n", rg );

	VecSet( rC, 1.0/(rg*ra) );

	Stg_VecDestroy(&rG );

	VecCopy( L1, R1 );
	VecCopy( L2, R2 );

	PetscFunctionReturn(0);
}

/*
Builds
  B Y1 X1 Bt
and creates a ksp when C=0, otherwise it builds
  B Y1 X1 Bt - C
*/
PetscErrorCode BSSCR_PCScGtKGUseStandardBBtOperator( PC pc )
{
	PC_SC_GtKG ctx = (PC_SC_GtKG)pc->data;
	PetscReal  fill;
	Mat        diag_mat,C;
	Vec        diag;
	PetscInt   M,N, m,n;
	MPI_Comm   comm;
	PetscInt   nnz_I, nnz_G;
	MatType    mtype;
	const char *prefix;
	Mat BBt;
	KSP ksp;
	PetscTruth ivalue, flg, has_cnst_nullsp;

	BSSCR_BSSCR_pc_error_ScGtKG( pc, __func__ );

	/* Assemble BBt */
	MatGetSize( ctx->Bt, &M, &N );
	MatGetLocalSize( ctx->Bt, &m, &n );

	MatGetVecs( ctx->Bt, PETSC_NULL, &diag );

	/* Define diagonal matrix Y1 X1 */
	VecPointwiseMult( diag, ctx->Y1, ctx->X1 );

	PetscObjectGetComm( (PetscObject)ctx->F, &comm );
	MatCreate( comm, &diag_mat );
	MatSetSizes( diag_mat, m,m , M, M );
#if (((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=3)) || (PETSC_VERSION_MAJOR>3) )
        MatSetUp(diag_mat);
#endif
	MatGetType( ctx->Bt, &mtype );
	MatSetType( diag_mat, mtype );

	MatDiagonalSet( diag_mat, diag, INSERT_VALUES );

	/* Build operator B Y1 X1 Bt */
	BSSCR_BSSCR_get_number_nonzeros_AIJ_ScGtKG( diag_mat, &nnz_I );
	BSSCR_BSSCR_get_number_nonzeros_AIJ_ScGtKG( ctx->Bt, &nnz_G );
	/*
	Not sure the best way to estimate the fill factor.
	BBt is a laplacian on the pressure space.
	This might tell us something useful...
	*/
	fill = (PetscReal)(nnz_G)/(PetscReal)( nnz_I );
	MatPtAP( diag_mat, ctx->Bt, MAT_INITIAL_MATRIX, fill, &BBt );

	Stg_MatDestroy(&diag_mat );
	Stg_VecDestroy(&diag );

	C = ctx->C;
	if( C !=PETSC_NULL ) {
		MatAXPY( BBt, -1.0, C, DIFFERENT_NONZERO_PATTERN );
	}

	/* Build the solver */
	KSPCreate( ((PetscObject)pc)->comm, &ksp );

	Stg_KSPSetOperators( ksp, BBt, BBt, SAME_NONZERO_PATTERN );

	PCGetOptionsPrefix( pc,&prefix );
	KSPSetOptionsPrefix( ksp, prefix );
	KSPAppendOptionsPrefix( ksp, "pc_gtkg_" ); /* -pc_GtKG_ksp_type <type>, -ksp_GtKG_pc_type <type> */

	BSSCR_PCScGtKGSetKSP( pc, ksp );

	BSSCR_MatContainsConstNullSpace( BBt, NULL, &has_cnst_nullsp );
	if( has_cnst_nullsp == PETSC_TRUE ) {
		BSSCR_PCScGtKGAttachNullSpace( pc );
	}

	PetscOptionsGetTruth( PETSC_NULL, "-pc_gtkg_monitor", &ivalue, &flg );
	BSSCR_PCScGtKGSetSubKSPMonitor( pc, ivalue );

	Stg_KSPDestroy(&ksp);
	Stg_MatDestroy(&BBt);

	PetscFunctionReturn(0);
}

#endif


text:
#pragma once

#include "SimpleEvaluator.h"
#include "BooleanDAG.h"
#include "Interface.h"
#include "SuggestionGenerator.h"
#include <map>
#include <set>

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif

class SimpleSuggestionGenerator: public SuggestionGenerator {
    vector<int> nodesToSuggest;
    Interface* interf;
    BooleanDAG* dag;
    map<string, int>& ctrls;
    SimpleEvaluator* seval;
public:
    SimpleSuggestionGenerator(BooleanDAG* _dag, Interface* _interface, map<string, int>& _ctrls): dag(_dag), interf(_interface), ctrls(_ctrls) {
        seval = new SimpleEvaluator(*dag, ctrls);
        seval->setInputs(interf);
        for (auto it = interf->varsMapping.begin(); it != interf->varsMapping.end(); it++) {
            int nodeid = it->first;
            bool_node* n = (*dag)[nodeid];
            if (n->type == bool_node::LT) {
                nodesToSuggest.push_back(nodeid);
            }
        }
    }

    pair<int, int> getSuggestion(const gsl_vector* state) {
        return make_pair(-1, 0);
    }
    pair<int, int> getUnsatSuggestion(const gsl_vector* state) {
        return make_pair(-1, 0);
    }
    virtual IClause* getConflictClause(int level, LocalState * state) {
        return NULL;
    }

    virtual SClause* getUnsatClause(const gsl_vector* state) {
        return NULL;
    }

    virtual SClause* getUnsatClause(const gsl_vector* state, const gsl_vector* prevState) {
        return NULL;
    }
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/* Copied from the article:
http://amitsaha.github.io/site/notes/articles/c_scientific/article.html
*/

/*Listing-1: gsl_vector.c*/

/* Simple demo of the vector support in GSL
 *  Also uses the random number generation feature
 */

#include <stdio.h>
#include <gsl/gsl_vector.h> /*For Vectors*/
#include <gsl/gsl_rng.h> /* For Random numbers*/

int main ()
{
  int i,n;
  /* Setup the Random number generator*/
  const gsl_rng_type * T;
  gsl_rng * r;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  n = 10;

  /* Allocate the vector of the specified size*/
  gsl_vector * v = gsl_vector_alloc (n);

  /* Set the elements to a uniform random number in [0,1]*/
  for (i = 0; i < n; i++)
    {
      gsl_vector_set (v, i, gsl_rng_uniform (r));
    }

  /* Print the vector*/
  for (i = 0; i < n; i++)
    {
      printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
    }

  gsl_vector_free (v);

  return 0;
}