text:
// The first test on 25/3/80 contains no information about
// infection because ewes had not been exposed to the donors at this time
// The first test date is end of March which we consider start of April

// Time is counted in months at the start of each month, so month 0 is 1 April 1980
// A seroconversion observed at the start of month m means the seroconversion
// occurred before the start of month m, eg in month m-1 or before

// Sheep are assumed to be born on the 1st April of each year, ie months 0, 12 and 24
// A sheep's last month is the last month it was tested. We assume that it is immediately removed
// This means that a +ve sheep is not infectious on its last_month
// We assume, however, that a sheep is potentially infected and infectious at most one month prior to
// its first positive month

// Donor sheep 72 is assumed infectious at the start, because its dam tests positive within 3 months of birth
// Donor 22 is assumed to be latent at the start due to it being removed after 3.5 years,
// so it may have just been infected at the start

#include <model.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include "sheep_data.h"

#define P_inv_gamma gsl_cdf_gamma_Pinv // invverse cumulative gamma dist
#define P_inc_gamma gsl_sf_gamma_inc_P // regularised incomplete Gamma function
#define lngamma gsl_sf_lngamma       // ln Gamma function
#define NUMBER_DONORS 2
#define NUMBER_SHEEP 45
#define INTRO_MONTHS 3
#define V0 11

const uint intro_months[] = {0, 12, 24};

typedef struct {
  double a;
  double *lambda;
  double *Lambda;
  double *b;
  double *B;
  double **W4;
  double *W2;
} parameters_t;

double u(int, int, double*);
double lngexpan(double, double, int*);
double ln_inc_gamma_star(double, double);
double w(int, int, int, int, parameters_t*);
double Prob_seroconvert(int, int, parameters_t*);
extern double R0(double, double, double, int, double*);

void Variables(struct variables *V)
{
  int i;
  Var(V,  0, "month");
  Var(V,  1, "expected number infected");
  Var(V,  2, "expected number infectious");
  Var(V,  3, "cumulative seroconversions");
  Var(V,  4, "seroconversion rate");
  Var(V,  5, "expected number seroconverted");
  Var(V,  6, "cumulative field infected");
  Var(V,  7, "cumulative housed infected");
  Var(V,  8, "cumulative poor condition infected");
  Var(V,  9, "cumulative total infected");
  Var(V, 10, "infections");
  for (i = 0; i < INTRO_MONTHS+1; i++)
    Var(V, V0+i, "P_seropos");
}

void Parameters(int trt, int ind, double *p, struct parameterset *Q)
{
  Par(Q,  0, Real, "beta_{field}", Exponential, 0.02);
  Par(Q,  1, Real, "beta_{housed}", Gamma, 2.0, 0.1);
  Par(Q,  2, Real, "a", Uniform, 1.e-3, 100.);
  Par(Q,  3, Real, "k1", Exponential, 24.0);
  Par(Q,  5, Real, "L", Uniform, 0., 60.);
  Par(Q,  7, Int,  "tau", Uniform, 1, 4);
  Par(Q,  8, Real, "beta_{pc}", Gamma, 2.0, 1.0);

  Par(Q, 10, Real, "beta_{housed}/beta_{field}", Der);
  Par(Q, 12, Real, "E[inf_{field}]", Der);
  Par(Q, 13, Real, "E[inf_{housed}]", Der);
  Par(Q, 14, Real, "E[inf]", Der);
  Par(Q, 16, Real, "s_{median,1}", Der);
  Par(Q, 17, Real, "s_{median,2}", Der);
  Par(Q, 18, Real, "s_{95,1}", Der);
  Par(Q, 19, Real, "s_{95,2}", Der);
}

int Model(const uint trt, const uint ind, const uint nmonths, double *p, double **V, TimePoints *Data)
{
  int i, j, n, m;
  double *Einf_field, *Einf_housed;
  double **p_seropos = V, pinf;

  double beta_field = p[0];
  double beta_housed = p[1];
  double beta_pc = p[8];
  double S_a = p[2];
  double S_b1 = S_a/p[3];
  double S_b2 = S_a/p[3];
  double L = p[5];
  int tau = lrint(p[7]);

  double *num_infectious = doublevector(nmonths);
  double *num_infected = doublevector(nmonths);
  // double *num_seroconverted = doublevector(nmonths);
  double *L_cdf = doublevector(nmonths);
  // double *S_cdf = doublevector(nmonths);

  parameters_t pars;
  pars.Lambda = doublevector(nmonths);
  pars.lambda = doublevector(nmonths);
  pars.B = doublevector(nmonths);
  pars.b = doublevector(nmonths);
  pars.W4 = doublematrix(nmonths, nmonths);
  pars.W2 = doublevector(nmonths);
  pars.a = S_a;

  // cumulative distribution of latent period
  for (m = 0; m < nmonths; m++) {
    // L_cdf[m] = gsl_cdf_gamma_P((double)m, L_a, 1./L_b);
    L_cdf[m] = (m < L? 0: 1);
    // S_cdf[m] = gsl_cdf_gamma_P((double)m, pars.a, 1./pars.b);
  }

  // calculate expected number of infectious ewes and force of infection, lambda
  // the two donor ewes and the lamb are infected at start of experiment
  int donors_lamb[] = {0, 1, 44};
  for (i = 0; i < NUMBER_DONORS+1; i++) {
    j = donors_lamb[i];
    for (m = sheep_intro_month[j]; m < sheep_last_month[j]; m++) {
      num_infected[m]++;
      if (j == 0)
        num_infectious[m]++;
      else
        num_infectious[m] += L_cdf[m-sheep_intro_month[j]];
    }
  }

  for (m = 0; m < nmonths; m++) {
    // calculate force of infection
    // ewes are housed for tau months (1 to 4) ending in March
    if (21 <= m && m < 24)
      // transmission rate for poor conditioned housed ewes in winter 1982
      pars.lambda[m] = num_infectious[m]*beta_pc;
    else if (4-tau <= (m+4)%12 && (m+4)%12 <= 3)
      // transmission rate for housed ewes in winter
      pars.lambda[m] = num_infectious[m]*beta_housed;
    else
      // field transmission rate
      pars.lambda[m] = num_infectious[m]*beta_field;

    // cumulative lambda used for prob of uninfected
    if (m > 1)
        pars.Lambda[m] = pars.Lambda[m-1]+pars.lambda[m-1];

    for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
      if (use_sheep[j] && negatives[j] > ind && sheep_intro_month[j] <= m && m < sheep_last_month[j]) {
        // probability this sheep is infected during this month m
        pinf = (1.-exp(-pars.lambda[m]))*u(sheep_intro_month[j], m, pars.Lambda);
        for (n = m+1; n < sheep_last_month[j]; n++) {
            num_infected[n] += pinf;
            num_infectious[n] += pinf*L_cdf[n-1-m];
            // num_seroconverted[i] += pinf*S_cdf[n-1-m];
        }
      }
  }

  // monthly seroconversion rates
  for (m = 0; m < nmonths; m++)
    if (m == 22)
      // increased seroconversion rate in winter 1982 due to poor condition and
      pars.b[m] = S_b2;
    else
      pars.b[m] = S_b1;
  // cumulative seroconversion rate
  for (m = 1; m < nmonths; m++)
    pars.B[m] = pars.B[m-1] + pars.b[m-1];

  // W4 and W2 per month
  for (m = 0; m < nmonths; m++) {
    for (n = 0; n < m; n++)
      pars.W4[n][m] = w(n, m, n+1, m+1, &pars)
                    - w(n, m, n+1, m,   &pars)
                    - w(n, m, n,   m+1, &pars)
                    + w(n, m, n,   m,   &pars);
    pars.W2[m] = w(m, m, m, m+1, &pars);
  }

  // calculate the probability that a randomly selected sheep seroconverts in month m
  for (j = 0; j < INTRO_MONTHS; j++)
    for (m = intro_months[j]; m < nmonths; m++)
      p_seropos[m][V0+j] = Prob_seroconvert(intro_months[j], m, &pars);

  // 3 sheep were born in month 12, were removed in month 20 and retested in month 36, all were negative
  j = 13;
  // save lambda
  double *lambda_s = doublevector(nmonths);
  double *Lambda_s = doublevector(nmonths);
  for (m = sheep_last_month[j]; m < sheep_last_test_month[j]; m++) {
    lambda_s[m] = pars.lambda[m];
    Lambda_s[m] = pars.Lambda[m];
  }

  // set lambda=0 for the months from removal (last_month) until testing and recalculate Lambda
  for (m = sheep_last_month[j]; m < sheep_last_test_month[j]; m++)
    pars.lambda[m] = 0;
  for (m = sheep_last_month[j]+1; m < sheep_last_test_month[j]; m++)
    pars.Lambda[m] = pars.Lambda[m-1];

  // recalculate prob of seroconverting for these 3 sheep
  for (m = sheep_last_month[j]; m < sheep_last_test_month[j]; m++) {
    for (n = sheep_last_month[j]; n < m; n++)
      pars.W4[n][m] = w(n, m, n+1, m+1, &pars)
                    - w(n, m, n+1, m,   &pars)
                    - w(n, m, n,   m+1, &pars)
                    + w(n, m, n,   m,   &pars);
    pars.W2[m] = w(m, m, m, m+1, &pars);
  }
  for (m = sheep_intro_month[j]; m < sheep_last_test_month[j]; m++)
    p_seropos[m][V0+sheep_var[j]] = Prob_seroconvert(sheep_intro_month[j], m, &pars);

  // restore lambda
  for (m = sheep_last_month[j]; m < sheep_last_test_month[j]; m++) {
    pars.lambda[m] = lambda_s[m];
    pars.Lambda[m] = Lambda_s[m];
  }
  free(lambda_s);
  free(Lambda_s);

  /******************** OUTPUT ***********************/
  if (Data->mode == OUTPUT) {
    for (m = 0; m < nmonths; m++)
      V[m][0] = m;

    // expected number of infected sheep in month m
    for (m = 0; m < nmonths-1; m++)
      V[m][1] = num_infected[m];
    V[nmonths-1][1] = V[nmonths-2][1];

    // expected number of infectious sheep in month m
    for (m = 0; m < nmonths-1; m++)
      V[m][2] = num_infectious[m];
    V[nmonths-1][2] = V[nmonths-2][2];

    // expected number of seroconverted sheep in month m
    // for (m = 0; m < nmonths-1; m++)
    //   V[m][5] = num_seroconverted[m];
    // V[nmonths-1][5] = V[nmonths-2][5];

    // expected number of seroconversions by the end of each month m
    // (which is the start of month m+1)
    // done this way because earlier observed seroconversions are recorded at start of month
    V[0][4] = 0;
    for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
      if (use_sheep[j] && negatives[j] > ind)
        for (m = sheep_intro_month[j]; m < sheep_last_month[j]; m++)
          V[m+1][4] += p_seropos[m][V0+sheep_var[j]];

    // cumulative expected number of seroconversions by end of month m
    // (which is the start of month m+1)
    V[0][3] = 0;
    for (m = 0; m <nmonths-1; m++)
      V[m+1][3] = V[m][3]+V[m+1][4];

    // expected number of infection events in each month m partitioned by route
    Einf_field = doublevector(nmonths);
    Einf_housed = doublevector(nmonths);
    for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
      if (use_sheep[j] && negatives[j] > ind) {
        for (m = sheep_intro_month[j]; m < sheep_last_month[j]; m++)
          if (4-tau <= (m+4)%12 && (m+4)%12 <= 3)
            Einf_housed[m] += (1.-exp(-pars.lambda[m]))*u(sheep_intro_month[j], m , pars.Lambda);
          else
            Einf_field[m] += (1.-exp(-pars.lambda[m]))*u(sheep_intro_month[j], m, pars.Lambda);
      }
    for (m = 0; m < nmonths; m++)
      V[m][10] = Einf_field[m]+Einf_housed[m];

    for (m = 1; m < nmonths; m++) {
      Einf_field[m] += Einf_field[m-1];
      Einf_housed[m] += Einf_housed[m-1];
    }
    for (m = 0; m < nmonths; m++) {
      V[m][6] = Einf_field[m];
      V[m][7] = Einf_housed[m];
      V[m][9] = V[m][6]+V[m][7];
    }
    free(Einf_field);
    free(Einf_housed);
  }

  free(L_cdf);
  // free(S_cdf);
  free(num_infected);
  free(num_infectious);
  // free(num_seroconverted);
  free(pars.lambda);
  free(pars.Lambda);
  free(pars.b);
  free(pars.B);
  free(pars.W4[0]);
  free(pars.W4);
  free(pars.W2);
  return SUCCESS;
}

double P_pos(int j, int s1, int s2, double **V)
{
  int m;
  double sum = 0;
  for (m = s1; m < s2; m++)
    sum += V[m][V0+sheep_var[j]];
  return sum;
}

double logLikelihood(const uint trt, const uint ind, const double *p, double **V, const TimePoints *Data)
{
  int j, pos_idx;
  double lnL = 0;

  for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
    if (use_sheep[j] && negatives[j] > ind) {
      pos_idx = sheep_pos_index[j];
      if (pos_idx > 0)
        // ewe is first seropositive at start of month m, therefore it seroconverted
        // sometime in the months from the last test to the month prior to seroconverting
        lnL += log(P_pos(j, lrint(Data->t[pos_idx-1][0]), sheep_pos_month[j], V));
      else
        // ewe never seroconverted
        lnL += log(1. - P_pos(j, sheep_intro_month[j], sheep_last_test_month[j], V));
    }
  return lnL;
}

void OutputModel(int trt, int ind, double *Output, double *V)
{
  int i;
  for (i = 1; i < V0+INTRO_MONTHS+1; i++)
    Output[i] = V[i];
}

void OutputData(int trt, int ind, double *Output, double *Var, double *Data, uint *value, int index)
{
  Output[3] = Data[47];
  Output[4] = Data[48];
  Output[5] = Data[49];
  // remove dam to lamb infection from data to match model output
  if (lrint(Data[0]) == 3)
    Output[4] = 0;
  if (lrint(Data[0]) >= 3)
    Output[3] -= 1;
}

int f0(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 2 month housed
  int i, n = 2;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = 100*i+1;
    y[i] = R0(x[i], 2., 0.1, 5, p);
  }
  return n;
}

int f1(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 0 months housed
  int i, n = 2;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = 100*i+1;
    y[i] = R0(x[i], 0., 0.1, 5, p);
  }
  return n;
}

int f2(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // R0 for 25 sheep housed for different durations
  int i, n = 14;
  if (x == NULL) return n;
  for (i = 0; i < n; i++) {
    x[i] = i;
    y[i] = R0(25., x[i]*12./365., 0.1, 5, p);
  }
  return n;
}

int f3(int trt, int ind, double *x, double *y, double *p, TimePoints *TP)
{
  // CDF of normal seroconversion period
  int i, n = 1000;
  if (x == NULL) return n;
  double S_a = p[2];
  double S_k1 = p[3];

  for (i = 0; i < n; i++) {
    x[i] = (double)i/20.;
    // y[i] = gsl_ran_gamma_pdf(x[i], a, 1./b);
    y[i] = gsl_cdf_gamma_P(x[i], S_a, S_k1);
  }
  return n;
}

void function_list(functions_t *F)
{
  F->n_func = 3;
  F->function_list = (function_ptr*) malloc(F->n_func*sizeof(function_ptr));
  F->function_list[0] = f0;
  F->function_list[1] = f1;
  F->function_list[2] = f2;
  // F->function_list[3] = f3;
}

double u(int m0, int m, double *Lambda)
{
  // The probability a sheep is uninfected at the start of month m
  // given it was introduced into the flock in month m0
  if (m == m0)
    return 1;
  else
    return exp(-Lambda[m]+Lambda[m0]);
}

double lngseries(double a, double x, int *ierr)
{
  double giant = HUGE_VAL/1000., eps = 1e-15;
  double t = 1./a, v = t;
  double p, q, lnigam;
  int k = 0;
  *ierr = 0;

  while ((fabs(t/v) > eps) && *ierr == 0) {
    p = (a+k)/(a+k+1);
    q = k+1;
    t *= -x*p/q;
    v += t;
    k += 1;
    if (t > giant)
      *ierr = 1;
  }
  if (*ierr == 0)
    if (lngamma(a) < log(giant))
      lnigam = log(v)-lngamma(a);
    else {
      lnigam = 0;
      *ierr = 1;
    }
  else {
    lnigam = 0;
    *ierr = 1;
  }
  return lnigam;
}

double lngexpan(double a, double x, int *ierr)
{
  double giant = HUGE_VAL/1000., eps = 1e-15;
  double t = 1, v = t;
  double p, lnigam;
  int k = 0;
  *ierr = 0;

  if (x > log(giant)) {
      *ierr = 1;
      return 0;
  }
  while ((fabs(t/v) > eps) && *ierr == 0) {
    p = 1-a+k;
    t *= p/x;
    v += t;
    k += 1;
    if (t > giant)
      *ierr = 1;
  }
  if (*ierr == 0)
    if (lngamma(a) < log(giant))
      lnigam = log(v)+x-log(x)-lngamma(a);
    else {
      lnigam = 0;
      *ierr = 1;
    }
  else {
    lnigam = 0;
    *ierr = 1;
  }
  return lnigam;
}

double ln_inc_gamma_star(double a, double z)
{
  int ierr;
  double lnigam;

  if (z < -50) {
    lnigam = lngexpan(a, -z, &ierr);
    if (ierr != 0)
      printf("error1\n");
  }
  else {
    lnigam = lngseries(a, z, &ierr);
    if (ierr != 0)
      printf("error2\n");
  }
  return lnigam;
}

double w(int n, int m, int t, int s, parameters_t *p)
{
  double a = p->a;
  double l = p->lambda[n];
  double b = p->b[n];
  double r = (b-l)/b;
  double Z, z, v;

  if (m-n <= 1 && t == m && s == m)
    return 0;

  Z = p->B[m] - p->B[n] + b*(n-t) + p->b[m]*(s-m);
  if (Z == 0)
    return 0;

  z = r*Z;
  if (z <= 0)
    v = exp(a*log(Z) - l*Z/b + ln_inc_gamma_star(a, z));
  else
    v = exp(-a*log(r) - l*Z/b + log(P_inc_gamma(a, z)));

  return exp(-l*(t-n))*(v - P_inc_gamma(a, Z));
}

double Prob_seroconvert(int m0, int m, parameters_t *p)
{
  int n;
  double P_m = 0;

  for (n = m0; n < m; n++)
      P_m += u(m0, n, p->Lambda)*p->W4[n][m];
  P_m -= u(m0, m, p->Lambda)*p->W2[m];

  return min(max(0, P_m), 1);
}

void Residual(int trt, int ind, double *R, double *V, double *p, double *Data, uint *value)
{
  R[3] = V[3]-Data[47];
  R[4] = V[4]-Data[48];
}

void DerivedParameters(int trt, int ind, int nmonths, double *p, double **V, TimePoints *Data)
{
  // int i, j, m;
  // double pinf;
  double beta_field = p[0];
  double beta_housed = p[1];
  double S_a = p[2];
  double S_b1 = S_a/p[3];
  double S_b2 = S_a/p[3];
  // double L_mu = p[5];
  // double L_sd = p[6];
  // double *num_infectious = doublevector(nmonths);
  // double *Lambda = doublevector(nmonths);
  // double *lambda = doublevector(nmonths);
  // double *B = doublevector(nmonths);
  // double *b = doublevector(nmonths);
  // double *L_cdf = doublevector(nmonths);

  p[10] = beta_housed/beta_field;
  p[16] = P_inv_gamma(0.5, S_a, 1./S_b1);
  p[17] = P_inv_gamma(0.5, S_a, 1./S_b2);
  p[18] = P_inv_gamma(0.95, S_a, 1./S_b1);
  p[19] = P_inv_gamma(0.95, S_a, 1./S_b2);
  return;

  // p[6] = gsl_cdf_exponential_Pinv(0.5, mu);
  // p[22] = gsl_cdf_exponential_Pinv(0.5, mu2);
  // p[15] = gsl_cdf_exponential_Pinv(0.025, mu);
  // p[16] = gsl_cdf_exponential_Pinv(0.975, mu);
  // p[20] = 6./16.*(1.-exp(-(0.75*beta_field)*(fmax(0, 4*12-L_mu))));
  // p[23] = 6./16.*(1.-exp(-(0.25*beta_housed)*(fmax(0, 4*12-L_mu))));

  // for (m = 0; m < nmonths; m++)
  //    L_cdf[m] = (m < L_mu? 0: 1);

  // // calculate expected number of infectious ewes and force of infection, lambda
  // // the two donor sheep are infectious at start of experiment
  // for (j = 0; j < 2; j++)
  //   for (m = 0; m < sheep_last_month[j]; m++)
  //     num_infectious[m]++;
  // // ewe 44 - infected by dam
  // for (m = sheep_intro_month[44]; m < sheep_last_month[44]; m++)
  //   num_infectious[m] += L_cdf[m-sheep_intro_month[44]];

  // for (m = 0; m < nmonths; m++) {
  //   // calculate force of infection
  //   if ((m+3)%12 < 3)
  //     // transmission rate for housed ewes (in winter Jan-Mar)
  //     lambda[m] = num_infectious[m]*beta_housed;
  //   else
  //     // field transmission rate
  //     lambda[m] = num_infectious[m]*beta_field;

  //   // cumulative lambda for prob of uninfected
  //   if (m > 1)
  //       Lambda[m] = Lambda[m-1]+lambda[m-1];

  //   for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
  //     if (use_sheep[j] && sheep_intro_month[j] <= m && m < sheep_last_month[j]) {
  //       pinf = (1.-exp(-lambda[m]))*u(sheep_intro_month[j], m, Lambda);
  //       for (i = m; i < sheep_last_month[j]; i++)
  //           num_infectious[i] += pinf*L_cdf[i-m];
  //     }
  // }

  // p[12] = 0;
  // p[13] = 0;
  // for (j = NUMBER_DONORS; j < NUMBER_SHEEP; j++)
  //   if (use_sheep[j]) {
  //     for (m = sheep_intro_month[j]; m < sheep_last_month[j]; m++)
  //       if ((m+3)%12 < 3)
  //         p[13] += (1.-exp(-lambda[m]))*u(sheep_intro_month[j], m, Lambda);
  //       else
  //         p[12] += (1.-exp(-lambda[m]))*u(sheep_intro_month[j], m, Lambda);
  //   }
  // p[14] = p[12]+p[13];

  // free(lambda);
  // free(Lambda);
  // free(B);
  // free(b);
  // free(num_infectious);
  // free(L_cdf);
}

void WAIC(int trt, int ind, double **lnL, double **V, double *p, const TimePoints *Data) {}
void PredictData(int trt, int ind, double *Output, double *V, double *p, gsl_rng *stream) {}
void SimulateData(int trt, int ind, double *Output, double *V, double *p, double *Data, uint *value, gsl_rng *stream) {}
double timestep(void) {return 1;}
void GlobalParameters() {}
double UserPDF(double x) {return 0;}
void HyperParameters(Treatments T, Hyperparameters *H) {}
void SaturatedModel(int trt, int ind, double **V, double *p, const TimePoints *Data) {}

max_stars_repo_name: nicksavill/maedi-visna-epidemiology

max_issues_repo_name: nicksavill/maedi-visna-epidemiology

max_forks_repo_name: nicksavill/maedi-visna-epidemiology


text:
// SPDX-License-Identifier: MIT
// The MIT License (MIT)
//
// Copyright (c) 2014-2018, Institute for Software & Systems Engineering
// Copyright (c) 2018-2019, Johannes Leupolz
//
// 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.

#ifndef PEMC_EXECUTABLE_MODEL_TEMPORARY_STATE_STORAGE_H_
#define PEMC_EXECUTABLE_MODEL_TEMPORARY_STATE_STORAGE_H_

#include <vector>
#include <gsl/span>
#include <cstdint>
#include <atomic>
#include <boost/align/aligned_allocator.hpp>

#include "pemc/basic/tsc_index.h"
#include "pemc/basic/probability.h"
#include "pemc/basic/label.h"
#include "pemc/basic/model_capacity.h"
#include "pemc/basic/raw_memory.h"
#include "pemc/formula/formula.h"

namespace pemc {

  ///   We store states in a contiguous array, indexed by a continuous variable.
  ///   The TemporaryStateStorage is not thread safe.
  class TemporaryStateStorage {
  private:

      // The length in bytes of a state vector required for the analysis model.
      int32_t modelStateVectorSize = 0;
      // Extra bytes in state vector for preStateStorage modifiers.
      int32_t preStateStorageModifierStateVectorSize = 0;
      // The length in bytes of the state vector of the analysis model with the extra bytes
      // required for the preStateStorage modifiers
      int32_t stateVectorSize = 0;

      // The number of saved states
      StateIndex savedStates = 0;
  	  // The number of states that can be cached and the number of reserved states.
      StateIndex totalCapacity;

      // the memory that contains the serialized states
      std::vector<gsl::byte> stateMemory;

      void resizeStateBuffer();

  public:
      TemporaryStateStorage(StateIndex _capacity);

      StateIndex getNumberOfSavedStates();

      gsl::span<gsl::byte> operator [](size_t idx);

      StateIndex getFreshStateIndex();

      void setStateVectorSize(int32_t _modelStateVectorSize, int32_t _preStateStorageModifierStateVectorSize);

      void clear();
  };

}

#endif  // PEMC_EXECUTABLE_MODEL_TEMPORARY_STATE_STORAGE_H_

ext: h

max_stars_repo_name: joleuger/pemc

max_issues_repo_name: joleuger/pemc

max_forks_repo_name: joleuger/pemc


text:
/* movstat/medacc.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * Original copyright notice:
 * Copyright (c) 2011 ashelly.myopenid.com under <http://www.opensource.org/licenses/mit-license>
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>

#define ItemLess(a,b)  ((a)<(b))
#define ItemMean(a,b)  (((a)+(b))/2)

#define minCt(m) (((m)->ct-1)/2) /* count of items in minheap */
#define maxCt(m) (((m)->ct)/2)   /* count of items in maxheap */

typedef double medacc_type_t;

typedef struct
{
  int n;                /* window size */
  int idx;              /* position in circular queue */
  int ct;               /* count of items in queue */
  medacc_type_t *data;  /* circular queue of values, size k */
  int *pos;             /* index into `heap` for each value, size 2*k */
  int *heap;            /* max/median/min heap holding indices into `data` */
} medacc_state_t;

static size_t medacc_size(const size_t n);
static int medacc_init(const size_t n, void * vstate);
static int medacc_insert(const medacc_type_t x, void * vstate);
static int medacc_delete(void * vstate);
static int medacc_get(void * params, medacc_type_t * result, const void * vstate);

static int mmless(const medacc_state_t * state, const int i, const int j);
static int mmexchange(medacc_state_t * state, const int i, const int j);
static int mmCmpExch(medacc_state_t * state, const int i, const int j);
static void minSortDown(medacc_state_t * state, int i);
static void maxSortDown(medacc_state_t * state, int i);
static int minSortUp(medacc_state_t * state, int i);
static int maxSortUp(medacc_state_t * state, int i);

static size_t
medacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(medacc_state_t);
  size += n * sizeof(medacc_type_t);
  size += 2 * n * sizeof(int);

  return size;
}

static int
medacc_init(const size_t n, void * vstate)
{
  medacc_state_t * state = (medacc_state_t *) vstate;
  int k = (int) n;

  state->n = n;
  state->ct = 0;
  state->idx = 0;

  state->data = (medacc_type_t *) ((unsigned char *) vstate + sizeof(medacc_state_t));
  state->pos = (int *) ((unsigned char *) state->data + n * sizeof(medacc_type_t));
  state->heap = state->pos + n + (n/2); /* points to middle of storage */

  /* set up initial heap fill pattern: median,max,min,max,... */
  while (k--)
    {
      state->pos[k] = ((k + 1)/2) * ((k & 1) ? -1 : 1);
      state->heap[state->pos[k]] = k;
    }

  return GSL_SUCCESS;
}

static int
medacc_insert(const medacc_type_t x, void * vstate)
{
  medacc_state_t * state = (medacc_state_t *) vstate;
  int isNew = (state->ct < (int) state->n);
  int p = state->pos[state->idx];
  medacc_type_t old = state->data[state->idx];

  state->data[state->idx] = x;
  state->idx = (state->idx + 1) % state->n;
  state->ct += isNew;

  if (p > 0)       /* new item is in minHeap */
    {
      if (!isNew && ItemLess(old, x))
        minSortDown(state, p * 2);
      else if (minSortUp(state, p))
        maxSortDown(state, -1);
    }
  else if (p < 0)  /* new item is in maxHeap */
    {
      if (!isNew && ItemLess(x, old))
        maxSortDown(state, p * 2);
      else if (maxSortUp(state, p))
        minSortDown(state, 1);
    }
  else             /* new item is at median */
    {
      if (maxCt(state))
        maxSortDown(state, -1);

      if (minCt(state))
        minSortDown(state, 1);
    }

  return GSL_SUCCESS;
}

static int
medacc_delete(void * vstate)
{
  medacc_state_t * state = (medacc_state_t *) vstate;

  if (state->ct > 0)
    {
      int p = state->pos[(state->idx - state->ct + state->n) % state->n];

      if (p > 0)        /* oldest item is in minHeap */
        {
          mmexchange(state, p, minCt(state));
          --(state->ct);
          minSortDown(state, 2 * p);
        }
      else if (p < 0)   /* oldest item is in maxHeap */
        {
          mmexchange(state, p, maxCt(state));
          --(state->ct);
          maxSortDown(state, 2 * p);
        }
      else if (p == 0)  /* oldest item is at median */
        {
        }
    }

  return GSL_SUCCESS;
}

/* returns median (or average of 2 when item count is even) */
static int
medacc_get(void * params, medacc_type_t * result, const void * vstate)
{
  const medacc_state_t * state = (const medacc_state_t *) vstate;
  medacc_type_t median = state->data[state->heap[0]];

  (void) params;

  if ((state->ct & 1) == 0)
    median = ItemMean(median, state->data[state->heap[-1]]);

  *result = median;

  return GSL_SUCCESS;
}

/* returns 1 if heap[i] < heap[j] */
static int
mmless(const medacc_state_t * state, const int i, const int j)
{
  return ItemLess(state->data[state->heap[i]], state->data[state->heap[j]]);
}

/* swaps items i and j in heap, maintains indexes */
static int
mmexchange(medacc_state_t * state, const int i, const int j)
{
  int t = state->heap[i];
  state->heap[i] = state->heap[j];
  state->heap[j] = t;
  state->pos[state->heap[i]] = i;
  state->pos[state->heap[j]] = j;
  return 1;
}

/* swaps items i and j if i < j; returns true if swapped */
static int
mmCmpExch(medacc_state_t * state, const int i, const int j)
{
  return (mmless(state, i, j) && mmexchange(state , i, j));
}

/* maintains minheap property for all items below i/2. */
static void
minSortDown(medacc_state_t * state, int i)
{
  for (; i <= minCt(state); i *= 2)
    {
      if (i > 1 && i < minCt(state) && mmless(state, i + 1, i))
        ++i;

      if (!mmCmpExch(state, i, i / 2))
        break;
   }
}

/* maintains maxheap property for all items below i/2. (negative indexes) */
static void
maxSortDown(medacc_state_t * state, int i)
{
  for (; i >= -maxCt(state); i *= 2)
    {
      if (i < -1 && i > -maxCt(state) && mmless(state, i, i - 1))
        --i;

      if (!mmCmpExch(state, i / 2, i))
        break;
   }
}

/* maintains minheap property for all items above i, including median
   returns true if median changed */
static int
minSortUp(medacc_state_t * state, int i)
{
  while (i > 0 && mmCmpExch(state, i, i / 2))
    i /= 2;

  return (i == 0);
}

/* maintains maxheap property for all items above i, including median
   returns true if median changed */
static int
maxSortUp(medacc_state_t * state, int i)
{
  while (i<0 && mmCmpExch(state, i / 2, i))
    i /= 2;

  return (i == 0);
}

static const gsl_movstat_accum median_accum_type =
{
  medacc_size,
  medacc_init,
  medacc_insert,
  NULL, /* XXX FIXME */
  medacc_get
};

const gsl_movstat_accum *gsl_movstat_accum_median = &median_accum_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* examples/C/ssmfe/precond_expert.c */
/* Laplacian on a square grid (using SPRAL_SSMFE_EXPERT routines) */
#include "spral.h"
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>
#include <math.h>

/* Header that implements Laplacian and preconditioners */
#include "laplace2d.h"
#include "ldltf.h"

int test_core(void);
int test_core_d(int, int);
int test_core_z(int, int);
int test_expert(void);
int test_expert_d(int);
int test_expert_z(int);
int test_ssmfe(void);
int test_ssmfe_d(int);
int test_ssmfe_z(int);

void zdcopy( int n, double complex* x, double* y );
void dzcopy( int n, double* y, double complex* x );
void vector_operations_d
  ( struct spral_ssmfe_rcid rci,
    int n, int m, int kw, int ind[m],
    double W[kw][m][n],
    double rr[3][2*m][2*m],
    double U[m]);

void vector_operations_z
  ( struct spral_ssmfe_rciz rci,
    int n, int m, int kw, int ind[m],
    double complex W[kw][m][n],
    double complex rr[3][2*m][2*m],
    double complex U[m]);

int main(void) {

  int errors = 0;
  int err;

  fprintf(stdout, "testing ssmfe_core...\n");
  err = test_core();
  errors += err;
  fprintf(stdout, "%d errors\n", err);

  fprintf(stdout, "testing ssmfe_expert...\n");
  err = test_expert();
  errors += err;
  fprintf(stdout, "%d errors\n", err);

  fprintf(stdout, "testing ssmfe...\n");
  err = test_ssmfe();
  errors += err;
  fprintf(stdout, "%d errors\n", err);

  fprintf(stdout, "=============================\n");
  fprintf(stdout, "Total number of errors = %d\n", errors);

  return errors;
}

int test_core(void) {

  int errors = 0;
  int err;

  fprintf(stdout, "testing standard_double...");
  err = test_core_d(0, 0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double...");
  err = test_core_d(1, 0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing largest_double...");
  err = test_core_d(0, 1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_double_complex...");
  err = test_core_z(0, 0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double_complex...");
  err = test_core_z(1, 0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing largest_double_complex...");
  err = test_core_z(0, 1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  return errors;

}

int test_core_d(int problem, int largest) {
   const int ngrid = 20;          /* grid points along each side */
   const int n     = ngrid*ngrid; /* problem size */
   const int nep   = 5;           /* eigenpairs wanted */
   const int m     = 3;           /* dimension of the iterated subspace */

   const int ncon_x = 5;
   const double lambda_x[] = { /* expected eigenvalues */
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01
   };
   const double lambda_lx[] = { /* expected eigenvalues in 'largest' mode */
      -7.9553233,
      -7.8888073,
      -7.8888073,
      -7.8222912,
      -7.7795994
   };
   int errors;

   int state = SPRAL_RANDOM_INITIAL_SEED; /* PRNG state */

   int ind[m];                    /* permutation index */
   double lambda[n];              /* eigenvalues storage */
   double X[n][n];                /* eigenvectors storage */
   /* Work arrays */
   double lmd[m];
   double rr[3][2*m][2*m];
   double W[8][m][n];
   double U[m][n];
   double V[m];

   /* Derived types */
   struct spral_ssmfe_rcid rci;            /* reverse communication data */
   struct spral_ssmfe_core_options options;/* options */
   void *keep;                             /* private data */
   struct spral_ssmfe_inform inform;       /* information */

   /* Initialize options to default values */
   spral_ssmfe_core_default_options(&options);

   /* Initialize W to lin indep vectors by randomizing */
   for(int i=0; i<n; i++)
   for(int j=0; j<m; j++)
      W[0][j][i] = spral_random_real(&state, true);

   int ncon = 0;   /* number of converged eigenpairs */
   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      if ( largest != 0 )
        spral_ssmfe_largest_double(&rci, problem, nep, m, lmd,
           &rr[0][0][0], ind, &keep, &options, &inform);
      else
        spral_ssmfe_double(&rci, problem, nep, 0, m, lmd,
           &rr[0][0][0], ind, &keep, &options, &inform);
      switch ( rci.job ) {
      case 1:
         apply_laplacian(
            ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
            );
         if ( largest != 0 )
            cblas_dscal( n*rci.nx, -1.0, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 2:
         if ( largest != 0 )
            cblas_dcopy(
               n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
               );
         else
            apply_gauss_seidel_step (
               ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
               );
         break;
      case 3:
         cblas_dcopy(
            n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
            );
         break;
      case 4:
         for ( int i = 0; i < m; i++ ) {
            if ( inform.converged[i] != 0 )
               continue;
            if ( inform.err_X[i] > 0 && inform.err_X[i] < 1e-6 )
               inform.converged[i] = 1;
         }
         break;
      case 5:
         if ( rci.i < 0 ) continue;
         for(int k=0; k<rci.nx; k++) {
           int j = ncon + k;
           lambda[j] = lmd[rci.jx + k];
           cblas_dcopy( n, &W[0][rci.jx + k][0], 1, &X[j][0], 1 );
         }
         ncon += rci.nx;
         if ( ncon >= nep || inform.iteration > 300 )
            goto finished;
         break;
      case 11:
      case 12:
      case 13:
      case 14:
      case 15:
      case 16:
      case 17:
        vector_operations_d( rci, n, m, 8, ind, W, rr, V );
         break;
      case 21: // Fall through to 22
      case 22:
         if( ncon > 0 ) {
            cblas_dgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, ncon, rci.nx, n,
               1.0, &X[0][0], n, &W[rci.ky][rci.jy][0], n, 0.0, &U[0][0], n
               );
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
               -1.0, &X[0][0], n, &U[0][0], n, 1.0, &W[rci.kx][rci.jx][0], n
               );
            if ( problem != 0 )
              cblas_dgemm(
                 CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
                 -1.0, &X[0][0], n, &U[0][0], n, 1.0, &W[rci.ky][rci.jy][0], n
                 );
         }
         break;
      case 999:
         if( rci.k == 0 ) {
            if( rci.jx > 1 ) {
               for(int j=0; j<rci.jx; j++)
               for(int i=0; i<n; i++)
                  W[0][j][i] = spral_random_real(&state, true);
            }
         }
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( ncon != ncon_x ) errors++;
   if ( largest != 0 ) {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_lx[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   spral_ssmfe_core_free(&keep, &inform);

   return errors;
}

int test_core_z(int problem, int largest) {
   const int ngrid = 20;          /* grid points along each side */
   const int n     = ngrid*ngrid; /* problem size */
   const int nep   = 5;           /* eigenpairs wanted */
   const int m     = 3;           /* dimension of the iterated subspace */

   const int ncon_x = 5;       /* expected number of converged eigenpairs */
   const double lambda_x[] = { /* expected eigenvalues */
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01
   };
   const double lambda_lx[] = { /* expected eigenvalues in 'largest' mode */
      -7.9553233,
      -7.8888073,
      -7.8888073,
      -7.8222912,
      -7.7795994
   };

   const double complex ZERO = 0.0;
   const double complex ONE = 1.0;
   const double complex MINUS_ONE = -1.0;

   int errors;

   int state = SPRAL_RANDOM_INITIAL_SEED; /* PRNG state */

   int ind[m];                    /* permutation index */
   double lambda[n];              /* eigenvalues storage */
   double complex X[n][n];                /* eigenvectors storage */
   /* Work arrays */
   double lmd[m];
   double complex rr[3][2*m][2*m];
   double complex W[8][m][n];
   double complex U[m][n];
   double complex V[m];

   /* Derived types */
   struct spral_ssmfe_rciz rci;            /* reverse communication data */
   struct spral_ssmfe_core_options options;/* options */
   void *keep;                             /* private data */
   struct spral_ssmfe_inform inform;       /* information */

   /* Initialize options to default values */
   spral_ssmfe_core_default_options(&options);

   /* Initialize W to lin indep vectors by randomizing */
   for(int i=0; i<n; i++)
   for(int j=0; j<m; j++)
      W[0][j][i] = spral_random_real(&state, true);

   int ncon = 0;   /* number of converged eigenpairs */
   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      if ( largest != 0 )
        spral_ssmfe_largest_double_complex(&rci, problem, nep, m, lmd,
           &rr[0][0][0], ind, &keep, &options, &inform);
      else
        spral_ssmfe_double_complex(&rci, problem, nep, 0, m, lmd,
           &rr[0][0][0], ind, &keep, &options, &inform);
      switch ( rci.job ) {
      case 1:
         apply_laplacian_z(
            ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
            );
         if ( largest != 0 )
            cblas_zscal( n*rci.nx, &MINUS_ONE, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 2:
         if ( largest != 0 )
            cblas_zcopy(
               n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
               );
         else
            apply_gauss_seidel_step_z (
               ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
               );
         break;
      case 3:
         cblas_zcopy(
            n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
            );
         break;
      case 4:
         for ( int i = 0; i < m; i++ ) {
            if ( inform.converged[i] != 0 )
               continue;
            if ( inform.err_X[i] > 0 && inform.err_X[i] < 1e-6 )
               inform.converged[i] = 1;
         }
         break;
      case 5:
         if ( rci.i < 0 ) continue;
         for(int k=0; k<rci.nx; k++) {
           int j = ncon + k;
           lambda[j] = lmd[rci.jx + k];
           cblas_zcopy( n, &W[0][rci.jx + k][0], 1, &X[j][0], 1 );
         }
         ncon += rci.nx;
         if ( ncon >= nep || inform.iteration > 300 )
            goto finished;
         break;
      case 11:
      case 12:
      case 13:
      case 14:
      case 15:
      case 16:
      case 17:
        vector_operations_z( rci, n, m, 8, ind, W, rr, V );
         break;
      case 21: // Fall through to 22
      case 22:
         if( ncon > 0 ) {
            cblas_zgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, ncon, rci.nx, n,
               &ONE, &X[0][0], n, &W[rci.ky][rci.jy][0], n, &ZERO, &U[0][0], n
               );
            cblas_zgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
               &MINUS_ONE, &X[0][0], n, &U[0][0], n, &ONE,
               &W[rci.kx][rci.jx][0], n
               );
            if ( problem != 0 )
              cblas_zgemm(
                 CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
                 &MINUS_ONE, &X[0][0], n, &U[0][0], n, &ONE,
                 &W[rci.ky][rci.jy][0], n
                 );
         }
         break;
      case 999:
         if( rci.k == 0 ) {
            if( rci.jx > 1 ) {
               for(int j=0; j<rci.jx; j++)
               for(int i=0; i<n; i++)
                  W[0][j][i] = spral_random_real(&state, true);
            }
         }
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( ncon != ncon_x ) errors++;
   if ( largest != 0 ) {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_lx[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   spral_ssmfe_core_free(&keep, &inform);

   return errors;
}

int test_expert(void) {

  int errors = 0;
  int err;

  fprintf(stdout, "testing standard_double...");
  err = test_expert_d(0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double...");
  err = test_expert_d(1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_shift_double...");
  err = test_expert_d(2);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_shift_double...");
  err = test_expert_d(3);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing buckling_double...");
  err = test_expert_d(4);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_double_complex...");
  err = test_expert_z(0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double_complex...");
  err = test_expert_z(1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_shift_double_complex...");
  err = test_expert_z(2);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_shift_double_complex...");
  err = test_expert_z(3);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing buckling_double_complex...");
  err = test_expert_z(4);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  return errors;
}

int test_expert_d(int problem) {

   const int ngrid = 20;          /* grid points along each side */
   const int n     = ngrid*ngrid; /* problem size */
   const int m     = 3;           /* dimension of the iterated subspace */

   const int ncon_x = 6;
   const double lambda_x[] = {
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01,
      2.2040061e-01
   };
   const double lambda_six[] = {
      9.5108266e-01,
      9.5108266e-01,
      8.8141871e-01,
      8.8141871e-01,
      8.4187478e-01,
      8.4187478e-01
   };

   int errors;
   int nep = 5;

   int state = SPRAL_RANDOM_INITIAL_SEED; /* PRNG state */

   int ind[m];                 /* permutation index */
   int ipiv[n];               /* LDLT pivot index */
   double sigma;              /* shift */
   double A[n][n];            /* matrix */
   double LDLT[n][n];         /* factors */
   double lambda[n];          /* eigenvalues */
   double X[n][n];            /* eigenvectors */
   /* Work arrays */
   double work[n][n];         /* work array for dsytrf */
   double rr[3][2*m][2*m];
   double W[8][m][n];
   double U[m][n];
   double V[m];

   /* Derived types */
   struct spral_ssmfe_rcid rci;            /* reverse communication data */
   struct spral_ssmfe_options options;     /* options */
   void *keep;                             /* private data */
   struct spral_ssmfe_inform inform;       /* information */

   if ( problem > 1 ) {
     sigma = 1.0;
     set_laplacian_matrix(ngrid, ngrid, n, A);
     for(int j=0; j<n; j++)
     for(int i=0; i<n; i++)
        LDLT[j][i] = (i==j) ? A[j][i] - sigma
                            : A[j][i];
     cwrap_dsytrf('L', n, &LDLT[0][0], n, ipiv, &work[0][0], n*n);
     int i = num_neg_D(n, n, LDLT, ipiv);   /* all evalues to left of sigma */
     if ( i < nep ) nep = i;
   }

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);
   /* gap between the last converged eigenvalue and the rest of the spectrum
    * must be at least 0.1 times average gap between computed eigenvalues */
   options.left_gap = -0.1;
   /* block size is small, convergence may be slow, allow more iterations */
   options.max_iterations = 300;

   /* Initialize W to lin indep vectors by randomizing */
   for(int i=0; i<n; i++)
   for(int j=0; j<m; j++)
      W[0][j][i] = spral_random_real(&state, true);

   int ncon = 0;   /* number of converged eigenpairs */
   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      switch ( problem ) {
      case 0:
        spral_ssmfe_expert_standard_double(&rci, nep, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 1:
        spral_ssmfe_expert_generalized_double(&rci, nep, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 2:
        spral_ssmfe_expert_standard_shift_double(&rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 3:
        spral_ssmfe_expert_generalized_shift_double(
           &rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 4:
        spral_ssmfe_expert_buckling_double(
           &rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      }
      switch ( rci.job ) {
      case 1:
         if ( problem < 4 )
           apply_laplacian(
              ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
              );
         else
           cblas_dcopy
              ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 2:
         if ( problem < 2 )
           apply_gauss_seidel_step (
              ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
              );
         else
            cblas_dcopy
               ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 3:
         if ( problem < 4 )
           cblas_dcopy
              ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         else
           apply_laplacian(
              ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
              );
         break;
      case 5:
         if ( rci.i < 0 ) continue;
         for(int k=0; k<rci.nx; k++) {
           int j = ncon + k;
           cblas_dcopy( n, &W[0][rci.jx+k][0], 1, &X[j][0], 1 );
         }
         ncon += rci.nx;
         break;
      case 9:
         cblas_dcopy(
            n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
            );
         cwrap_dsytrs(
            'L', n, rci.nx, &LDLT[0][0], n, ipiv, &W[rci.ky][rci.jy][0], n
            );
         break;
      case 11:
      case 12:
      case 13:
      case 14:
      case 15:
      case 16:
      case 17:
        vector_operations_d( rci, n, m, 8, ind, W, rr, V );
         break;
      case 21: // Fall through to 22
      case 22:
         if( ncon > 0 ) {
            cblas_dgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, ncon, rci.nx, n,
               1.0, &X[0][0], n, &W[rci.ky][rci.jy][0], n, 0.0, &U[0][0], n
               );
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
               -1.0, &X[0][0], n, &U[0][0], n, 1.0, &W[rci.kx][rci.jx][0], n
               );
            if ( problem == 1 || problem == 3 )
              cblas_dgemm(
                 CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
                 -1.0, &X[0][0], n, &U[0][0], n, 1.0, &W[rci.ky][rci.jy][0], n
                 );
            else if ( problem == 4 )
               apply_laplacian(
                  ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0],
                  &W[rci.ky][rci.jy][0]
                  );
         }
         break;
      case 999:
         if( rci.k == 0 ) {
            if( rci.jx > 1 ) {
               for(int j=0; j<rci.jx; j++)
               for(int i=0; i<n; i++)
                  W[0][j][i] = spral_random_real(&state, true);
            }
         }
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( ncon != ncon_x ) errors++;
   if ( problem < 2 ) {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_six[i]) > 1e-6 ) errors++;
   }
   spral_ssmfe_expert_free(&keep, &inform);

   return errors;
}

int test_expert_z(int problem) {

   const int ngrid = 20;          /* grid points along each side */
   const int n     = ngrid*ngrid; /* problem size */
   const int m     = 3;           /* dimension of the iterated subspace */

   const int ncon_x = 6;
   const double lambda_x[] = {
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01,
      2.2040061e-01
   };
   const double lambda_six[] = {
      9.5108266e-01,
      9.5108266e-01,
      8.8141871e-01,
      8.8141871e-01,
      8.4187478e-01,
      8.4187478e-01
   };

   const double complex ZERO = 0.0;
   const double complex ONE = 1.0;
   const double complex NONE = -1.0;

   int errors;
   int nep = 5;

   int state = SPRAL_RANDOM_INITIAL_SEED; /* PRNG state */

   int ind[m];                    /* permutation index */
   int ipiv[n];                   /* LDLT pivot index */
   double sigma;                  /* shift */
   double lambda[n];              /* eigenvalues */
   double complex X[n][n];        /* eigenvectors */
   /* Work arrays */
   double B[n][n];
   double LDLT[n][n];
   double work[n][n];
   double complex rr[3][2*m][2*m];
   double complex W[8][m][n];
   double complex U[m][n];
   double complex V[m];

   /* Derived types */
   struct spral_ssmfe_rciz rci;            /* reverse communication data */
   struct spral_ssmfe_options options;     /* options */
   void *keep;                             /* private data */
   struct spral_ssmfe_inform inform;       /* information */

   if ( problem > 1 ) {
     sigma = 1.0;
//     set_laplacian_matrix_z(ngrid, ngrid, n, A);
     set_laplacian_matrix(ngrid, ngrid, n, B);
     for(int j=0; j<n; j++)
     for(int i=0; i<n; i++)
        LDLT[j][i] = (i==j) ? B[j][i] - sigma
                            : B[j][i];
//        LDLT[j][i] = (i==j) ? A[j][i] - sigma
//                            : A[j][i];
     cwrap_dsytrf('L', n, &LDLT[0][0], n, ipiv, &work[0][0], n*n);
   }

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);
   /* gap between the last converged eigenvalue and the rest of the spectrum
    * must be at least 0.1 times average gap between computed eigenvalues */
   options.left_gap = -0.1;
   /* block size is small, convergence may be slow, allow more iterations */
   options.max_iterations = 200;

   /* Initialize W to lin indep vectors by randomizing */
   for(int i=0; i<n; i++)
   for(int j=0; j<m; j++)
      W[0][j][i] = spral_random_real(&state, true);

   int ncon = 0;   /* number of converged eigenpairs */
   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      switch ( problem ) {
      case 0:
        spral_ssmfe_expert_standard_double_complex(&rci, nep, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 1:
        spral_ssmfe_expert_generalized_double_complex(&rci, nep, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 2:
        spral_ssmfe_expert_standard_shift_double_complex(
           &rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 3:
        spral_ssmfe_expert_generalized_shift_double_complex(
           &rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      case 4:
        spral_ssmfe_expert_buckling_double_complex(
           &rci, sigma, nep, 0, n, lambda,
           m, &rr[0][0][0], ind, &keep, &options, &inform);
        break;
      }
      switch ( rci.job ) {
      case 1:
         if ( problem < 4 )
           apply_laplacian_z(
              ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
              );
         else
            cblas_zcopy
               ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 2:
         if ( problem < 2 )
            apply_gauss_seidel_step_z (
               ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
               );
         else
            cblas_zcopy
               ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         break;
      case 3:
         if ( problem < 4 )
           cblas_zcopy
              ( n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1 );
         else
           apply_laplacian_z(
              ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
              );
         break;
      case 5:
         if ( rci.i < 0 ) continue;
         for(int k=0; k<rci.nx; k++) {
           int j = ncon + k;
           cblas_zcopy( n, &W[0][rci.jx+k][0], 1, &X[j][0], 1 );
         }
         ncon += rci.nx;
         break;
      case 9:
         for(int i=0; i<n; i++)
           for(int j=0; j<rci.nx; j++)
              work[j][i] = W[rci.kx][rci.jx + j][i];
         cwrap_dsytrs(
            'L', n, rci.nx, &LDLT[0][0], n, ipiv, &work[0][0], n
            );
         for(int i=0; i<n; i++)
           for(int j=0; j<rci.nx; j++)
              W[rci.ky][rci.jy + j][i] = work[j][i];
         break;
      case 11:
      case 12:
      case 13:
      case 14:
      case 15:
      case 16:
      case 17:
        vector_operations_z( rci, n, m, 8, ind, W, rr, V );
         break;
      case 21: // Fall through to 22
      case 22:
         if( ncon > 0 ) {
            cblas_zgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, ncon, rci.nx, n,
               &ONE, &X[0][0], n, &W[rci.ky][rci.jy][0], n, &ZERO, &U[0][0], n
               );
            cblas_zgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
               &NONE, &X[0][0], n, &U[0][0], n, &ONE, &W[rci.kx][rci.jx][0], n
               );
            if ( problem == 1 || problem == 3 )
              cblas_zgemm(
                 CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
                 &NONE, &X[0][0], n, &U[0][0], n, &ONE, &W[rci.ky][rci.jy][0], n
                 );
            else if ( problem == 4 )
               apply_laplacian_z(
                  ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0],
                  &W[rci.ky][rci.jy][0]
                  );
         }
         break;
      case 999:
         if( rci.k == 0 ) {
            if( rci.jx > 1 ) {
               for(int j=0; j<rci.jx; j++)
               for(int i=0; i<n; i++)
                  W[0][j][i] = spral_random_real(&state, true);
            }
         }
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( ncon != ncon_x ) errors++;
   if ( problem < 2 ) {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < ncon && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_six[i]) > 1e-6 ) errors++;
   }
   spral_ssmfe_expert_free(&keep, &inform);

   return errors;
}

int test_ssmfe(void) {

  int errors = 0;
  int err;

  fprintf(stdout, "testing standard_double...");
  err = test_ssmfe_d(0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double...");
  err = test_ssmfe_d(1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_shift_double...");
  err = test_ssmfe_d(2);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_shift_double...");
  err = test_ssmfe_d(3);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing buckling_double...");
  err = test_ssmfe_d(4);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_double_complex...");
  err = test_ssmfe_z(0);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_double_complex...");
  err = test_ssmfe_z(1);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing standard_shift_double_complex...");
  err = test_ssmfe_z(2);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing generalized_shift_double_complex...");
  err = test_ssmfe_z(3);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  fprintf(stdout, "testing buckling_double_complex...");
  err = test_ssmfe_z(4);
  if ( err != 0 )
    fprintf(stdout, "%d errors\n", err);
  else
    fprintf(stdout, "ok\n");
  errors += err;

  return errors;
}

int test_ssmfe_d(int problem) {

   const int m   = 20;     /* grid points along each side */
   const int n   = m*m;    /* problem size */

   const int ncon_x = 6;
   const double sigma = 1.0;  /* shift */
   const double lambda_x[] = {
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01,
      2.2040061e-01
   };
   const double lambda_six[] = {
      8.4187478e-01,
      8.4187478e-01,
      8.8141871e-01,
      8.8141871e-01,
      9.5108266e-01,
      9.5108266e-01
   };

   int errors;
   int nep = 5;

   int ipiv[n];                           /* LDLT pivot index */

   double lambda[n];                  /* eigenvalues */
   double X[n][n];                    /* eigenvectors */
   double A[n][n];                        /* matrix */
   double LDLT[n][n];                     /* factors */
   double work[n][n];                     /* work array for dsytrf */

   struct spral_ssmfe_rcid rci;           /* reverse communication data */
   struct spral_ssmfe_options options;    /* options */
   void *keep;                            /* private data */
   struct spral_ssmfe_inform inform;      /* information */

   if ( problem > 1 ) {
     set_laplacian_matrix( m, m, n, A );
     for(int j=0; j<n; j++)
     for(int i=0; i<n; i++)
        LDLT[j][i] = (i==j) ? A[j][i] - sigma
                            : A[j][i];
     cwrap_dsytrf( 'L', n, &LDLT[0][0], n, ipiv, &work[0][0], n*n );
   }

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);
   /* gap between the last converged eigenvalue and the rest of the spectrum
    * must be at least 0.1 times average gap between computed eigenvalues */
   options.left_gap = -0.1;

   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      switch ( problem ) {
      case 0:
         spral_ssmfe_standard_double( &rci, nep, 2*nep, lambda, n, &X[0][0], n,
            &keep, &options, &inform );
         break;
      case 1:
         spral_ssmfe_standard_double( &rci, nep, 2*nep, lambda, n, &X[0][0], n,
            &keep, &options, &inform );
         break;
      case 2:
         spral_ssmfe_standard_shift_double( &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      case 3:
         spral_ssmfe_generalized_shift_double( &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      case 4:
         spral_ssmfe_buckling_double( &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      }
      switch ( rci.job ) {
      case 1:
         if ( problem < 4 )
            apply_laplacian(m, m, rci.nx, rci.x, rci.y);
         else
            cblas_dcopy( n*rci.nx, rci.x, 1, rci.y, 1 );
         break;
      case 2:
         if ( problem < 2 )
            apply_gauss_seidel_step(m, m, rci.nx, rci.x, rci.y);
         break;
      case 3:
         if ( problem < 4 )
            cblas_dcopy( n*rci.nx, rci.x, 1, rci.y, 1 );
         else
            apply_laplacian(m, m, rci.nx, rci.x, rci.y);
         break;
      case 9:
         cblas_dcopy( n*rci.nx, rci.x, 1, rci.y, 1 );
         cwrap_dsytrs( 'L', n, rci.nx, &LDLT[0][0], n, ipiv, rci.y, n );
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( inform.left != ncon_x ) errors++;
   if ( problem < 2 ) {
     for ( int i = 0; i < inform.left && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < inform.left && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_six[i]) > 1e-6 ) errors++;
   }
   spral_ssmfe_free_double(&keep, &inform);
   return errors;
}

int test_ssmfe_z(int problem) {

   const int m   = 20;     /* grid points along each side */
   const int n   = m*m;    /* problem size */

   const int ncon_x = 6;
   const double sigma = 1.0;  /* shift */
   const double lambda_x[] = {
      4.4676695e-02,
      1.1119274e-01,
      1.1119274e-01,
      1.7770878e-01,
      2.2040061e-01,
      2.2040061e-01
   };
   const double lambda_six[] = {
      8.4187478e-01,
      8.4187478e-01,
      8.8141871e-01,
      8.8141871e-01,
      9.5108266e-01,
      9.5108266e-01
   };

   int errors;
   int nep = 5;

   int ipiv[n];                           /* LDLT pivot index */

   double lambda[n];                  /* eigenvalues */
   double complex X[n][n];                    /* eigenvectors */
   double A[n][n];                        /* matrix */
   double LDLT[n][n];                     /* factors */
   double work[n][n];                     /* work array for dsytrf */

   struct spral_ssmfe_rciz rci;           /* reverse communication data */
   struct spral_ssmfe_options options;    /* options */
   void *keep;                            /* private data */
   struct spral_ssmfe_inform inform;      /* information */

   if ( problem > 1 ) {
     set_laplacian_matrix( m, m, n, A );
     for(int j=0; j<n; j++)
     for(int i=0; i<n; i++)
        LDLT[j][i] = (i==j) ? A[j][i] - sigma
                            : A[j][i];
     cwrap_dsytrf( 'L', n, &LDLT[0][0], n, ipiv, &work[0][0], n*n );
   }

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);
   /* gap between the last converged eigenvalue and the rest of the spectrum
    * must be at least 0.1 times average gap between computed eigenvalues */
   options.left_gap = -0.1;

   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      switch ( problem ) {
      case 0:
         spral_ssmfe_standard_double_complex(
            &rci, nep, 2*nep, lambda, n, &X[0][0], n,
            &keep, &options, &inform );
         break;
      case 1:
         spral_ssmfe_standard_double_complex(
            &rci, nep, 2*nep, lambda, n, &X[0][0], n,
            &keep, &options, &inform );
         break;
      case 2:
         spral_ssmfe_standard_shift_double_complex(
            &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      case 3:
         spral_ssmfe_generalized_shift_double_complex(
            &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      case 4:
         spral_ssmfe_buckling_double_complex(
            &rci, sigma, nep, 0, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform );
         break;
      }
      switch ( rci.job ) {
      case 1:
         if ( problem < 4 )
            apply_laplacian_z( m, m, rci.nx, rci.x, rci.y );
         else
            cblas_zcopy( n*rci.nx, rci.x, 1, rci.y, 1 );
         break;
      case 2:
         if ( problem < 2 )
            apply_gauss_seidel_step_z( m, m, rci.nx, rci.x, rci.y );
         break;
      case 3:
         if ( problem < 4 )
            cblas_zcopy( n*rci.nx, rci.x, 1, rci.y, 1 );
         else
            apply_laplacian_z( m, m, rci.nx, rci.x, rci.y );
         break;
      case 9:
         zdcopy( n*rci.nx, rci.x, &work[0][0] );
         cwrap_dsytrs( 'L', n, rci.nx, &LDLT[0][0], n, ipiv, &work[0][0], n );
         dzcopy( n*rci.nx, &work[0][0], rci.y );
         break;
      default:
         goto finished;
      }
   }
finished:
   errors = 0;
   if ( inform.flag != 0 ) errors++;
   if ( inform.left != ncon_x ) errors++;
   if ( problem < 2 ) {
     for ( int i = 0; i < inform.left && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_x[i]) > 1e-6 ) errors++;
   }
   else {
     for ( int i = 0; i < inform.left && i < ncon_x; i++ )
        if ( fabs(lambda[i] - lambda_six[i]) > 1e-6 ) errors++;
   }
   if ( problem == -2 ) {
   printf("%d eigenpairs converged in %d iterations\n", inform.left, inform.iteration);
   for(int i=0; i<inform.left; i++)
      printf(" lambda[%1d] = %13.7e\n", i, lambda[i]);
   }
   spral_ssmfe_free_double_complex(&keep, &inform);
   return errors;
}

void zdcopy( int n, double complex* x, double* y ) {
   for ( int i = 0; i < n; i++ )
      y[i] = x[i];
}

void dzcopy( int n, double* x, double complex* y ) {
   for ( int i = 0; i < n; i++ )
      y[i] = x[i];
}

void vector_operations_d
  ( struct spral_ssmfe_rcid rci,
    int n, int m, int kw, int ind[m],
    double W[kw][m][n],
    double rr[3][2*m][2*m],
    double U[m]) {

   switch ( rci.job ) {
      case 11:
         if ( rci.i == 0 ) {
            if ( rci.kx != rci.ky || rci.jx > rci.jy ) {
               cblas_dcopy(
                  n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
                  );
            } else if ( rci.jx < rci.jy ) {
               for(int j=rci.nx-1; j>=0; j--)
                  cblas_dcopy(
                     n, &W[rci.kx][rci.jx+j][0], 1, &W[rci.ky][rci.jy+j][0], 1
                     );
            }
         } else {
            for(int i=0; i<n; i++) {
               for(int j=0; j<rci.nx; j++)
                  U[j] = W[rci.kx][ind[j]][i];
               for(int j=0; j<rci.nx; j++)
                  W[rci.kx][j][i] = U[j];
               if(rci.ky != rci.kx) {
                  for(int j=0; j<rci.nx; j++)
                    U[j] = W[rci.ky][ind[j]][i];
                  for(int j=0; j<rci.nx; j++)
                    W[rci.ky][j][i] = U[j];
               }
            }
         }
         break;
      case 12:
         for(int i=0; i<rci.nx; i++)
           rr[rci.k][rci.j+i][rci.i+i] =
             cblas_ddot(
                n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1
                );
         break;
      case 13:
         for(int i=0; i<rci.nx; i++) {
            if( rci.kx == rci.ky ) {
               double s = cblas_dnrm2(n, &W[rci.kx][rci.jx+i][0], 1);
               if( s > 0 )
                  cblas_dscal(n, 1/s, &W[rci.kx][rci.jx+i][0], 1);
            } else {
               double s = sqrt(fabs(cblas_ddot(
                  n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1)
                  ));
               if ( s > 0 ) {
                  cblas_dscal(n, 1/s, &W[rci.kx][rci.jx+i][0], 1);
                  cblas_dscal(n, 1/s, &W[rci.ky][rci.jy+i][0], 1);
               } else {
                  for(int j=0; j<n; j++)
                     W[rci.ky][rci.jy+i][j] = 0.0;
               }
            }
         }
         break;
      case 14:
         for(int i=0; i<rci.nx; i++) {
           double s = -rr[rci.k][rci.j+i][rci.i+i];
           cblas_daxpy(
              n, s, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1
              );
         }
         break;
      case 15:
         if ( rci.nx > 0 && rci.ny > 0 )
            cblas_dgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, rci.nx, rci.ny, n,
               rci.alpha, &W[rci.kx][rci.jx][0], n, &W[rci.ky][rci.jy][0], n,
               rci.beta, &rr[rci.k][rci.j][rci.i], 2*m
               );
         break;
      case 16: // Fall through to 17
      case 17:
         if( rci.job == 17 ) {
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
               1.0, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i], 2*m,
               0.0, &W[rci.ky][rci.jy][0], n
               );
            cblas_dcopy(
               n*rci.ny, &W[rci.ky][rci.jy][0], 1, &W[rci.kx][rci.jx][0], 1
               );
         } else {
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
               rci.alpha, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i],
               2*m, rci.beta, &W[rci.ky][rci.jy][0], n
               );
         }
         break;
   }
}

void vector_operations_z
  ( struct spral_ssmfe_rciz rci,
    int n, int m, int kw, int ind[m],
    double complex W[kw][m][n],
    double complex rr[3][2*m][2*m],
    double complex U[m]) {

   const double complex ZERO = 0.0;
   const double complex ONE = 1.0;

   double complex z;

    switch ( rci.job ) {
    case 11:
       if ( rci.i == 0 ) {
          if ( rci.kx != rci.ky || rci.jx > rci.jy ) {
             cblas_zcopy(
                n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1
                );
          } else if ( rci.jx < rci.jy ) {
             for(int j=rci.nx-1; j>=0; j--)
                cblas_zcopy(
                   n, &W[rci.kx][rci.jx+j][0], 1, &W[rci.ky][rci.jy+j][0], 1
                   );
          }
       } else {
          for(int i=0; i<n; i++) {
             for(int j=0; j<rci.nx; j++)
                U[j] = W[rci.kx][ind[j]][i];
             for(int j=0; j<rci.nx; j++)
                W[rci.kx][j][i] = U[j];
             if(rci.ky != rci.kx) {
                for(int j=0; j<rci.nx; j++)
                  U[j] = W[rci.ky][ind[j]][i];
                for(int j=0; j<rci.nx; j++)
                  W[rci.ky][j][i] = U[j];
             }
          }
       }
       break;
    case 12:
       for(int i=0; i<rci.nx; i++) {
          cblas_zdotc_sub(
              n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1, &z
              );
          rr[rci.k][rci.j+i][rci.i+i] = z;
       }
       break;
    case 13:
       for(int i=0; i<rci.nx; i++) {
          if( rci.kx == rci.ky ) {
             double s = cblas_dznrm2(n, &W[rci.kx][rci.jx+i][0], 1);
             if( s > 0 ) {
                z = ONE/s;
                cblas_zscal(n, &z, &W[rci.kx][rci.jx+i][0], 1);
             }
          } else {
             cblas_zdotc_sub(
                n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1, &z
                );
             double s = sqrt(fabs(z));
             if ( s > 0 ) {
                z = ONE/s;
                cblas_zscal(n, &z, &W[rci.kx][rci.jx+i][0], 1);
                cblas_zscal(n, &z, &W[rci.ky][rci.jy+i][0], 1);
             } else {
                for(int j=0; j<n; j++)
                   W[rci.ky][rci.jy+i][j] = 0.0;
             }
          }
       }
       break;
    case 14:
       for(int i=0; i<rci.nx; i++) {
         double s = -rr[rci.k][rci.j+i][rci.i+i];
         z = s;
         cblas_zaxpy(
            n, &z, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1
            );
       }
       break;
    case 15:
       if ( rci.nx > 0 && rci.ny > 0 )
          cblas_zgemm(
             CblasColMajor, CblasTrans, CblasNoTrans, rci.nx, rci.ny, n,
             &rci.alpha, &W[rci.kx][rci.jx][0], n, &W[rci.ky][rci.jy][0], n,
             &rci.beta, &rr[rci.k][rci.j][rci.i], 2*m
             );
       break;
    case 16: // Fall through to 17
    case 17:
       if( rci.job == 17 ) {
          cblas_zgemm(
             CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
             &ONE, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i], 2*m,
             &ZERO, &W[rci.ky][rci.jy][0], n
             );
          cblas_zcopy(
             n*rci.ny, &W[rci.ky][rci.jy][0], 1, &W[rci.kx][rci.jx][0], 1
             );
       } else {
          cblas_zgemm(
             CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
             &rci.alpha, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i],
             2*m, &rci.beta, &W[rci.ky][rci.jy][0], n
             );
       }
       break;
     }
}

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_interp.h>

#include "allvars.h"
#include "proto.h"

#define SNR (120.0)
int n_con_sim=301, n_line_sim=301;
double sim_error=10.0*1.0/SNR, sim_rel_error=1.0/SNR;

double *PQmat, *PSmat, *Peigens, *Peigens_mat, *Peigens_vecs;

gsl_rng_type const * gsl_T_sim;
gsl_rng * gsl_r_sim;
gsl_interp_accel *gsl_acc, *gsl_acc_error;
gsl_interp  *gsl_linear, *gsl_linear_error;

void simulate_con(double sigma, double tau, double alpha, double ave_con)
{
  FILE *fscon_out;

  double *Py, *Prandvec;
  int i, info;

  Py = array_malloc(n_con_sim);
  Prandvec = array_malloc(n_con_sim);

  simulate_con_init();

  set_covar_Simmat(sigma, tau, alpha);
  memcpy(PQmat, PSmat, n_con_sim*n_con_sim*sizeof(double));

/*  eigen_sym_mat(PQmat, n_con_sim, Peigens, &info);
  memcpy(Peigens_vecs, PQmat, n_con_sim*n_con_sim*sizeof(double));

  for(i=0; i<n_con_sim; i++)
  {
    for(j=0; j<n_con_sim; j++)
      Peigens_mat[i*n_con_sim+j] = 0.0;
    Peigens_mat[i*n_con_sim+i] = sqrt(Peigens[i]);
  }*/

  for(i=0;i<n_con_sim; i++)
  {
    Prandvec[i]  = gsl_ran_gaussian(gsl_r_sim, 1.0);
  }

  Chol_decomp_U(PQmat, n_con_sim, &info);
  multiply_matvec(PQmat, Prandvec, n_con_sim, Py);

//   multiply_matvec(Peigens_mat, Prandvec, n_con_sim, Pybuf);
//   multiply_matvec(Peigens_vecs, Pybuf, n_con_sim, Py);

  gsl_rng_set(gsl_r_sim, time(NULL));
  for(i=0;i<n_con_sim;i++)
  {
    Fcon[i] = Py[i] + ave_con + sim_error * gsl_ran_gaussian(gsl_r_sim, 1.0);
    Fcerrs[i] = sim_error;
  }

  fscon_out = fopen("data/sim_con.txt", "w");

  for(i=0;i<n_con_sim;i++)
  {
    if(Tcon[i]>150.0)
      fprintf(fscon_out, "%f\t%f\t%f\n", Tcon[i], Fcon[i], Fcerrs[i]);
  }

  gsl_interp_init(gsl_linear, Tcon, Fcon, n_con_sim);
  gsl_interp_init(gsl_linear_error, Tcon, Fcerrs, n_con_sim);

  fclose(fscon_out);
}

void simulate_con_init()
{
  int i;
  double Tcon_min, Tcon_max, Tline_min, Tline_max, dTs;

  PSmat = array_malloc(n_con_sim*n_con_sim);
  PQmat = array_malloc(n_con_sim*n_con_sim);
  Peigens = array_malloc(n_con_sim);
  Peigens_vecs = array_malloc(n_con_sim*n_con_sim);
  Peigens_mat = array_malloc(n_con_sim*n_con_sim);

  Tcon = array_malloc(n_con_sim);
  Fcon = array_malloc(n_con_sim);
  Fcerrs = array_malloc(n_con_sim);

  Tcon_min = 0.0;
  Tcon_max = 300.0;
  dTs = (Tcon_max - Tcon_min) / (n_con_sim-1.0);
  for(i=0; i<n_con_sim; i++)
  {
    Tcon[i] = dTs * i + Tcon_min;
  }

  Tline = array_malloc(n_line_sim);
  Fline = array_malloc(n_line_sim);
  Flerrs = array_malloc(n_line_sim);

  Tline_min = 0.0;
  Tline_max = 300.0;
  dTs = (Tline_max - Tline_min) / (n_line_sim-1.0);
  for(i=0; i<n_line_sim; i++)
  {
    Tline[i] = dTs * i + Tline_min;
  }

  gsl_acc = gsl_interp_accel_alloc();
  gsl_acc_error = gsl_interp_accel_alloc();

  gsl_linear = gsl_interp_alloc(gsl_interp_linear, n_con_sim);
  gsl_linear_error = gsl_interp_alloc(gsl_interp_linear, n_con_sim);

  gsl_T_sim = gsl_rng_default;
  gsl_r_sim = gsl_rng_alloc (gsl_T_sim);
  //gsl_rng_set(gsl_r_sim, 12);
  gsl_rng_set(gsl_r_sim, time(NULL));

  return;
}

void set_covar_Simmat(double sigma, double tau, double alpha)
{
  double t1, t2, nerr;
  int i, j;

  for(i=0; i<n_con_sim; i++)
  {
    t1 = Tcon[i];
    for(j=0; j<=i; j++)
    {
      t2 = Tcon[j];
      PSmat[i*n_con_sim+j] = sigma*sigma* exp (- pow (fabs(t1-t2) / tau, alpha));
      PSmat[j*n_con_sim+i] = PSmat[i*n_con_sim+j];
    }
    nerr = 0.0; //sim_error;
    PSmat[i*n_con_sim+i] += nerr * nerr;
  }
  return;
}

void simulate_line()
{
  int i, j;
  double flux, err, tline, tcon, taup, norm, fcon, fcon_err;
  double dtau, sig;
  int type=1;

  dtau = 30.0/(ntau-1.0);
  for(i=0; i<ntau; i++)
  {
    TF_tau[i] =  dtau * i;
    TF[i] = 0.0;
  }
// top-hat transfer function
  FILE *ftf;
  ftf=fopen("data/transfer_sim.txt", "w");
  norm = 0.0;
  for(i=0; i<ntau; i++)
  {
    if(type==0)  // two top-hats
    {
      if(TF_tau[i]<=5.0)
        TF[i] = 0.0;
      else if(TF_tau[i]<10.0)
        TF[i] = 1.0;
      else if(TF_tau[i]<15.0)
        TF[i] = 0.0;
      else if(TF_tau[i]<20.0)
        TF[i] = 2.0;
      else
        TF[i] = 0.0;

      norm += TF[i] * dtau;
    }
    if(type==1) // two Gaussians
    {
      sig = 2.5;
      TF[i] = exp(-0.5*pow(TF_tau[i]-10.0, 2.0)/sig/sig) + exp(-0.5*pow(TF_tau[i]-18.0, 2.0)/sig/sig);
      norm += TF[i] * dtau;
    }
  }
  for(i=0; i<ntau; i++)
  {
    TF[i] /= norm;
    fprintf(ftf, "%f %f\n", TF_tau[i], TF[i]);
  }
  fclose(ftf);

  FILE *fp;
  fp = fopen("data/sim_line.txt", "w");

  for(i=0; i<n_line_sim; i++)
  {
    flux = 0.0;
    err = 0.0;
    tline = Tline[i]; //4540.0 + 70.0/(nline_data-1.0) * i;
    for(j=0; j<ntau; j++)
    {
      taup = TF_tau[j];
      tcon = tline - taup;
      if(tcon >=Tcon[0] && tcon <= Tcon[n_con_sim-1])
      {
        fcon = gsl_interp_eval(gsl_linear, Tcon, Fcon, tcon, gsl_acc);
        fcon_err = gsl_interp_eval(gsl_linear_error, Tcon, Fcerrs, tcon, gsl_acc_error);
        flux += fcon * TF[j];
        err += pow(fcon_err *TF[j] * dtau, 2);
      }
    }
    flux *= dtau;
    err = sqrt(err);
    Fline[i] = flux;
    Flerrs[i] = sim_error;
  }

  gsl_rng_set(gsl_r_sim, time(NULL)); // reset the seed of the random generator
  for(i=0; i<n_line_sim; i++)
  {
    if(Tline[i]>170.0)
      fprintf(fp, "%f %f %f\n", Tline[i], Fline[i] + Flerrs[i]*gsl_ran_gaussian(gsl_r_sim, 1.0), Flerrs[i]);
  }
  fclose(fp);
}

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


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 rendersysd3d11_8555fd57_19e9_4747_97b5_8d3bfb0fb50f_h
#define rendersysd3d11_8555fd57_19e9_4747_97b5_8d3bfb0fb50f_h

#include <gslib/type.h>
#include <ariel/rendersys.h>

__ariel_begin__

class rendersys_d3d11:
    public rendersys
{
public:
    rendersys_d3d11();
    virtual ~rendersys_d3d11();
    virtual bool setup(uint hwnd, const configs& cfg) override;
    virtual void destroy() override;
    virtual void setup_pipeline_state() override;
    virtual render_blob* compile_shader_from_file(const gchar* file, const gchar* entry, const gchar* sm, render_include* inc) override;
    virtual render_blob* compile_shader_from_memory(const char* src, int len, const gchar* name, const gchar* entry, const gchar* sm, render_include* inc) override;
    virtual vertex_shader* create_vertex_shader(const void* ptr, size_t len) override;
    virtual pixel_shader* create_pixel_shader(const void* ptr, size_t len) override;
    virtual compute_shader* create_compute_shader(const void* ptr, size_t len) override;
    virtual geometry_shader* create_geometry_shader(const void* ptr, size_t len) override;
    virtual hull_shader* create_hull_shader(const void* ptr, size_t len) override;
    virtual domain_shader* create_domain_shader(const void* ptr, size_t len) override;
    virtual vertex_format* create_vertex_format(const void* ptr, size_t len, vertex_format_desc vfdesc[], uint n) override;
    virtual vertex_buffer* create_vertex_buffer(uint stride, uint count, bool read, bool write, uint usage, const void* ptr) override;
    virtual index_buffer* create_index_buffer(uint count, bool read, bool write, uint usage, const void* ptr) override;
    virtual constant_buffer* create_constant_buffer(uint stride, bool read, bool write, const void* ptr) override;
    virtual shader_resource_view* create_shader_resource_view(render_resource* res) override;
    virtual depth_stencil_view* create_depth_stencil_view(render_resource* res) override;
    virtual unordered_access_view* create_unordered_access_view(render_resource* res) override;
    virtual sampler_state* create_sampler_state(sampler_state_filter filter) override;
    virtual texture2d* create_texture2d(const image& img, uint mips, uint usage, uint bindflags, uint cpuflags, uint miscflags) override;
    virtual texture2d* create_texture2d(int width, int height, uint format, uint mips, uint usage, uint bindflags, uint cpuflags, uint miscflags) override;
    virtual void load_with_mips(texture2d* tex, const image& img) override;
    virtual void update_buffer(void* buf, int size, const void* ptr) override;
    virtual void set_vertex_format(vertex_format* vfmt) override;
    virtual void set_vertex_buffer(vertex_buffer* vb, uint stride, uint offset) override;
    virtual void set_index_buffer(index_buffer* ib, uint offset) override;
    virtual void begin_render() override;
    virtual void end_render() override;
    virtual void set_render_option(render_option opt, uint val) override;
    virtual void set_vertex_shader(vertex_shader* vs) override;
    virtual void set_pixel_shader(pixel_shader* ps) override;
    virtual void set_geometry_shader(geometry_shader* gs) override;
    virtual void set_viewport(const viewport& vp) override;
    virtual void set_constant_buffer(uint slot, constant_buffer* cb, shader_type st) override;
    virtual void set_sampler_state(uint slot, sampler_state* sstate, shader_type st) override;
    virtual void set_shader_resource(uint slot, shader_resource_view* srv, shader_type st) override;
    virtual void draw(uint count, uint start) override;
    virtual void draw_indexed(uint count, uint start, int base) override;
    virtual void capture_screen(image& img, const rectf& rc, int buff_id) override;
    virtual void enable_alpha_blend(bool b) override;
    virtual void enable_depth(bool b) override;

protected:
    D3D_DRIVER_TYPE         _drvtype        = D3D_DRIVER_TYPE_NULL;
    D3D_FEATURE_LEVEL       _level          = D3D_FEATURE_LEVEL_11_0;
    render_device*          _device         = nullptr;
    render_context*         _context        = nullptr;
    render_swap_chain*      _swapchain      = nullptr;
    render_target_view*     _rtview         = nullptr;
    render_blend_state*     _blendstate     = nullptr;
    depth_stencil_view*     _dsview         = nullptr;
    render_raster_state*    _rasterstate    = nullptr;
    render_depth_state*     _depthstate     = nullptr;
    uint                    _msaa_x         = 0;
    bool                    _vsync          = false;
    bool                    _fullscreen     = false;
    bool                    _msaa           = false;

protected:
    void install_configs(const configs& cfg);

public:
    render_device* get_device() const { return _device; }
    render_context* get_immediate_context() const { return _context; }
};

template<class res_class>
inline render_resource* convert_to_resource(res_class* p)
{ return static_cast<render_resource*>(p); }

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib