text:
/**
 *
 * @file qwrapper_dlange.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlange(Quark *quark, Quark_Task_Flags *task_flags,
                       int norm, int M, int N,
                       const double *A, int LDA, int szeA,
                       int szeW, double *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_dlange_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(double)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(double)*szeW,                 NULL,          SCRATCH,
        sizeof(double),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlange_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, int M, int N,
                          const double *A, int LDA, int szeA,
                          int szeW, double *result,
                          double *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_dlange_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(double)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_dlange_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(double)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT,
            sizeof(double)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlange_quark = PCORE_dlange_quark
#define CORE_dlange_quark PCORE_dlange_quark
#endif
void CORE_dlange_quark(Quark *quark)
{
    double *normA;
    int norm;
    int M;
    int N;
    double *A;
    int LDA;
    double *work;

    quark_unpack_args_7(quark, norm, M, N, A, LDA, work, normA);
    *normA = LAPACKE_dlange_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlange_f1_quark = PCORE_dlange_f1_quark
#define CORE_dlange_f1_quark PCORE_dlange_f1_quark
#endif
void CORE_dlange_f1_quark(Quark *quark)
{
    double *normA;
    int norm;
    int M;
    int N;
    double *A;
    int LDA;
    double *work;
    double *fake;

    quark_unpack_args_8(quark, norm, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_dlange_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        M, N, A, LDA, work);
}


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_laguerre.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_hyperg.h>

#include "error.h"
#include "hyperg.h"

#define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON)

#define SERIES_EVAL_OK(a,b,x) ((fabs(a) < 5 && b < 5 && x < 2.0) || (fabs(a) <  10 && b < 10 && x < 1.0))

#define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x))

/* Log[U(a,2a,x)]
 * [Abramowitz+stegun, 13.6.21]
 * Assumes x > 0, a > 1/2.
 */
static
int
hyperg_lnU_beq2a(const double a, const double x, gsl_sf_result * result)
{
  const double lx = log(x);
  const double nu = a - 0.5;
  const double lnpre = 0.5*(x - M_LNPI) - nu*lx;
  gsl_sf_result lnK;
  gsl_sf_bessel_lnKnu_e(nu, 0.5*x, &lnK);
  result->val  = lnpre + lnK.val;
  result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + 0.5*M_LNPI + fabs(nu*lx));
  result->err += lnK.err;
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return GSL_SUCCESS;
}

/* Evaluate u_{N+1}/u_N by Steed's continued fraction method.
 *
 * u_N := Gamma[a+N]/Gamma[a] U(a + N, b, x)
 *
 * u_{N+1}/u_N = (a+N) U(a+N+1,b,x)/U(a+N,b,x)
 */
static
int
hyperg_U_CF1(const double a, const double b, const int N, const double x,
             double * result, int * count)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 20000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = -(a + N);
  double b1 =  (b - 2.0*a - x - 2.0*(N+1));
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = -(a + N + n - b)*(a + N + n - 1.0);
    bn =  (b - 2.0*a - x - 2.0*(N+n));
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;

    if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  }

  *result = fn;
  *count  = n;

  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Large x asymptotic for  x^a U(a,b,x)
 * Based on SLATEC D9CHU() [W. Fullerton]
 *
 * Uses a rational approximation due to Luke.
 * See [Luke, Algorithms for the Computation of Special Functions, p. 252]
 *     [Luke, Utilitas Math. (1977)]
 *
 * z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
 *
 * This assumes that a is not a negative integer and
 * that 1+a-b is not a negative integer. If one of them
 * is, then the 2F0 actually terminates, the above
 * relation is an equality, and the sum should be
 * evaluated directly [see below].
 */
static
int
d9chu(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double EPS   = 8.0 * GSL_DBL_EPSILON;  /* EPS = 4.0D0*D1MACH(4)   */
  const int maxiter = 500;
  double aa[4], bb[4];
  int i;

  double bp = 1.0 + a - b;
  double ab = a*bp;
  double ct2 = 2.0 * (x - ab);
  double sab = a + bp;

  double ct3 = sab + 1.0 + ab;
  double anbn = ct3 + sab + 3.0;
  double ct1 = 1.0 + 2.0*x/anbn;

  bb[0] = 1.0;
  aa[0] = 1.0;

  bb[1] = 1.0 + 2.0*x/ct3;
  aa[1] = 1.0 + ct2/ct3;

  bb[2] = 1.0 + 6.0*ct1*x/ct3;
  aa[2] = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3;

  for(i=4; i<maxiter; i++) {
    int j;
    double c2;
    double d1z;
    double g1, g2, g3;
    double x2i1 = 2*i - 3;
    ct1   = x2i1/(x2i1-2.0);
    anbn += x2i1 + sab;
    ct2   = (x2i1 - 1.0)/anbn;
    c2    = x2i1*ct2 - 1.0;
    d1z   = 2.0*x2i1*x/anbn;

    ct3 = sab*ct2;
    g1  = d1z + ct1*(c2+ct3);
    g2  = d1z - c2;
    g3  = ct1*(1.0 - ct3 - 2.0*ct2);

    bb[3] = g1*bb[2] + g2*bb[1] + g3*bb[0];
    aa[3] = g1*aa[2] + g2*aa[1] + g3*aa[0];

    if(fabs(aa[3]*bb[0]-aa[0]*bb[3]) < EPS*fabs(bb[3]*bb[0])) break;

    for(j=0; j<3; j++) {
      aa[j] = aa[j+1];
      bb[j] = bb[j+1];
    }
  }

  result->val = aa[3]/bb[3];
  result->err = 8.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(i == maxiter) {
    GSL_ERROR ("error", GSL_EMAXITER);
  }
  else {
    return GSL_SUCCESS;
  }
}

/* Evaluate asymptotic for z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
 * We check for termination of the 2F0 as a special case.
 * Assumes x > 0.
 * Also assumes a,b are not too large compared to x.
 */
static
int
hyperg_zaU_asymp(const double a, const double b, const double x, gsl_sf_result *result)
{
  const double ap = a;
  const double bp = 1.0 + a - b;
  const double rintap = floor(ap + 0.5);
  const double rintbp = floor(bp + 0.5);
  const int ap_neg_int = ( ap < 0.0 && fabs(ap - rintap) < INT_THRESHOLD );
  const int bp_neg_int = ( bp < 0.0 && fabs(bp - rintbp) < INT_THRESHOLD );

  if(ap_neg_int || bp_neg_int) {
    /* Evaluate 2F0 polynomial.
     */
    double mxi = -1.0/x;
    double nmax = -(int)(GSL_MIN(ap,bp) - 0.1);
    double tn  = 1.0;
    double sum = 1.0;
    double n   = 1.0;
    double sum_err = 0.0;
    while(n <= nmax) {
      double apn = (ap+n-1.0);
      double bpn = (bp+n-1.0);
      tn  *= ((apn/n)*mxi)*bpn;
      sum += tn;
      sum_err += 2.0 * GSL_DBL_EPSILON * fabs(tn);
      n += 1.0;
    }
    result->val  = sum;
    result->err  = sum_err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(nmax)+1.0) * fabs(sum);
    return GSL_SUCCESS;
  }
  else {
    return d9chu(a,b,x,result);
  }
}

/* Evaluate finite sum which appears below.
 */
static
int
hyperg_U_finite_sum(int N, double a, double b, double x, double xeps,
                    gsl_sf_result * result)
{
  int i;
  double sum_val;
  double sum_err;

  if(N <= 0) {
    double t_val = 1.0;
    double t_err = 0.0;
    gsl_sf_result poch;
    int stat_poch;

    sum_val = 1.0;
    sum_err = 0.0;
    for(i=1; i<= -N; i++) {
      const double xi1  = i - 1;
      const double mult = (a+xi1)*x/((b+xi1)*(xi1+1.0));
      t_val *= mult;
      t_err += fabs(mult) * t_err + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
      sum_val += t_val;
      sum_err += t_err;
    }

    stat_poch = gsl_sf_poch_e(1.0+a-b, -a, &poch);

    result->val  = sum_val * poch.val;
    result->err  = fabs(sum_val) * poch.err + sum_err * fabs(poch.val);
    result->err += fabs(poch.val) * (fabs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
    return stat_poch;
  }
  else {
    const int M = N - 2;
    if(M < 0) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else {
      gsl_sf_result gbm1;
      gsl_sf_result gamr;
      int stat_gbm1;
      int stat_gamr;
      double t_val = 1.0;
      double t_err = 0.0;

      sum_val = 1.0;
      sum_err = 0.0;
      for(i=1; i<=M; i++) {
        const double mult = (a-b+i)*x/((1.0-b+i)*i);
        t_val *= mult;
        t_err += t_err * fabs(mult) + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
        sum_val += t_val;
        sum_err += t_err;
      }

      stat_gbm1 = gsl_sf_gamma_e(b-1.0, &gbm1);
      stat_gamr = gsl_sf_gammainv_e(a,  &gamr);

      if(stat_gbm1 == GSL_SUCCESS) {
        gsl_sf_result powx1N;
        int stat_p = gsl_sf_pow_int_e(x, 1-N, &powx1N);
        double pe_val = powx1N.val * xeps;
        double pe_err = powx1N.err * fabs(xeps) + 2.0 * GSL_DBL_EPSILON * fabs(pe_val);
        double coeff_val = gbm1.val * gamr.val * pe_val;
        double coeff_err = gbm1.err * fabs(gamr.val * pe_val)
                         + gamr.err * fabs(gbm1.val * pe_val)
                         + fabs(gbm1.val * gamr.val) * pe_err
                         + 2.0 * GSL_DBL_EPSILON * fabs(coeff_val);

        result->val  = sum_val * coeff_val;
        result->err  = fabs(sum_val) * coeff_err + sum_err * fabs(coeff_val);
        result->err += 2.0 * GSL_DBL_EPSILON * (M+2.0) * fabs(result->val);
        result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
        return stat_p;
      }
      else {
        result->val = 0.0;
        result->err = 0.0;
        return stat_gbm1;
      }
    }
  }
}

/* Based on SLATEC DCHU() [W. Fullerton]
 * Assumes x > 0.
 * This is just a series summation method, and
 * it is not good for large a.
 *
 * I patched up the window for 1+a-b near zero. [GJ]
 */
static
int
hyperg_U_series(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */
  const double SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON;

  if(fabs(1.0 + a - b) < SQRT_EPS) {
    /* Original Comment: ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X
     */
    /* We can however do the following:
     * U(a,b,x) = U(a,a+1,x) when 1+a-b=0
     * and U(a,a+1,x) = x^(-a).
     */
    double lnr = -a * log(x);
    int stat_e =  gsl_sf_exp_e(lnr, result);
    result->err += 2.0 * SQRT_EPS * fabs(result->val);
    return stat_e;
  }
  else {
    double aintb = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) );
    double beps  = b - aintb;
    int N = aintb;

    double lnx  = log(x);
    double xeps = exp(-beps*lnx);

    /* Evaluate finite sum.
     */
    gsl_sf_result sum;
    int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum);

    /* Evaluate infinite sum. */

    int istrt = ( N < 1 ? 1-N : 0 );
    double xi = istrt;

    gsl_sf_result gamr;
    gsl_sf_result powx;
    int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr);
    int stat_powx = gsl_sf_pow_int_e(x, istrt, &powx);
    double sarg   = beps*M_PI;
    double sfact  = ( sarg != 0.0 ? sarg/sin(sarg) : 1.0 );
    double factor_val = sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * gamr.val * powx.val;
    double factor_err = fabs(gamr.val) * powx.err + fabs(powx.val) * gamr.err
                      + 2.0 * GSL_DBL_EPSILON * fabs(factor_val);

    gsl_sf_result pochai;
    gsl_sf_result gamri1;
    gsl_sf_result gamrni;
    int stat_pochai = gsl_sf_poch_e(a, xi, &pochai);
    int stat_gamri1 = gsl_sf_gammainv_e(xi + 1.0, &gamri1);
    int stat_gamrni = gsl_sf_gammainv_e(aintb + xi, &gamrni);
    int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);
    int stat_gamall = GSL_ERROR_SELECT_4(stat_sum, stat_gam123, stat_pochai, stat_powx);

    gsl_sf_result pochaxibeps;
    gsl_sf_result gamrxi1beps;
    int stat_pochaxibeps = gsl_sf_poch_e(a, xi-beps, &pochaxibeps);
    int stat_gamrxi1beps = gsl_sf_gammainv_e(xi + 1.0 - beps, &gamrxi1beps);

    int stat_all = GSL_ERROR_SELECT_3(stat_gamall, stat_pochaxibeps, stat_gamrxi1beps);

    double b0_val = factor_val * pochaxibeps.val * gamrni.val * gamrxi1beps.val;
    double b0_err =  fabs(factor_val * pochaxibeps.val * gamrni.val) * gamrxi1beps.err
                   + fabs(factor_val * pochaxibeps.val * gamrxi1beps.val) * gamrni.err
                   + fabs(factor_val * gamrni.val * gamrxi1beps.val) * pochaxibeps.err
                   + fabs(pochaxibeps.val * gamrni.val * gamrxi1beps.val) * factor_err
                   + 2.0 * GSL_DBL_EPSILON * fabs(b0_val);

    if(fabs(xeps-1.0) < 0.5) {
      /*
       C  X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE
       C  CAREFUL IN EVALUATING THE DIFFERENCES.
       */
      int i;
      gsl_sf_result pch1ai;
      gsl_sf_result pch1i;
      gsl_sf_result poch1bxibeps;
      int stat_pch1ai = gsl_sf_pochrel_e(a + xi, -beps, &pch1ai);
      int stat_pch1i  = gsl_sf_pochrel_e(xi + 1.0 - beps, beps, &pch1i);
      int stat_poch1bxibeps = gsl_sf_pochrel_e(b+xi, -beps, &poch1bxibeps);
      double c0_t1_val = beps*pch1ai.val*pch1i.val;
      double c0_t1_err = fabs(beps) * fabs(pch1ai.val) * pch1i.err
                       + fabs(beps) * fabs(pch1i.val)  * pch1ai.err
                       + 2.0 * GSL_DBL_EPSILON * fabs(c0_t1_val);
      double c0_t2_val = -poch1bxibeps.val + pch1ai.val - pch1i.val + c0_t1_val;
      double c0_t2_err =  poch1bxibeps.err + pch1ai.err + pch1i.err + c0_t1_err
                       + 2.0 * GSL_DBL_EPSILON * fabs(c0_t2_val);
      double c0_val = factor_val * pochai.val * gamrni.val * gamri1.val * c0_t2_val;
      double c0_err =  fabs(factor_val * pochai.val * gamrni.val * gamri1.val) * c0_t2_err
                     + fabs(factor_val * pochai.val * gamrni.val * c0_t2_val) * gamri1.err
                     + fabs(factor_val * pochai.val * gamri1.val * c0_t2_val) * gamrni.err
                     + fabs(factor_val * gamrni.val * gamri1.val * c0_t2_val) * pochai.err
                     + fabs(pochai.val * gamrni.val * gamri1.val * c0_t2_val) * factor_err
                     + 2.0 * GSL_DBL_EPSILON * fabs(c0_val);
      /*
       C  XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS)
       */
      gsl_sf_result dexprl;
      int stat_dexprl = gsl_sf_exprel_e(-beps*lnx, &dexprl);
      double xeps1_val = lnx * dexprl.val;
      double xeps1_err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(beps*lnx)) * fabs(dexprl.val)
                       + fabs(lnx) * dexprl.err
                       + 2.0 * GSL_DBL_EPSILON * fabs(xeps1_val);
      double dchu_val = sum.val + c0_val + xeps1_val*b0_val;
      double dchu_err = sum.err + c0_err
                      + fabs(xeps1_val)*b0_err + xeps1_err * fabs(b0_val)
                      + fabs(b0_val*lnx)*dexprl.err
                      + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(c0_val) + fabs(xeps1_val*b0_val));
      double xn = N;
      double t_val;
      double t_err;

      stat_all = GSL_ERROR_SELECT_5(stat_all, stat_dexprl, stat_poch1bxibeps, stat_pch1i, stat_pch1ai);

      for(i=1; i<2000; i++) {
        const double xi  = istrt + i;
        const double xi1 = istrt + i - 1;
        const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps);
        const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps));
        const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi);
        const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps));
        b0_val *= b0_multiplier;
        b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
        c0_val  = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val;
        c0_err  =  fabs(c0_multiplier_1) * c0_err
                 + fabs(c0_multiplier_2) * b0_err
                 + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON
                 + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON;
        t_val  = c0_val + xeps1_val*b0_val;
        t_err  = c0_err + fabs(xeps1_val)*b0_err;
        t_err += fabs(b0_val*lnx) * dexprl.err;
        t_err += fabs(b0_val)*xeps1_err;
        dchu_val += t_val;
        dchu_err += t_err;
        if(fabs(t_val) < EPS*fabs(dchu_val)) break;
      }

      result->val  = dchu_val;
      result->err  = 2.0 * dchu_err;
      result->err += 2.0 * fabs(t_val);
      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {
        GSL_ERROR ("error", GSL_EMAXITER);
      }
      else {
        return stat_all;
      }
    }
    else {
      /*
       C  X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE
       C  STRAIGHTFORWARD FORMULATION IS STABLE.
       */
      int i;
      double dchu_val;
      double dchu_err;
      double t_val;
      double t_err;
      gsl_sf_result dgamrbxi;
      int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);
      double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps;
      double a0_err =  fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err
                     + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err
                     + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err
                     + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err
                     + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);
      stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi);

      b0_val = xeps * b0_val / beps;
      b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);
      dchu_val = sum.val + a0_val - b0_val;
      dchu_err = sum.err + a0_err + b0_err
        + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));

      for(i=1; i<2000; i++) {
        double xi = istrt + i;
        double xi1 = istrt + i - 1;
        double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);
        double b0_multiplier = (a+xi1-beps)*x/((aintb+xi1)*(xi-beps));
        a0_val *= a0_multiplier;
        a0_err += fabs(a0_multiplier) * a0_err;
        b0_val *= b0_multiplier;
        b0_err += fabs(b0_multiplier) * b0_err;
        t_val = a0_val - b0_val;
        t_err = a0_err + b0_err;
        dchu_val += t_val;
        dchu_err += t_err;
        if(fabs(t_val) < EPS*fabs(dchu_val)) break;
      }

      result->val  = dchu_val;
      result->err  = 2.0 * dchu_err;
      result->err += 2.0 * fabs(t_val);
      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {
        GSL_ERROR ("error", GSL_EMAXITER);
      }
      else {
        return stat_all;
      }
    }
  }
}

/* Assumes b > 0 and x > 0.
 */
static
int
hyperg_U_small_ab(const double a, const double b, const double x, gsl_sf_result * result)
{
  if(a == -1.0) {
    /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x
     */
    result->val  = -b + x;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(a == 0.0) {
    /* U(0,c+1,x) = Laguerre[c,0,x] = 1
     */
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ASYMP_EVAL_OK(a,b,x)) {
    double p = pow(x, -a);
    gsl_sf_result asymp;
    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
    result->val  = asymp.val * p;
    result->err  = asymp.err * p;
    result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_asymp;
  }
  else {
    return hyperg_U_series(a, b, x, result);
  }
}

/* Assumes b > 0 and x > 0.
 */
static
int
hyperg_U_small_a_bgt0(const double a, const double b, const double x,
                      gsl_sf_result * result,
                      double * ln_multiplier
                      )
{
  if(a == 0.0) {
    result->val = 1.0;
    result->err = 1.0;
    *ln_multiplier = 0.0;
    return GSL_SUCCESS;
  }
  else if(   (b > 5000.0 && x < 0.90 * fabs(b))
          || (b >  500.0 && x < 0.50 * fabs(b))
    ) {
    int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier);
    if(stat == GSL_EOVRFLW)
      return GSL_SUCCESS;
    else
      return stat;
  }
  else if(b > 15.0) {
    /* Recurse up from b near 1.
     */
    double eps = b - floor(b);
    double b0  = 1.0 + eps;
    gsl_sf_result r_Ubm1;
    gsl_sf_result r_Ub;
    int stat_0 = hyperg_U_small_ab(a, b0,     x, &r_Ubm1);
    int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub);
    double Ubm1 = r_Ubm1.val;
    double Ub   = r_Ub.val;
    double Ubp1;
    double bp;

    for(bp = b0+1.0; bp<b-0.1; bp += 1.0) {
      Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
      Ubm1 = Ub;
      Ub   = Ubp1;
    }
    result->val  = Ub;
    result->err  = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub);
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub);
    *ln_multiplier = 0.0;
    return GSL_ERROR_SELECT_2(stat_0, stat_1);
  }
  else {
    *ln_multiplier = 0.0;
    return hyperg_U_small_ab(a, b, x, result);
  }
}

/* We use this to keep track of large
 * dynamic ranges in the recursions.
 * This can be important because sometimes
 * we want to calculate a very large and
 * a very small number and the answer is
 * the product, of order 1. This happens,
 * for instance, when we apply a Kummer
 * transform to make b positive and
 * both x and b are large.
 */
#define RESCALE_2(u0,u1,factor,count)      \
do {                                       \
  double au0 = fabs(u0);                   \
  if(au0 > factor) {                       \
    u0 /= factor;                          \
    u1 /= factor;                          \
    count++;                               \
  }                                        \
  else if(au0 < 1.0/factor) {              \
    u0 *= factor;                          \
    u1 *= factor;                          \
    count--;                               \
  }                                        \
} while (0)

/* Specialization to b >= 1, for integer parameters.
 * Assumes x > 0.
 */
static
int
hyperg_U_int_bge1(const int a, const int b, const double x,
                  gsl_sf_result_e10 * result)
{
  if(a == 0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a == -1) {
    result->val  = -b + x;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->e10  = 0;
    return GSL_SUCCESS;
  }
  else if(b == a + 1) {
    /* U(a,a+1,x) = x^(-a)
     */
    return gsl_sf_exp_e10_e(-a*log(x), result);
  }
  else if(ASYMP_EVAL_OK(a,b,x)) {
    const double ln_pre_val = -a*log(x);
    const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
    gsl_sf_result asymp;
    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
    int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
                                              asymp.val, asymp.err,
                                              result);
    return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
  }
  else if(SERIES_EVAL_OK(a,b,x)) {
    gsl_sf_result ser;
    const int stat_ser = hyperg_U_series(a, b, x, &ser);
    result->val = ser.val;
    result->err = ser.err;
    result->e10 = 0;
    return stat_ser;
  }
  else if(a < 0) {
    /* Recurse backward from a = -1,0.
     */
    int scale_count = 0;
    const double scale_factor = GSL_SQRT_DBL_MAX;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double lnscale;
    double Uap1 = 1.0;     /* U(0,b,x)  */
    double Ua   = -b + x;  /* U(-1,b,x) */
    double Uam1;
    int ap;

    for(ap=-1; ap>a; ap--) {
      Uam1 = ap*(b-ap-1.0)*Uap1 + (x+2.0*ap-b)*Ua;
      Uap1 = Ua;
      Ua   = Uam1;
      RESCALE_2(Ua,Uap1,scale_factor,scale_count);
    }

    lnscale = log(scale_factor);
    lnm.val = scale_count*lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val);
    y.val = Ua;
    y.err = 4.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(Ua);
    return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  }
  else if(b >= 2.0*a + x) {
    /* Recurse forward from a = 0,1.
     */
    int scale_count = 0;
    const double scale_factor = GSL_SQRT_DBL_MAX;
    gsl_sf_result r_Ua;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double lnscale;
    double lm;
    int stat_1 = hyperg_U_small_a_bgt0(1.0, b, x, &r_Ua, &lm);  /* U(1,b,x) */
    int stat_e;
    double Uam1 = 1.0;  /* U(0,b,x) */
    double Ua   = r_Ua.val;
    double Uap1;
    int ap;

    Uam1 *= exp(-lm);

    for(ap=1; ap<a; ap++) {
      Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
      Uam1 = Ua;
      Ua   = Uap1;
      RESCALE_2(Ua,Uam1,scale_factor,scale_count);
    }

    lnscale = log(scale_factor);
    lnm.val = lm + scale_count * lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale));
    y.val  = Ua;
    y.err  = fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
    y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ua);
    stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    return GSL_ERROR_SELECT_2(stat_e, stat_1);
  }
  else {
    if(b <= x) {
      /* Recurse backward either to the b=a+1 line
       * or to a=0, whichever we hit.
       */
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count = 0;
      int stat_CF1;
      double ru;
      int CF1_count;
      int a_target;
      double lnU_target;
      double Ua;
      double Uap1;
      double Uam1;
      int ap;

      if(b < a + 1) {
        a_target = b-1;
        lnU_target = -a_target*log(x);
      }
      else {
        a_target = 0;
        lnU_target = 0.0;
      }

      stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);

      Ua   = 1.0;
      Uap1 = ru/a * Ua;
      for(ap=a; ap>a_target; ap--) {
        Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
        Uap1 = Ua;
        Ua   = Uam1;
        RESCALE_2(Ua,Uap1,scale_factor,scale_count);
      }

      if(Ua == 0.0) {
        result->val = 0.0;
        result->err = 0.0;
        result->e10 = 0;
        GSL_ERROR ("error", GSL_EZERODIV);
      }
      else {
        double lnscl = -scale_count*log(scale_factor);
        double lnpre_val = lnU_target + lnscl;
        double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl));
        double oUa_err   = 2.0 * (fabs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua);
        int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err,
                                                  1.0/Ua, oUa_err,
                                                  result);
        return GSL_ERROR_SELECT_2(stat_e, stat_CF1);
      }
    }
    else {
      /* Recurse backward to near the b=2a+x line, then
       * determine normalization by either direct evaluation
       * or by a forward recursion. The direct evaluation
       * is needed when x is small (which is precisely
       * when it is easy to do).
       */
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count_for = 0;
      int scale_count_bck = 0;
      int a0 = 1;
      int a1 = a0 + ceil(0.5*(b-x) - a0);
      double Ua1_bck_val;
      double Ua1_bck_err;
      double Ua1_for_val;
      double Ua1_for_err;
      int stat_for;
      int stat_bck;
      gsl_sf_result lm_for;

      {
        /* Recurse back to determine U(a1,b), sans normalization.
         */
        double ru;
        int CF1_count;
        int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
        double Ua   = 1.0;
        double Uap1 = ru/a * Ua;
        double Uam1;
        int ap;
        for(ap=a; ap>a1; ap--) {
          Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
          Uap1 = Ua;
          Ua   = Uam1;
          RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
        }
        Ua1_bck_val = Ua;
        Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs(a1-a)+CF1_count+1.0) * fabs(Ua);
        stat_bck = stat_CF1;
      }

      if(b == 2*a1 && a1 > 1) {
        /* This can happen when x is small, which is
         * precisely when we need to be careful with
         * this evaluation.
         */
        hyperg_lnU_beq2a((double)a1, x, &lm_for);
        Ua1_for_val = 1.0;
        Ua1_for_err = 0.0;
        stat_for = GSL_SUCCESS;
      }
      else if(b == 2*a1 - 1 && a1 > 1) {
        /* Similar to the above. Happens when x is small.
         * Use
         *   U(a,2a-1) = (x U(a,2a) - U(a-1,2(a-1))) / (2a - 2)
         */
        gsl_sf_result lnU00, lnU12;
        gsl_sf_result U00, U12;
        hyperg_lnU_beq2a(a1-1.0, x, &lnU00);
        hyperg_lnU_beq2a(a1,     x, &lnU12);
        if(lnU00.val > lnU12.val) {
          lm_for.val = lnU00.val;
          lm_for.err = lnU00.err;
          U00.val = 1.0;
          U00.err = 0.0;
          gsl_sf_exp_err_e(lnU12.val - lm_for.val, lnU12.err + lm_for.err, &U12);
        }
        else {
          lm_for.val = lnU12.val;
          lm_for.err = lnU12.err;
          U12.val = 1.0;
          U12.err = 0.0;
          gsl_sf_exp_err_e(lnU00.val - lm_for.val, lnU00.err + lm_for.err, &U00);
        }
        Ua1_for_val  = (x * U12.val - U00.val) / (2.0*a1 - 2.0);
        Ua1_for_err  = (fabs(x)*U12.err + U00.err) / fabs(2.0*a1 - 2.0);
        Ua1_for_err += 2.0 * GSL_DBL_EPSILON * fabs(Ua1_for_val);
        stat_for = GSL_SUCCESS;
      }
      else {
        /* Recurse forward to determine U(a1,b) with
         * absolute normalization.
         */
        gsl_sf_result r_Ua;
        double Uam1 = 1.0;  /* U(a0-1,b,x) = U(0,b,x) */
        double Ua;
        double Uap1;
        int ap;
        double lm_for_local;
        stat_for = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_for_local); /* U(1,b,x) */
        Ua = r_Ua.val;
        Uam1 *= exp(-lm_for_local);
        lm_for.val = lm_for_local;
        lm_for.err = 0.0;

        for(ap=a0; ap<a1; ap++) {
          Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
          Uam1 = Ua;
          Ua   = Uap1;
          RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
        }
        Ua1_for_val  = Ua;
        Ua1_for_err  = fabs(Ua) * fabs(r_Ua.err/r_Ua.val);
        Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs(a1-a0)+1.0) * fabs(Ua1_for_val);
      }

      /* Now do the matching to produce the final result.
       */
      if(Ua1_bck_val == 0.0) {
        result->val = 0.0;
        result->err = 0.0;
        result->e10 = 0;
        GSL_ERROR ("error", GSL_EZERODIV);
      }
      else if(Ua1_for_val == 0.0) {
        /* Should never happen. */
        UNDERFLOW_ERROR_E10(result);
      }
      else {
        double lns = (scale_count_for - scale_count_bck)*log(scale_factor);
        double ln_for_val = log(fabs(Ua1_for_val));
        double ln_for_err = GSL_DBL_EPSILON + fabs(Ua1_for_err/Ua1_for_val);
        double ln_bck_val = log(fabs(Ua1_bck_val));
        double ln_bck_err = GSL_DBL_EPSILON + fabs(Ua1_bck_err/Ua1_bck_val);
        double lnr_val = lm_for.val + ln_for_val - ln_bck_val + lns;
        double lnr_err = lm_for.err + ln_for_err + ln_bck_err
          + 2.0 * GSL_DBL_EPSILON * (fabs(lm_for.val) + fabs(ln_for_val) + fabs(ln_bck_val) + fabs(lns));
        double sgn = GSL_SIGN(Ua1_for_val) * GSL_SIGN(Ua1_bck_val);
        int stat_e = gsl_sf_exp_err_e10_e(lnr_val, lnr_err, result);
        result->val *= sgn;
        return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
      }
    }
  }
}

/* Handle b >= 1 for generic a,b values.
 */
static
int
hyperg_U_bge1(const double a, const double b, const double x,
              gsl_sf_result_e10 * result)
{
  const double rinta = floor(a+0.5);
  const int a_neg_integer = (a < 0.0 && fabs(a - rinta) < INT_THRESHOLD);

  if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a_neg_integer && fabs(rinta) < INT_MAX) {
    /* U(-n,b,x) = (-1)^n n! Laguerre[n,b-1,x]
     */
    const int n = -(int)rinta;
    const double sgn = (GSL_IS_ODD(n) ? -1.0 : 1.0);
    gsl_sf_result lnfact;
    gsl_sf_result L;
    const int stat_L = gsl_sf_laguerre_n_e(n, b-1.0, x, &L);
    gsl_sf_lnfact_e(n, &lnfact);
    {
      const int stat_e = gsl_sf_exp_mult_err_e10_e(lnfact.val, lnfact.err,
                                                      sgn*L.val, L.err,
                                                      result);
      return GSL_ERROR_SELECT_2(stat_e, stat_L);
    }
  }
  else if(ASYMP_EVAL_OK(a,b,x)) {
    const double ln_pre_val = -a*log(x);
    const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
    gsl_sf_result asymp;
    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
    int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
                                              asymp.val, asymp.err,
                                              result);
    return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
  }
  else if(fabs(a) <= 1.0) {
    gsl_sf_result rU;
    double ln_multiplier;
    int stat_U = hyperg_U_small_a_bgt0(a, b, x, &rU, &ln_multiplier);
    int stat_e = gsl_sf_exp_mult_err_e10_e(ln_multiplier, 2.0*GSL_DBL_EPSILON*fabs(ln_multiplier), rU.val, rU.err, result);
    return GSL_ERROR_SELECT_2(stat_U, stat_e);
  }
  else if(SERIES_EVAL_OK(a,b,x)) {
    gsl_sf_result ser;
    const int stat_ser = hyperg_U_series(a, b, x, &ser);
    result->val = ser.val;
    result->err = ser.err;
    result->e10 = 0;
    return stat_ser;
  }
  else if(a < 0.0) {
    /* Recurse backward on a and then upward on b.
     */
    const double scale_factor = GSL_SQRT_DBL_MAX;
    const double a0 = a - floor(a) - 1.0;
    const double b0 = b - floor(b) + 1.0;
    int scale_count = 0;
    double lm_0, lm_1;
    double lm_max;
    gsl_sf_result r_Uap1;
    gsl_sf_result r_Ua;
    int stat_0 = hyperg_U_small_a_bgt0(a0+1.0, b0, x, &r_Uap1, &lm_0);
    int stat_1 = hyperg_U_small_a_bgt0(a0,     b0, x, &r_Ua,   &lm_1);
    int stat_e;
    double Uap1 = r_Uap1.val;
    double Ua   = r_Ua.val;
    double Uam1;
    double ap;
    lm_max = GSL_MAX(lm_0, lm_1);
    Uap1 *= exp(lm_0-lm_max);
    Ua   *= exp(lm_1-lm_max);

    /* Downward recursion on a.
     */
    for(ap=a0; ap>a+0.1; ap -= 1.0) {
      Uam1 = ap*(b0-ap-1.0)*Uap1 + (x+2.0*ap-b0)*Ua;
      Uap1 = Ua;
      Ua   = Uam1;
      RESCALE_2(Ua,Uap1,scale_factor,scale_count);
    }

    if(b < 2.0) {
      /* b == b0, so no recursion necessary
       */
      const double lnscale = log(scale_factor);
      gsl_sf_result lnm;
      gsl_sf_result y;
      lnm.val = lm_max + scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + scale_count * fabs(lnscale));
      y.val  = Ua;
      y.err  = fabs(r_Uap1.err/r_Uap1.val) * fabs(Ua);
      y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
      y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
      y.err *= fabs(lm_0-lm_max) + 1.0;
      y.err *= fabs(lm_1-lm_max) + 1.0;
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    }
    else {
      /* Upward recursion on b.
       */
      const double err_mult = fabs(b-b0) + fabs(a-a0) + 1.0;
      const double lnscale = log(scale_factor);
      gsl_sf_result lnm;
      gsl_sf_result y;

      double Ubm1 = Ua;                                 /* U(a,b0)   */
      double Ub   = (a*(b0-a-1.0)*Uap1 + (a+x)*Ua)/x;   /* U(a,b0+1) */
      double Ubp1;
      double bp;
      for(bp=b0+1.0; bp<b-0.1; bp += 1.0) {
        Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
        Ubm1 = Ub;
        Ub   = Ubp1;
        RESCALE_2(Ub,Ubm1,scale_factor,scale_count);
      }

      lnm.val = lm_max + scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
      y.val = Ub;
      y.err  = 2.0 * err_mult * fabs(r_Uap1.err/r_Uap1.val) * fabs(Ub);
      y.err += 2.0 * err_mult * fabs(r_Ua.err/r_Ua.val) * fabs(Ub);
      y.err += 2.0 * GSL_DBL_EPSILON * err_mult * fabs(Ub);
      y.err *= fabs(lm_0-lm_max) + 1.0;
      y.err *= fabs(lm_1-lm_max) + 1.0;
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    }
    return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  }
  else if(b >= 2*a + x) {
    /* Recurse forward from a near zero.
     * Note that we cannot cross the singularity at
     * the line b=a+1, because the only way we could
     * be in that little wedge is if a < 1. But we
     * have already dealt with the small a case.
     */
    int scale_count = 0;
    const double a0 = a - floor(a);
    const double scale_factor = GSL_SQRT_DBL_MAX;
    double lnscale;
    double lm_0, lm_1, lm_max;
    gsl_sf_result r_Uam1;
    gsl_sf_result r_Ua;
    int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
    int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
    int stat_e;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double Uam1 = r_Uam1.val;
    double Ua   = r_Ua.val;
    double Uap1;
    double ap;
    lm_max = GSL_MAX(lm_0, lm_1);
    Uam1 *= exp(lm_0-lm_max);
    Ua   *= exp(lm_1-lm_max);

    for(ap=a0; ap<a-0.1; ap += 1.0) {
      Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
      Uam1 = Ua;
      Ua   = Uap1;
      RESCALE_2(Ua,Uam1,scale_factor,scale_count);
    }

    lnscale = log(scale_factor);
    lnm.val = lm_max + scale_count * lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
    y.val  = Ua;
    y.err  = fabs(r_Uam1.err/r_Uam1.val) * fabs(Ua);
    y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
    y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
    y.err *= fabs(lm_0-lm_max) + 1.0;
    y.err *= fabs(lm_1-lm_max) + 1.0;
    stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  }
  else {
    if(b <= x) {
      /* Recurse backward to a near zero.
       */
      const double a0 = a - floor(a);
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count = 0;
      gsl_sf_result lnm;
      gsl_sf_result y;
      double lnscale;
      double lm_0;
      double Uap1;
      double Ua;
      double Uam1;
      gsl_sf_result U0;
      double ap;
      double ru;
      double r;
      int CF1_count;
      int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
      int stat_U0;
      int stat_e;
      r = ru/a;
      Ua   = GSL_SQRT_DBL_MIN;
      Uap1 = r * Ua;
      for(ap=a; ap>a0+0.1; ap -= 1.0) {
        Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
        Uap1 = Ua;
        Ua   = Uam1;
        RESCALE_2(Ua,Uap1,scale_factor,scale_count);
      }

      stat_U0 = hyperg_U_small_a_bgt0(a0, b, x, &U0, &lm_0);

      lnscale = log(scale_factor);
      lnm.val = lm_0 - scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_0) + fabs(scale_count * lnscale));
      y.val  = GSL_SQRT_DBL_MIN*(U0.val/Ua);
      y.err  = GSL_SQRT_DBL_MIN*(U0.err/fabs(Ua));
      y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a0-a) + CF1_count + 1.0) * fabs(y.val);
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
      return GSL_ERROR_SELECT_3(stat_e, stat_U0, stat_CF1);
    }
    else {
      /* Recurse backward to near the b=2a+x line, then
       * forward from a near zero to get the normalization.
       */
      int scale_count_for = 0;
      int scale_count_bck = 0;
      const double scale_factor = GSL_SQRT_DBL_MAX;
      const double eps = a - floor(a);
      const double a0 = ( eps == 0.0 ? 1.0 : eps );
      const double a1 = a0 + ceil(0.5*(b-x) - a0);
      gsl_sf_result lnm;
      gsl_sf_result y;
      double lm_for;
      double lnscale;
      double Ua1_bck;
      double Ua1_for;
      int stat_for;
      int stat_bck;
      int stat_e;
      int CF1_count;

      {
        /* Recurse back to determine U(a1,b), sans normalization.
         */
        double Uap1;
        double Ua;
        double Uam1;
        double ap;
        double ru;
        double r;
        int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
        r = ru/a;
        Ua   = GSL_SQRT_DBL_MIN;
        Uap1 = r * Ua;
        for(ap=a; ap>a1+0.1; ap -= 1.0) {
          Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
          Uap1 = Ua;
          Ua   = Uam1;
          RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
        }
        Ua1_bck = Ua;
        stat_bck = stat_CF1;
      }
      {
        /* Recurse forward to determine U(a1,b) with
         * absolute normalization.
         */
        gsl_sf_result r_Uam1;
        gsl_sf_result r_Ua;
        double lm_0, lm_1;
        int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
        int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
        double Uam1 = r_Uam1.val;
        double Ua   = r_Ua.val;
        double Uap1;
        double ap;

        lm_for = GSL_MAX(lm_0, lm_1);
        Uam1 *= exp(lm_0 - lm_for);
        Ua   *= exp(lm_1 - lm_for);

        for(ap=a0; ap<a1-0.1; ap += 1.0) {
          Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
          Uam1 = Ua;
          Ua   = Uap1;
          RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
        }
        Ua1_for = Ua;
        stat_for = GSL_ERROR_SELECT_2(stat_0, stat_1);
      }

      lnscale = log(scale_factor);
      lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs(scale_count_for - scale_count_bck)*fabs(lnscale));
      y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck;
      y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val);
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
      return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
    }
  }
}

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

int
gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x,
                             gsl_sf_result_e10 * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR_E10(result);
  }
  else {
    if(b >= 1) {
      return hyperg_U_int_bge1(a, b, x, result);
    }
    else {
      /* Use the reflection formula
       * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
       */
      gsl_sf_result_e10 U;
      double ln_x = log(x);
      int ap = 1 + a - b;
      int bp = 2 - b;
      int stat_e;
      int stat_U = hyperg_U_int_bge1(ap, bp, x, &U);
      double ln_pre_val = (1.0-b)*ln_x;
      double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs(b)+1.0) * fabs(ln_x);
      ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */
      stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
                                            U.val, U.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_U);
    }
  }
}

int
gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x,
                         gsl_sf_result_e10 * result)
{
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const int a_integer = ( fabs(a - rinta) < INT_THRESHOLD );
  const int b_integer = ( fabs(b - rintb) < INT_THRESHOLD );

  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR_E10(result);
  }
  else if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a_integer && b_integer) {
    return gsl_sf_hyperg_U_int_e10_e(rinta, rintb, x, result);
  }
  else {
    if(b >= 1.0) {
      /* Use b >= 1 function.
       */
      return hyperg_U_bge1(a, b, x, result);
    }
    else {
      /* Use the reflection formula
       * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
       */
      const double lnx = log(x);
      const double ln_pre_val = (1.0-b)*lnx;
      const double ln_pre_err = fabs(lnx) * 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(b));
      const double ap = 1.0 + a - b;
      const double bp = 2.0 - b;
      gsl_sf_result_e10 U;
      int stat_U = hyperg_U_bge1(ap, bp, x, &U);
      int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
                                            U.val, U.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_U);
    }
  }
}

int
gsl_sf_hyperg_U_int_e(const int a, const int b, const double x, gsl_sf_result * result)
{
  gsl_sf_result_e10 re;
  int stat_U = gsl_sf_hyperg_U_int_e10_e(a, b, x, &re);
  int stat_c = gsl_sf_result_smash_e(&re, result);
  return GSL_ERROR_SELECT_2(stat_c, stat_U);
}

int
gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result)
{
  gsl_sf_result_e10 re;
  int stat_U = gsl_sf_hyperg_U_e10_e(a, b, x, &re);
  int stat_c = gsl_sf_result_smash_e(&re, result);
  return GSL_ERROR_SELECT_2(stat_c, stat_U);
}

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

#include "eval.h"

double gsl_sf_hyperg_U_int(const int a, const int b, const double x)
{
  EVAL_RESULT(gsl_sf_hyperg_U_int_e(a, b, x, &result));
}

double gsl_sf_hyperg_U(const double a, const double b, const double x)
{
  EVAL_RESULT(gsl_sf_hyperg_U_e(a, b, x, &result));
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef MODULE_WASM_LINEAR_MEMORY_H
#define MODULE_WASM_LINEAR_MEMORY_H
#include <array>
#include <vector>
#include <algorithm>
#include "wasm_base.h"
#include "wasm_value.h"
#include <stdexcept>
#include <iostream>
#include <cstring>
#include <gsl/span>
#include "utilities/endianness.h"

namespace wasm {

struct WasmLinearMemory
{
	static constexpr const std::size_t page_size = 65536;
	using page_type = alignas(page_size) char[page_size];
private:
	using vector_type = wasm::SimpleVector<page_type>;
public:
	using wasm_external_kind_type = parse::Memory;

	using page_iterator = page_type*;
	using const_page_iterator = const page_type*;
	using page_pointer = page_type*;
	using const_page_pointer = const page_type*;
	using page_reference = page_type&;
	using const_page_reference = const page_type&;

	using page_span = gsl::span<const page_type>;
	using const_page_span = gsl::span<const page_type>;

	using value_type = char;
	using iterator = value_type*;
	using const_iterator = const value_type*;
	using pointer = value_type*;
	using const_pointer = const value_type*;
	using reference = value_type&;
	using const_reference = const value_type&;

	using raw_span = gsl::span<const value_type>;
	using const_raw_span = gsl::span<const value_type>;

	using size_type = vector_type::size_type;
	using difference_type = vector_type::difference_type;

	WasmLinearMemory(const parse::Memory& def, const parse::DataSegment& seg);

	friend const_page_span pages(const WasmLinearMemory& self)
	{ return const_page_span(self.memory_.data(), self.memory_.size()); }

	friend page_span pages(WasmLinearMemory& self)
	{ return page_span(self.memory_.data(), self.memory_.size()); }

	friend const_raw_span data(const WasmLinearMemory& self)
	{
		return raw_span(
			static_cast<const_pointer>(self.vec_.data()),
			page_size * self.vec_.size()
		);
	}

	friend raw_span data(WasmLinearMemory& self)
	{
		return raw_span(
			static_cast<pointer>(self.vec_.data()),
			page_size * self.vec_.size()
		);
	}

	friend wasm_sint32_t grow_memory(WasmLinearMemory& self, wasm_uint32_t delta)
	{
		wasm_uint32_t prev = static_cast<wasm_uint32_t>(pages(self).size());
		assert(prev == pages(self).size());
		if(delta > 0u)
		{
			std::size_t new_size = self.vec_.size() + delta;
			if(new_size > maximum_)
				return -1;
			try
			{
				self.resize(self.vec_.size() + delta);
			}
			catch(std::bad_alloc& e)
			{
				return -1;
			}
		}
		return reinterpret_cast<const wasm_sint32_t&>(prev);
	}

	friend bool matches(const WasmLinearMemory& self, const parse::Memory& tp)
	{
		return self.memory_.size() >= tp.size() and (
			self.maximum_ == tp.maximum.value_or(std::numeric_limits<std::size_t>::max())
		);
	}

private:
	void resize(std::size_t n)
	{
		assert(n <= maximum_);
		using std::swap;
		vector_type tmp(n);
		swap(memory_, tmp);
		std::memcpy(v.data(), tmp.data(), tmp.size() * sizeof(tmp.front()));
		std::memset(v.data() + tmp.size(), 0, v.size() - tmp.size());
	}
	vector_type memory_;
	const std::size_t maximum_ = std::numeric_limits<std::size_t>::max();
};

LanguageType index_type(const WasmLinearMemory& self)
{
	return
}

WasmLinearMemory::size_type current_memory(const WasmLinearMemory& self)
{ return pages(self).size(); }

WasmLinearMemory::size_type page_count(const WasmLinearMemory& self)
{ return current_memory(self); }

WasmLinearMemory::size_type size(const WasmLinearMemory& self)
{ return data(self).size(); }

WasmLinearMemory::WasmLinearMemory(const parse::Memory& def):
	memory_(def.initial),
	maximum_(def.maximum.value_or(std::numeric_limits<std::size_t>::max())
{
	auto bytes = data(*this);
	std::fill(bytes.begin(), bytes.end(), 0);
}

gsl::span<const char> compute_effective_address(
	const WasmLinearMemory& self,
	wasm_uint32_t base,
	wasm_uint32_t offset,
	std::size_t size
)
{
	auto mem = data(self);
	if(mem.size() <= base)
		throw std::out_of_range("Base address is too large while computing linear memory effective address.");
	std::size_t effective_address = base;
	effective_address += offset;
	if(mem.size() <= effective_address)
		throw std::out_of_range("Computed effective address is too large for linear memory.");
	if(std::size_t remaining = mem.size() - effective_address; remaining < sizeof(Type))
		throw std::out_of_range("Linear memory access partially acesses out-of-bounds memory.");
	return mem.subspan(effective_address, size);
}

gsl::span<char> compute_effective_address(
	WasmLinearMemory& self,
	wasm_uint32_t base,
	wasm_uint32_t offset,
	std::size_t size
)
{
	auto addr = compute_effective_address(std::as_const(self), base, offset, size);
	const char* data = addr.data();
	auto size = addr.data();
	return gsl::span<char>(const_cast<char*>(data), size);
}

template <class Type>
Type load_little_endian(const WasmLinearMemory& self, wasm_uint32_t base, wasm_uint32_t offset)
{
	static_assert(std::is_trivially_copyable_v<Type>);
	static_assert(std::is_arithmetic_v<Type>);
	Type result;
	auto addr = compute_effective_address(self, base, offset, sizeof(Type));
	assert(addr.size() == sizeof(result));
	std::memcpy(&result, addr.data(), sizeof(result));
	if(not system_is_little_endian())
	{
		static_assert(std::is_same_v<decltype(byte_swap(value)), decltype(value)>);
		result = byte_swap(result);
	}
	return result;
}

template <class Type, class PassedType>
void store_little_endian(const WasmLinearMemory& self, wasm_uint32_t base, wasm_uint32_t offset, PassedType value)
{
	static_assert(std::is_trivially_copyable_v<Type>);
	static_assert(std::is_arithmetic_v<Type>);
	static_assert(std::is_same_v<Type, PassedType>);
	auto pos = compute_effective_address(self, base, offset, sizeof(Type));
	assert(addr.size() == sizeof(value));
	if(not system_is_little_endian())
	{
		static_assert(std::is_same_v<decltype(byte_swap(value)), decltype(value)>);
		value = byte_swap(value);
	}
	std::memcpy(pos, &value, sizeof(value));
}

} /* namespace wasm */

#endif /* MODULE_WASM_LINEAR_MEMORY_H */

ext: h

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
// Copyright (C) 2021 Intel Corporation
// SPDX-License-Identifier: Apache-2.0

#ifndef _HEBench_ClearText_LogReg_H_7e5fa8c2415240ea93eff148ed73539b
#define _HEBench_ClearText_LogReg_H_7e5fa8c2415240ea93eff148ed73539b

#include <gsl/gsl>

#include "hebench/api_bridge/cpp/hebench.hpp"

#include "clear_benchmark.h"

template <class T>
class LogReg_Benchmark : public ClearTextBenchmark
{
private:
    HEBERROR_DECLARE_CLASS_NAME(LogReg_Benchmark)

public:
    LogReg_Benchmark(hebench::cpp::BaseEngine &engine,
                     const hebench::APIBridge::BenchmarkDescriptor &bench_desc,
                     const hebench::APIBridge::WorkloadParams &bench_params);
    ~LogReg_Benchmark() override;

    hebench::APIBridge::Handle encode(const hebench::APIBridge::DataPackCollection *p_parameters) override;
    void decode(hebench::APIBridge::Handle encoded_data, hebench::APIBridge::DataPackCollection *p_native) override;

    hebench::APIBridge::Handle load(const hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;
    void store(hebench::APIBridge::Handle remote_data,
               hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;

    hebench::APIBridge::Handle operate(hebench::APIBridge::Handle h_remote_packed,
                                       const hebench::APIBridge::ParameterIndexer *p_param_indexers) override;

protected:
    std::uint64_t m_vector_size;
    hebench::APIBridge::Workload m_workload;

private:
    static void logReg(gsl::span<T> &result,
                       const gsl::span<const T> &V0, const gsl::span<const T> &V1, const gsl::span<const T> &V2,
                       std::size_t element_count, hebench::APIBridge::Workload workload);
};

#include "inl/bench_logreg.inl"

#endif // defined _HEBench_ClearText_LogReg_H_7e5fa8c2415240ea93eff148ed73539b

max_stars_repo_name: hebench/backend-cpu-cleartext

max_issues_repo_name: hebench/backend-cpu-cleartext

max_forks_repo_name: hebench/backend-cpu-cleartext


text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <string.h>
#include <fftw3.h>

#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>

#include "../../../cosmolike_core/theory/basics.c"
#include "tinker_emulator.h"
#include "tinker_emulator.c"

int main(){
  //Test kernel
  double x0[7] = {0.0043*pow(0.7,2), 0.286*pow(0.7,2), -1, 1, 3.0, 70, 3.0};
  double x1[7] = {0.0043*pow(0.7,2), 0.286*pow(0.7,2), -1, 1, 3.0, 70, 5.0};
  printf("kernel : %e \n", kernel(x0, x1, emu_tinker_bias_log_lambda[0], tinkerEmuParam.tinker_bias_ncosmo));
  //Test bias emulator
  double emu_tinker_bias_param[6];
  predict_tinker_bias_parameters(0.0043*pow(0.7,2), 0.286*pow(0.7,2), -1, 1, 3.0, 70, 3.0, 0.3, emu_tinker_bias_param);
  double bias_ref[6] = {4.432419358797392, 0.32058653615384625, 1.3586846172228393, 1.347376928846154, 0.728081326923077, -1.0947315110189246};
  printf("bias result:\n");
  printf("item c ref\n");
  for(int i=0; i<tinkerEmuParam.tinker_bias_nparam_redshift; ++i) printf("%d, %e, %e\n", i,emu_tinker_bias_param[i], bias_ref[i]);
  printf("####################\n");

  //Test hmf emulator
  double emu_tinker_hmf_param[4];
  double hmf_ref[4] = { 2.5198116919192444, 0.8454997909106351, 0.5441108055593771, 1.3001900391246013};
  predict_tinker_hmf_parameters(0.0043*pow(0.7,2), 0.286*pow(0.7,2), -1, 1, 3.0, 70, 3.0, 0.3, emu_tinker_hmf_param);
  printf("hmf result:\n");
  printf("item c ref\n");
  for(int i=0; i<tinkerEmuParam.tinker_hmf_nparam_redshift; ++i) printf("%d, %e, %e\n", i,emu_tinker_hmf_param[i], hmf_ref[i]);
  printf("####################\n");

}

max_stars_repo_name: tmcclintock/TinkerEmulator

max_issues_repo_name: tmcclintock/TinkerEmulator

max_forks_repo_name: tmcclintock/TinkerEmulator


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
/*
 *
 *
 *       Filename:  total_current.c
 *
 *    Description:  Total current ( AC + DC )
 *
 *        Version:  1.0
 *        Created:  26/06/2014 14:59:12
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include "funcs.h"
#include <gsl/gsl_vector.h>

/*
 *      FUNCTION
 *         Name:  tot_current
 *  Description:
 *
 */
double tot_current ( const gsl_vector* rho, void* params )
{
	struct f_params* pars = (struct f_params*) params ;
	double Omega = pars->Omega ;
	double o_1   = pars->omega_1 ;
	double D     = sqrt(o_1*o_1 - Omega*Omega) ;
	double curr = -((VECTOR(rho, 2)*D + VECTOR(rho, 3)*Omega)/o_1) ;
	return (curr) ;
}		/* -----  end of function tot_current  ----- */


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void generate_MCMC_proposition(MCMC_info *MCMC, int flag_chain_init) {
    // Set the parameters for first trial parameters for the proposition
    switch(flag_chain_init) {
        case GBP_TRUE:
            memcpy(MCMC->P_new, MCMC->P_init, (size_t)MCMC->n_P * sizeof(double));
            memcpy(MCMC->P_last, MCMC->P_init, (size_t)MCMC->n_P * sizeof(double));
            memcpy(MCMC->P_chain, MCMC->P_init, (size_t)MCMC->n_P * sizeof(double));
            MCMC->flag_init_chain = GBP_FALSE;
            MCMC->n_success       = 0;
            MCMC->n_fail          = 0;
            MCMC->n_propositions  = 0;
            break;
        case GBP_FALSE:
            generate_MCMC_parameters(MCMC);
            break;
    }

    // Keep generating new parameter sets until the mapping function is satisfied
    while(MCMC->map_P_to_M(MCMC->P_new, MCMC, MCMC->M_new)) {
        MCMC->n_map_calls++;
        MCMC->first_map_call = GBP_FALSE;
        generate_MCMC_parameters(MCMC);
    }
    MCMC->first_map_call = GBP_FALSE;
    MCMC->n_map_calls++;

    // Produce likelihood for this proposition
    MCMC->compute_MCMC_ln_likelihood(
        MCMC, MCMC->M_new, MCMC->P_new, MCMC->ln_likelihood_DS, MCMC->n_DoF_DS, &(MCMC->ln_likelihood_new), &(MCMC->n_DoF));
    if(!SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_PARALLEL)) {
        SID_Bcast(MCMC->ln_likelihood_DS, MCMC->n_DS, SID_DOUBLE, SID_MASTER_RANK, MCMC->comm);
        SID_Bcast(MCMC->n_DoF_DS, MCMC->n_DS, SID_INT, SID_MASTER_RANK, MCMC->comm);
        SID_Bcast(&(MCMC->ln_likelihood_new), 1, SID_DOUBLE, SID_MASTER_RANK, MCMC->comm);
        SID_Bcast(&(MCMC->n_DoF), 1, SID_INT, SID_MASTER_RANK, MCMC->comm);
    }

    MCMC->first_likelihood_call = GBP_FALSE;
    if(flag_chain_init) {
        MCMC->ln_likelihood_chain = MCMC->ln_likelihood_new;
        MCMC->ln_likelihood_best  = MCMC->ln_likelihood_new;
        memcpy(MCMC->P_best, MCMC->P_new, (size_t)MCMC->n_P * sizeof(double));
        MCMC->ln_Pr_chain = 0.;
        MCMC->ln_Pr_new   = 0.;
    }
    // Keep track of the best proposition
    else if(MCMC->ln_likelihood_new > MCMC->ln_likelihood_best) {
        MCMC->ln_likelihood_best = MCMC->ln_likelihood_new;
        memcpy(MCMC->P_best, MCMC->P_new, (size_t)MCMC->n_P * sizeof(double));
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
 least_square_solver.c
 written by Eunseok Lee

 function: obtain the least-square solution of a linear equation system

 v1: Feb 2, 2018
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <cem.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

void least_square_solver(double *matA, int n_row, int n_col, double *y, double *sol_x) {

    int i, j, k, tmp;
    double sum;
    gsl_matrix * A = gsl_matrix_alloc (n_row,n_col);
    gsl_matrix * V = gsl_matrix_alloc (n_col,n_col);
    gsl_vector * S = gsl_vector_alloc (n_col);
    gsl_vector * work = gsl_vector_alloc (n_col);
//    gsl_matrix * U = gsl_matrix_alloc (n_row-1,n_col_ids);
    gsl_vector * b = gsl_vector_alloc (n_row);
    gsl_vector * x = gsl_vector_alloc (n_col);

    // need to copy matA to A, y to b!
    for (i=0;i<n_row;i++) {
        for (j=0; j<n_col; j++) {
            gsl_matrix_set (A, i, j, *(matA+i*n_col+j));
        }
        gsl_vector_set (b, i, y[i]);
    }
    gsl_linalg_SV_decomp(A, V, S, work);    // on output, A is replaced by U!
    gsl_linalg_SV_solve (A, V, S, b, x);    // So, A should be used instead of U, here.

    for (i=0;i<n_col;i++) {
        sol_x[i] = gsl_vector_get(x,i);
    }

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(x);

}


text:
#ifndef CONVOLUTION_H
#define CONVOLUTION_H

#include <cstddef>

#include <cstring>

#include <iostream>

extern "C" {
  #include <cblas.h>
}

using namespace std;

template<typename T>
constexpr T relu(T x) {
  return (x > 0) ? x : 0;
}

template<typename T>
void mul_mat(T *A, double *B, double *C, size_t m, size_t k, size_t n) {
  for(size_t r = 0; r < m; r++) {
    for(size_t c = 0; c < n; c++) {
      double x = 0.0;
      for(size_t z = 0; z < k; z++) {
        x += A[r * k + z] * B[z * n + c];
      }
      C[r * n + c] = x;
    }
  }
}

template<size_t in_maps, size_t in_height, size_t in_width,
         size_t conv_height, size_t conv_width,
         size_t vert_stride, size_t horiz_stride,
         size_t out_maps>
struct ConvolutionLayer {
  static constexpr size_t conv_size = conv_height * conv_width;
  static constexpr size_t in_map_size = in_height * in_width;
  static constexpr size_t inputs_no = in_maps * in_height * in_width;
  static constexpr size_t parameters_no = out_maps * (conv_width * conv_height * in_maps + 1);
  static constexpr size_t bias_idx = conv_width * conv_height * in_maps;
  static constexpr size_t out_height = (in_height - conv_height) / vert_stride +1;
  static constexpr size_t out_width = (in_width - conv_width) / horiz_stride + 1;
  static constexpr size_t out_map_size = out_height * out_width;
  static constexpr size_t internal_no = out_maps * out_height * out_width;
  static constexpr size_t outputs_no = out_maps * out_height * out_width;

  static constexpr size_t p_idx(size_t out_map, size_t in_map, size_t conv_r, size_t conv_c) {
    return out_map * (in_maps * conv_height * conv_width + 1) + in_map * (conv_height * conv_width)
      + conv_r * conv_width + conv_c;
  }

  static void forward(double* inputs, double * parameters, double* outputs) {
    for (size_t j = 0; j < out_maps; j++) {
      double *pj = &(parameters[j*(conv_height * conv_width * in_maps + 1)]);

      for (size_t out_r = 0; out_r < out_height; out_r++) {
        for (size_t out_c = 0; out_c < out_width; out_c++) {
          const size_t out_idx = j * (out_map_size) + out_r * out_width + out_c;
          double a = 0;
          for (size_t i = 0; i < in_maps; i++) {
            const size_t i_idx = i * conv_size;
            const size_t i2_idx = i * in_map_size;
            for (size_t conv_r = 0; conv_r < conv_height; conv_r++) {
              for (size_t conv_c = 0; conv_c < conv_width; conv_c++) {
                size_t p2_idx =  i_idx + conv_r * conv_width + conv_c;
                size_t in_idx =
                  i2_idx + (out_r * vert_stride + conv_r) * in_width + out_c * horiz_stride + conv_c;
                a += pj[p2_idx] * inputs[in_idx];
              }
            }
          }
          a += pj[bias_idx];
          outputs[out_idx] = relu(a);
        }
      }
    }
  }
};

template<size_t in_maps, size_t in_height, size_t in_width,
         size_t conv_height, size_t conv_width,
         size_t vert_stride, size_t horiz_stride,
         size_t out_maps>
struct ConvolutionLayer2 {
  static constexpr size_t in_map_size = in_height * in_width;
  static constexpr size_t inputs_no = in_maps * in_map_size;
  static constexpr size_t out_height = (in_height - conv_height) / vert_stride +1;
  static constexpr size_t out_width = (in_width - conv_width) / horiz_stride + 1;
  static constexpr size_t out_map_size = out_height * out_width;
  static constexpr size_t internal_no = out_maps * out_map_size;
  static constexpr size_t outputs_no = out_maps * out_map_size;
  static constexpr size_t conv_size = conv_height * conv_width;
  static constexpr size_t parameters_no = out_maps * (in_maps * conv_size + 1);

  static constexpr size_t p_idx(size_t out_map, size_t in_map, size_t conv_r, size_t conv_c) {
    return out_map * (in_maps * conv_height * conv_width + 1) + in_map * (conv_height * conv_width)
      + conv_r * conv_width + conv_c;
  }

  static void forward(double* inputs, double *parameters, double* outputs) {
    double* _in = (double*) malloc(sizeof(double) * in_maps * conv_size * out_map_size);

    for (size_t out_r = 0; out_r < out_height; out_r++) {
      for (size_t out_c = 0; out_c < out_width; out_c++) {
        for (size_t i = 0; i < in_maps; i++) {
          for (size_t conv_r = 0; conv_r < conv_height; conv_r++) {
            for (size_t conv_c = 0; conv_c < conv_width; conv_c++) {
              size_t in_idx =
                i * in_map_size +
                (out_r * vert_stride + conv_r) * in_width +
                out_c * horiz_stride + conv_c;
              size_t row_idx = i * conv_size + conv_r * conv_width + conv_c;
              size_t _in_idx = row_idx * out_map_size + out_r * out_width + out_c;
              _in[_in_idx] = inputs[in_idx];
            }
          }
        }
      }
    }

    for (size_t j = 0; j < out_maps;j++) {
      double bias = parameters[(j+1) * (in_maps * conv_size + 1) - 1];
      for (size_t i = 0; i < out_map_size; i++) {
        outputs[j * out_map_size + i] = bias;
      }
    }

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                out_maps, out_map_size, in_maps * conv_size,
                1.0, parameters, in_maps * conv_size + 1,
                _in, out_map_size,
                1.0, outputs, out_map_size);

    for(size_t j = 0; j < out_maps * out_map_size; j++) {
      outputs[j] = relu(outputs[j]);
    }
  }
};

#endif

max_stars_repo_name: tudor-berariu/cc-hacks

max_issues_repo_name: tudor-berariu/cc-hacks

max_forks_repo_name: tudor-berariu/cc-hacks


text:
#ifndef DERIVATIVES_H_5CHQ89V7
#define DERIVATIVES_H_5CHQ89V7

#include <gsl/gsl>
#include <tuple>
#include <type_traits>

namespace sens_loc::math {

/// Calculate the first derivate with the central differential quotient.
/// \tparam Real precision of the calculation
/// \param y__1 \f$y_{i-1}\f$
/// \param y_1 \f$y_{i+1}\f$
/// \param dx \f$2. * dx\f$
/// \returns first derivative at this point of order \f$\mathcal{O}(dx^2)\f$
template <typename Real>
inline Real first_derivative_central(Real y__1, Real y_1, Real dx) noexcept {
    static_assert(std::is_floating_point_v<Real>);

    Expects(dx > Real(0.));
    // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
    return (y_1 - y__1) / (Real(2.) * dx);
}

/// Calculate the second derivate with the central differential quotient.
/// \tparam Real precision of the calculation
/// \param y__1 \f$y_{i-1}\f$
/// \param y_0 \f$y_{i}\f$
/// \param y_1 \f$y_{i+1}\f$
/// \param dx \f$dx\f$
/// \returns second derivative at this point of order \f$\mathcal{O}(dx^2)\f$
template <typename Real>
inline Real
second_derivative_central(Real y__1, Real y_0, Real y_1, Real dx) noexcept {
    static_assert(std::is_floating_point_v<Real>);

    Expects(dx > Real(0.));
    // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
    return (y_1 + y__1 - Real(2.) * y_0) / (dx * dx);
}

/// Calculate the derivatives for a surface patch.
///
/// Index convention:
/// \f$d\_\_1 == d_{-1}\f$
/// \f$d\_\_0 == d_{0}\f$
/// \f$d\_1   == d_{1}\f$
///
/// Angle convention:
/// \f$\varphi\f$ -> u direction
/// \f$\theta\f$ -> v direction
///
/// \tparam Real precision of the calculation
/// \param d__1__1,d__1__0,d__1_1 neighbours "above" central pixel
/// \param d__0__1,d__0__0,d__0_1 same row as the central pixel
/// \param d_1__1,d_1__0,d_1_1 row "after" the central pixel
/// \param d_phi angle between rays in x direction \f$(u - 1, u + 1)\f$
/// \param d_theta angle between rays in y direction \f$(v - 1, v + 1)\f$
/// \param d_phi_theta angle between rays in diagonal direction
/// \f$(u - 1, v + 1)\f$
/// \returns partial derivatives \f$(f_u, f_v, f_uu, f_vv, f_uv)\f$
/// \pre the depth values shuold be positive, as they encode depth values
/// \pre \p / d_phi, \p d_theta, \p d_phi_theta are all positive angles
// clang-format off
template <typename Real = float>
inline std::tuple<Real, Real, Real, Real, Real>
derivatives(Real d__1__1, Real d__1__0, Real d__1_1,
            Real d__0__1, Real d__0__0, Real d__0_1,
            Real d_1__1, Real d_1__0, Real d_1_1,
            Real d_phi, Real d_theta, Real d_phi_theta) noexcept {
    static_assert(std::is_floating_point_v<Real>);

    Expects(d_phi > 0.);
    Expects(d_theta > 0.);
    Expects(d_phi_theta > 0.);

    (void)d__1__1;
    (void)d__1_1;
    (void)d_1__1;
    (void)d_1_1;

    // clang-format on
    const Real f_u = math::first_derivative_central(d__0__1, d__0_1, d_phi);
    const Real f_v = math::first_derivative_central(d__1__0, d_1__0, d_theta);
    const Real f_uu =
        math::second_derivative_central(d__0__1, d__0__0, d__0_1, d_phi);
    const Real f_vv =
        math::second_derivative_central(d__1__0, d__0__0, d_1__0, d_theta);
    const Real f_uv =
        math::second_derivative_central(d__1__0, d__0__0, d_1__0, d_phi_theta);

    return std::make_tuple(f_u, f_v, f_uu, f_vv, f_uv);
};
// clang-format on
}  // namespace sens_loc::math

#endif /* end of include guard: DERIVATIVES_H_5CHQ89V7 */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#include <slepcmfn.h>
#include "cs.h"
#include <gsl/gsl_rng.h>
//for use in programs like im_clam
typedef struct clam_lik_params{
	int snpNumber, maxSampleSize,nstates,n1,n2,nParams;
	gsl_vector *sampleSizeVector,*rates, *stateVec, *resVec ;
	gsl_vector *paramVector,*paramCIUpper,*paramCILower, *mlParams;
	//For straight MLE the paramVector takes the form [N2,NA,m1,m2,t]
	gsl_matrix *log_afs, *transMat,*expAFS, *expAFS2, *invMatGSL, *invMatGSLA, *obsData;
	afsStateSpace *stateSpace ,*reducedStateSpace;
	gsl_rng *rng;
	struct site *data;
	double *cs_work,*new,*topA,*resultTmp,**invMat,*b,*res, *st, *top;
	int *dim1, *dim2, *moveA, *move, *dim1A, *dim2A; //these pointers ease the construction of sparse matrices
	int nnz, nnzA;
	double maxLik;
	int *map, *reverseMap; //these provide mapping of between state spaces of two phases
	int optimized;
	int rank;
	int Na,fEvals;
	double *expoArray;
	double meanTreeLength;
	cs_di *spMat, *eye, *eyeAnc, *triplet, *tmpMat, *tmpMat2;

	Vec  x_0, x, u,bP, u2, bP2, v, y, v_seq, *x2, ancStateVec, ancResVec;
	Vec xInv,bInv;
	Mat            ident,C, C2, C_transpose, D, D_copy, F, denseIdent, denseMat1, denseMat2, denseMat3;
	VecScatter	ctx;
	MFN                mfn;
	PetscInt        iStart,iEnd;
      /* linear solver context */
	PC             pc;           /* preconditioner context */

}clam_lik_params;

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
#ifndef ALLVARS_H
#define ALLVARS_H

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "core_simulation.h"

#define ABORT(sigterm)                                                  \
do {                                                                \
  printf("Error in file: %s\tfunc: %s\tline: %i\n", __FILE__, __FUNCTION__, __LINE__); \
  myexit(sigterm);                                                \
} while(0)

#define  STEPS 10         // Number of integration intervals between two snapshots
#define  MAXGALFAC 1
#define  ALLOCPARAMETER 10.0
#define  MAX_NODE_NAME_LEN 50
#define  ABSOLUTEMAXSNAPS 1000

#define  GRAVITY     6.672e-8
#define  SOLAR_MASS  1.989e33
#define  SOLAR_LUM   3.826e33
#define  RAD_CONST   7.565e-15
#define  AVOGADRO    6.0222e23
#define  BOLTZMANN   1.3806e-16
#define  GAS_CONST   8.31425e7
#define  C           2.9979e10
#define  PLANCK      6.6262e-27
#define  CM_PER_MPC  3.085678e24
#define  PROTONMASS  1.6726e-24
#define  HUBBLE      3.2407789e-18   /* in h/sec */

#define  SEC_PER_MEGAYEAR   3.155e13
#define  SEC_PER_YEAR       3.155e7

// This structure contains the properties that are output
struct GALAXY_OUTPUT
{
  int   SnapNum;
  int   Type;

  long long   GalaxyIndex;
  long long   CentralGalaxyIndex;
  int   SAGEHaloIndex;
  int   SAGETreeIndex;
  long long   SimulationHaloIndex;

  int   mergeType;  //0=none; 1=minor merger; 2=major merger; 3=disk instability; 4=disrupt to ICS
  int   mergeIntoID;
  int   mergeIntoSnapNum;
  float dT;

  // (sub)halo properties
  float Pos[3];
  float Vel[3];
  float Spin[3];
  int   Len;
  float Mvir;
  float CentralMvir;
  float Rvir;
  float Vvir;
  float Vmax;
  float VelDisp;

  // baryonic reservoirs
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float BlackHoleMass;
  float ICS;

  // metals
  float MetalsColdGas;
  float MetalsStellarMass;
  float MetalsBulgeMass;
  float MetalsHotGas;
  float MetalsEjectedMass;
  float MetalsICS;

  // to calculate magnitudes
  float SfrDisk;
  float SfrBulge;
  float SfrDiskZ;
  float SfrBulgeZ;

  // misc
  float DiskScaleRadius;
  float Cooling;
  float Heating;
  float QuasarModeBHaccretionMass;
  float TimeOfLastMajorMerger;
  float TimeOfLastMinorMerger;
  float OutflowRate;

  // infall properties
  float infallMvir;
  float infallVvir;
  float infallVmax;
};

// This structure contains the properties used within the code
struct GALAXY
{
  int   SnapNum;
  int   Type;

  int   GalaxyNr;
  int   CentralGal;
  int   HaloNr;
  long long MostBoundID;

  int   mergeType;  //0=none; 1=minor merger; 2=major merger; 3=disk instability; 4=disrupt to ICS
  int   mergeIntoID;
  int   mergeIntoSnapNum;
  float dT;

  // (sub)halo properties
  float Pos[3];
  float Vel[3];
  int   Len;
  float Mvir;
  float deltaMvir;
  float CentralMvir;
  float Rvir;
  float Vvir;
  float Vmax;

  // baryonic reservoirs
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float BlackHoleMass;
  float ICS;

  // metals
  float MetalsColdGas;
  float MetalsStellarMass;
  float MetalsBulgeMass;
  float MetalsHotGas;
  float MetalsEjectedMass;
  float MetalsICS;

  // to calculate magnitudes
  float SfrDisk[STEPS];
  float SfrBulge[STEPS];
  float SfrDiskColdGas[STEPS];
  float SfrDiskColdGasMetals[STEPS];
  float SfrBulgeColdGas[STEPS];
  float SfrBulgeColdGasMetals[STEPS];

  // misc
  float DiskScaleRadius;
  float MergTime;
  double Cooling;
  double Heating;
  float r_heat;
  float QuasarModeBHaccretionMass;
  float TimeOfLastMajorMerger;
  float TimeOfLastMinorMerger;
  float OutflowRate;
  float TotalSatelliteBaryons;

  // infall properties
  float infallMvir;
  float infallVvir;
  float infallVmax;
}
*Gal, *HaloGal;

// auxiliary halo data
struct halo_aux_data
{
  int DoneFlag;
  int HaloFlag;
  int NGalaxies;
  int FirstGalaxy;
}
*HaloAux;

extern int    FirstFile;    // first and last file for processing
extern int    LastFile;

extern int    Ntrees;      // number of trees in current file
extern int    NumGals;     // Total number of galaxies stored for current tree
extern int    MaxGals;     // Maximum number of galaxies allowed for current tree
extern int    FoF_MaxGals;

extern int    GalaxyCounter;     // unique galaxy ID for main progenitor line in tree

extern int    LastSnapShotNr;

extern char   OutputDir[512];
extern char   FileNameGalaxies[512];
extern char   TreeName[512];
extern char   SimulationDir[512];
extern char   FileWithSnapList[512];

extern int    TotHalos;
extern int    TotGalaxies[ABSOLUTEMAXSNAPS];
extern int    *TreeNgals[ABSOLUTEMAXSNAPS];

extern int    *FirstHaloInSnap;

extern int    *TreeNHalos;
extern int    *TreeFirstHalo;

#ifdef MPI
extern int ThisTask, NTask, nodeNameLen;
extern char *ThisNode;
#endif

extern double Omega;
extern double OmegaLambda;
extern double PartMass;
extern double Hubble_h;
extern double EnergySNcode, EnergySN;
extern double EtaSNcode, EtaSN;

// recipe flags
extern int    ReionizationOn;
extern int    SupernovaRecipeOn;
extern int    DiskInstabilityOn;
extern int    AGNrecipeOn;
extern int    SFprescription;

// recipe parameters
extern double RecycleFraction;
extern double Yield;
extern double FracZleaveDisk;
extern double ReIncorporationFactor;
extern double ThreshMajorMerger;
extern double BaryonFrac;
extern double SfrEfficiency;
extern double FeedbackReheatingEpsilon;
extern double FeedbackEjectionEfficiency;
extern double RadioModeEfficiency;
extern double QuasarModeEfficiency;
extern double BlackHoleGrowthRate;
extern double Reionization_z0;
extern double Reionization_zr;
extern double ThresholdSatDisruption;

extern double UnitLength_in_cm,
  UnitTime_in_s,
  UnitVelocity_in_cm_per_s,
  UnitMass_in_g,
  RhoCrit,
  UnitPressure_in_cgs,
  UnitDensity_in_cgs,
  UnitCoolingRate_in_cgs,
  UnitEnergy_in_cgs,
  UnitTime_in_Megayears,
  G,
  Hubble,
  a0, ar;

extern int    ListOutputSnaps[ABSOLUTEMAXSNAPS];

extern double ZZ[ABSOLUTEMAXSNAPS];
extern double AA[ABSOLUTEMAXSNAPS];
extern double *Age;

extern int    MAXSNAPS;
extern int    NOUT;
extern int    Snaplistlen;

extern gsl_rng *random_generator;

extern int TreeID;
extern int FileNum;

#endif  // #ifndef ALLVARS_H

max_stars_repo_name: jacobseiler/SAGE

max_issues_repo_name: jacobseiler/SAGE

max_forks_repo_name: jacobseiler/SAGE


text:
#include <petsc.h>
#include <petscts2.h>

#if PETSC_VERSION_LE(3,3,0)
#define TSRegister(s,f) TSRegister(s,0,0,f)
#endif

typedef struct {
  PetscReal Omega; /* frequency */
  PetscReal Xi;    /* damping   */
} UserParams;

#undef  __FUNCT__
#define __FUNCT__ "Residual1"
PetscErrorCode Residual1(TS ts,PetscReal t,Vec X,Vec A,Vec R,void *ctx)
{
  UserParams *user = (UserParams *)ctx;
  PetscReal Omega = user->Omega;
  PetscScalar *x,*a,*r;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(A,&a);CHKERRQ(ierr);
  ierr = VecGetArray(R,&r);CHKERRQ(ierr);

  r[0] = a[0] + (Omega*Omega)*x[0];

  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(A,&a);CHKERRQ(ierr);
  ierr = VecRestoreArray(R,&r);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(R);CHKERRQ(ierr);
  ierr = VecAssemblyEnd  (R);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "Tangent1"
PetscErrorCode Tangent1(TS ts,PetscReal t,Vec X,Vec A,PetscReal shiftA,Mat J,Mat P,void *ctx)
{
  UserParams *user = (UserParams *)ctx;
  PetscReal Omega = user->Omega;
  PetscReal T = 0;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  T = shiftA + (Omega*Omega);

  ierr = MatSetValue(P,0,0,T,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd  (P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (J != P) {
    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}
#if PETSC_VERSION_LT(3,5,0)
PetscErrorCode Tangent1_Legacy(TS ts,PetscReal t,Vec U,Vec V,PetscReal shift,Mat *J,Mat *P,MatStructure *m,void *ctx)
{*m = SAME_NONZERO_PATTERN; return Tangent1(ts,t,U,V,shift,*J,*P,ctx);}
#define Tangent1 Tangent1_Legacy
#endif

#undef  __FUNCT__
#define __FUNCT__ "Residual2"
PetscErrorCode Residual2(TS ts,PetscReal t,Vec X,Vec V,Vec A,Vec R,void *ctx)
{
  UserParams *user = (UserParams *)ctx;
  PetscReal Omega = user->Omega, Xi = user->Xi;
  PetscScalar *x,*v,*a,*r;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(V,&v);CHKERRQ(ierr);
  ierr = VecGetArray(A,&a);CHKERRQ(ierr);
  ierr = VecGetArray(R,&r);CHKERRQ(ierr);

  r[0] = a[0] + (2*Xi*Omega)*v[0] + (Omega*Omega)*x[0];

  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(V,&v);CHKERRQ(ierr);
  ierr = VecRestoreArray(A,&a);CHKERRQ(ierr);
  ierr = VecRestoreArray(R,&r);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(R);CHKERRQ(ierr);
  ierr = VecAssemblyEnd  (R);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "Tangent2"
PetscErrorCode Tangent2(TS ts,PetscReal t,Vec X,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P,void *ctx)
{
  UserParams *user = (UserParams *)ctx;
  PetscReal Omega = user->Omega, Xi = user->Xi;
  PetscReal T = 0;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  T = shiftA + shiftV * (2*Xi*Omega) + (Omega*Omega);

  ierr = MatSetValue(P,0,0,T,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd  (P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (J != P) {
    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "Monitor"
PetscErrorCode Monitor(TS ts,PetscInt i,PetscReal t,Vec U,void *ctx)
{
  const char *filename = (const char *)ctx;
  static FILE *fp = 0;
  Vec X,V;
  const PetscScalar *x,*v;
  TSConvergedReason reason;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  if (i<0) PetscFunctionReturn(0); /*XXX temporary petsc-dev fix */

  if (!fp) {ierr = PetscFOpen(PETSC_COMM_SELF,filename,"w",&fp);CHKERRQ(ierr);}
  ierr = TSGetSolution2(ts,&X,&V);CHKERRQ(ierr);
  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecGetArrayRead(V,&v);CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%g %g %g\n",(double)t,(double)x[0],(double)v[0]);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(V,&v);CHKERRQ(ierr);
  ierr = TSGetConvergedReason(ts,&reason); CHKERRQ(ierr);
  if (reason) {ierr = PetscFClose(PETSC_COMM_SELF,fp);CHKERRQ(ierr); fp=0;}

  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
PetscErrorCode TSCreate_Alpha2(TS);
EXTERN_C_END

#undef  __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char *argv[]) {

  TS             ts;
  Vec            R;
  Mat            J;
  Vec            X,V;
  PetscScalar    *x,*v;
  UserParams     user;
  PetscBool      out;
  char           output[PETSC_MAX_PATH_LEN] = {0};
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
  ierr = TSRegister(TSALPHA2,TSCreate_Alpha2);CHKERRQ(ierr);

  user.Omega = 1.0;
  user.Xi    = 0.0;
  ierr = PetscOptionsBegin(PETSC_COMM_SELF,"","Oscillator Options","TS");CHKERRQ(ierr);
  ierr = PetscOptionsReal("-frequency","Frequency",__FILE__,user.Omega,&user.Omega,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-damping",  "Damping",  __FILE__,user.Xi,   &user.Xi,   NULL);CHKERRQ(ierr);
  ierr = PetscOptionsString("-output","Output",__FILE__,output,output,sizeof(output),&out);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  if (out && !output[0]) {ierr = PetscStrcpy(output,"Oscillator.out");CHKERRQ(ierr);}

  ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr);
  ierr = TSSetType(ts,TSALPHA2);CHKERRQ(ierr);
  ierr = TSSetDuration(ts,PETSC_MAX_INT,2*M_PI * 5);CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);

  ierr = VecCreateSeq(PETSC_COMM_SELF,1,&R);CHKERRQ(ierr);
  ierr = VecSetUp(R);CHKERRQ(ierr);
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,1,1,NULL,&J);CHKERRQ(ierr);
  ierr = MatSetUp(J);CHKERRQ(ierr);
  if (user.Xi <= 0.0) {
    ierr = TSSetIFunction(ts,R,Residual1,&user);CHKERRQ(ierr);
    ierr = TSSetIJacobian(ts,J,J,Tangent1,&user);CHKERRQ(ierr);
  } else {
    ierr = TSSetIFunction2(ts,R,Residual2,&user);CHKERRQ(ierr);
    ierr = TSSetIJacobian2(ts,J,J,Tangent2,&user);CHKERRQ(ierr);
  }
  ierr = VecDestroy(&R);CHKERRQ(ierr);
  ierr = MatDestroy(&J);CHKERRQ(ierr);

  if (output[0]) {
    ierr = TSMonitorSet(ts,Monitor,output,NULL);CHKERRQ(ierr);
  }

  ierr = VecCreateSeq(PETSC_COMM_SELF,1,&X);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,1,&V);CHKERRQ(ierr);
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  ierr = VecGetArray(V,&v);CHKERRQ(ierr);
  x[0] = 1.0;
  v[0] = 0.0;
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(V,&v);CHKERRQ(ierr);

  ierr = TSSetSolution2(ts,X,V);CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
  ierr = TSSolve2(ts,X,V);CHKERRQ(ierr);

  ierr = VecDestroy(&X);CHKERRQ(ierr);
  ierr = VecDestroy(&V);CHKERRQ(ierr);
  ierr = TSDestroy(&ts);CHKERRQ(ierr);
  ierr = PetscFinalize();CHKERRQ(ierr);

  return 0;
}

max_stars_repo_name: otherlab/petiga

max_issues_repo_name: otherlab/petiga

max_forks_repo_name: otherlab/petiga


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef treeop_00fc8cde_c6b6_464d_b51d_71eaad31b9ef_h
#define treeop_00fc8cde_c6b6_464d_b51d_71eaad31b9ef_h

#include <gslib/tree.h>

__gslib_begin__

/*
 * The original tree was created from top to bottom, to implement the new operation it might be needed
 * to create the tree from bottom to top.
 * These classes were designed to achieve that.
 * Caution:
 * You'd better NOT include me in headers.
 * You'd better use def_tr$ macros in a function scope.
 * with argumented def_tr$ versions, IE. like def_tr$_v1, def_tr$_v2, and so on...
 * the tree should NOT be a copy wrapper one.
 */

template<class _tree>
class treeop_helper
{
public:
    typedef _tree tree;
    typedef typename _tree::value value;
    typedef typename _tree::wrapper wrapper;
    typedef typename _tree::alloc alloc;
    typedef typename _tree::const_iterator const_iterator;
    typedef typename _tree::iterator iterator;
    typedef typename wrapper::children children;
};

template<class _tree>
class treeop_set
{
public:
    typedef treeop_helper<_tree> helper;
    typedef typename helper::tree tree;
    typedef typename helper::value value;
    typedef typename helper::wrapper wrapper;
    typedef typename helper::alloc alloc;
    typedef typename helper::const_iterator const_iterator;
    typedef typename helper::iterator iterator;
    typedef typename helper::children children;
    typedef children opset;

public:
    opset           _opset;

public:
    treeop_set() {}
    treeop_set(wrapper* w) { _opset.init(w); }
    treeop_set(treeop_set& rhs) { _opset.swap(rhs.dump()); }
    bool empty() const { return _opset.empty(); }
    bool unary() const { return _opset.size() == 1; }
    opset& dump() { return _opset; }
    treeop_set& operator+ (treeop_set& rhs) { return connect(rhs); }
    void acquire(treeop_set& ops)
    {
        assert(unary());
        wrapper* m = _opset.front();
        assert(m);
        m->acquire_children(ops.dump());
    }
    treeop_set& sub(treeop_set& ops)
    {
        acquire(ops);
        return *this;
    }
    treeop_set& connect(treeop_set& ops)
    {
        _opset.connect(ops.dump());
        return *this;
    }
    wrapper* detach()
    {
        assert(unary());
        wrapper* w = _opset.front();
        _opset.reset();
        return w;
    }
};

template<class _tree>
class treeop
{
public:
    typedef treeop_helper<_tree> helper;
    typedef typename helper::tree tree;
    typedef typename helper::value value;
    typedef typename helper::wrapper wrapper;
    typedef typename helper::alloc alloc;
    typedef typename helper::const_iterator const_iterator;
    typedef typename helper::iterator iterator;
    typedef typename helper::children children;
    typedef treeop_set<_tree> tropset;
    typedef children opset;

public:
    treeop(_tree& tr) { _tr = &tr; }
    void create(tropset& ops) { _tr->set_root(ops.detach()); }

protected:
    _tree*          _tr;
};

/*
 * Usage:
 * For example,
 * 1.if we wanted to create a tree in the following structure,
 *   { a = 1, b = "hello"; }
 *        { a = 2, b = "bingo"; }
 *            { a = 3, b = "bravo"; }
 *        { a = 6, b = "amigo"; }
 *   and the node was for: struct node { int a; string b; };
 *   We could write it like,
 *   typedef tree<node, _treenode_wrapper<node> > mytree;
 *   def_tr$(node1, mytree, ([&](){ a = 1, b = "hello"; }));
 *   def_tr$(node2, mytree, ([&](){ a = 2, b = "bingo"; }));
 *   def_tr$(node3, mytree, ([&](){ a = 3, b = "bravo"; }));
 *   def_tr$(node4, mytree, ([&](){ a = 6, b = "amigo"; }));
 *   tr$(node1).sub(
 *      tr$(node2).sub(
 *          tr$(node3)
 *          )
 *      + tr$(node4)
 *      );
 * 2.if your compiler doesn't support c++0x style of lambda,
 *   you may modify the def_tr$ micro a little bit, like assignments;
 *   then you could temporally write it like:
 *   typedef tree<node, _treenode_wrapper<node> > mytree;
 *   def_tr$(node1, mytree, (a = 1, b = "hello";));
 *   it works, just maybe unsafe.
 */

#define def_tr$(name, treeop, prototype, assignments) \
struct tag##name \
{ \
    typedef treeop::wrapper wrapper; \
    typedef treeop::value value; \
    typedef treeop::alloc alloc; \
    typedef treeop::tropset tropset; \
    struct constructor: \
        public prototype \
    { \
        constructor() { (assignments)(); } \
    }; \
    tropset create() \
    { \
        wrapper* w = alloc::born(); \
        w->born<constructor>(); \
        return tropset(w); \
    } \
} name;

#define def_tr$_v1(name, treeop, prototype, arg1, assignments) \
struct tag##name \
{ \
    typedef treeop::wrapper wrapper; \
    typedef treeop::value value; \
    typedef treeop::alloc alloc; \
    typedef treeop::tropset tropset; \
    typedef decltype(arg1) type1; \
    type1& ref_arg1; \
    tag##name(type1& arg1): ref_arg1(arg1) {} \
    struct constructor: \
        public prototype \
    { \
        constructor(type1& arg1) { (assignments)(arg1); } \
    }; \
    tropset create() \
    { \
        wrapper* w = alloc::born(); \
        w->born<constructor>(ref_arg1); \
        return tropset(w); \
    } \
} name(arg1);

#define def_tr$_v2(name, treeop, prototype, arg1, arg2, assignments) \
struct tag##name \
{ \
    typedef treeop::wrapper wrapper; \
    typedef treeop::value value; \
    typedef treeop::alloc alloc; \
    typedef treeop::tropset tropset; \
    typedef decltype(arg1) type1; \
    typedef decltype(arg2) type2; \
    type1& ref_arg1; \
    type2& ref_arg2; \
    tag##name(type1& arg1, type2& arg2): ref_arg1(arg1), ref_arg2(arg2) {} \
    struct constructor: \
        public prototype \
    { \
        constructor(type1& arg1, type2& arg2) { (assignments)(arg1, arg2); } \
    }; \
    tropset create() \
    { \
        wrapper* w = alloc::born(); \
        w->born<constructor>(ref_arg1, ref_arg2); \
        return tropset(w); \
    } \
} name(arg1, arg2);

#define def_tr$_v3(name, treeop, prototype, arg1, arg2, arg3, assignments) \
struct tag##name \
{ \
    typedef treeop::wrapper wrapper; \
    typedef treeop::value value; \
    typedef treeop::alloc alloc; \
    typedef treeop::tropset tropset; \
    typedef decltype(arg1) type1; \
    typedef decltype(arg2) type2; \
    typedef decltype(arg3) type3; \
    type1& ref_arg1; \
    type2& ref_arg2; \
    type3& ref_arg3; \
    tag##name(type1& arg1, type2& arg2, type3& arg3): ref_arg1(arg1), ref_arg2(arg2), ref_arg3(arg3) {} \
    struct constructor: \
        public prototype \
    { \
        constructor(type1& arg1, type2& arg2, type3& arg3) { (assignments)(arg1, arg2, arg3); } \
    }; \
    tropset create() \
    { \
        wrapper* w = alloc::born(); \
        w->born<constructor>(ref_arg1, ref_arg2, ref_arg3); \
        return tropset(w); \
    } \
} name(arg1, arg2, arg3);

#define tr$(name) (name.create())

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#pragma once

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

class MaxOptimizationWrapper {
public:
	virtual gsl_vector* getMinState() = 0;
	virtual double getObjectiveVal() = 0;
	virtual bool maximize(Interface* inputs, const gsl_vector* initState, const gsl_vector* initDir, const set<int>& assertConstraints, int minimizeNode, float beta, int level, int id) = 0;
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend