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_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 <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/parameters.c"
#include "../cosmolike_core/emu17/P_cb/emu.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift_spline.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/pt.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"
#include "../cosmolike_core/theory/IA.c"
#include "../cosmolike_core/theory/cluster.c"
#include "../cosmolike_core/theory/BAO.c"
#include "../cosmolike_core/theory/external_prior.c"
#include "../cosmolike_core/theory/covariances_3D.c"
#include "../cosmolike_core/theory/covariances_fourier.c"
#include "../cosmolike_core/theory/covariances_cluster.c"
#include "init_emu.c"

void run_cov_N_N (char *OUTFILE, char *PATH, int nzc1, int nzc2,int start);
void run_cov_cgl_N (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int nzc2, int start);
void run_cov_cgl_cgl (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int start);
void run_cov_cgl_cgl_all (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster);
void run_cov_shear_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_shear_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);
void run_cov_ggl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_ggl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);
void run_cov_cl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_cl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);

void run_cov_ggl_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_shear_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);

void run_cov_N_N (char *OUTFILE, char *PATH, int nzc1, int nzc2,int start)
{
  int nN1, nN2,i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      i += Cluster.N200_Nbin*nzc1+nN1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;

      cov =cov_N_N(nzc1,nN1, nzc2, nN2);
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,0.0,0.0, nzc1, nN1, nzc2, nN2,cov,0.0);
    }
  }
  fclose(F1);
}

void run_cov_cgl_N (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int nzc2, int start)
{
  int nN1, nN2, nl1, nzc1, nzs1,i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc1 = ZC(N1);
  nzs1 = ZSC(N1);
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
     for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
       i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
       i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
       j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
       j += Cluster.N200_Nbin*nzc2+nN2;

       cov =cov_cgl_N(ell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2);
       fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell_Cluster[nl1], 0., nzc1, nzs1, nzc2, nN2,cov,0.);
     }
   }
 }
 fclose(F1);
}

void run_cov_cgl_cgl (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int start)
{
  int nN1, nN2, nl1, nzc1, nzs1, nl2, nzc2, nzs2,i,j;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc1 = ZC(N1);
  nzs1 = ZSC(N1);
  nzc2 = ZC(N2);
  nzs2 = ZSC(N2);
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
      for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
        for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
          i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
          i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
          j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
          j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

          c_g = 0;
          c_ng = cov_NG_cgl_cgl(ell_Cluster[nl1],ell_Cluster[nl2],nzc1,nN1, nzs1, nzc2, nN2,nzs2);
          if (nl2 == nl1){c_g =cov_G_cgl_cgl(ell_Cluster[nl1],dell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2,nzs2);}
          fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell_Cluster[nl1],ell_Cluster[nl2], nzc1, nzs1, nzc2, nzs2,c_g, c_ng);
        }
      }
    }
  }
  fclose(F1);
}

void run_cov_cgl_cgl_all (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster)
{
  int nN1, nN2, nl1, nzc1, nzs1, nl2, nzc2, nzs2,i,j, N1,N2;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s",PATH,OUTFILE);
  F1 =fopen(filename,"w");
  for (N1 = 0; N1 < tomo.cgl_Npowerspectra; N1 ++){
    for (N2 = 0; N2 < tomo.cgl_Npowerspectra; N2 ++){
      nzc1 = ZC(N1);
      nzs1 = ZSC(N1);
      nzc2 = ZC(N2);
      nzs2 = ZSC(N2);
      for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
        for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
          for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
            for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
              i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
              i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
              j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
              j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;
              c_g = 0;
              c_ng = cov_NG_cgl_cgl(ell_Cluster[nl1],ell_Cluster[nl2],nzc1,nN1, nzs1, nzc2, nN2,nzs2);
              if (nl2 == nl1){c_g =cov_G_cgl_cgl(ell_Cluster[nl1],dell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2,nzs2);}
              fprintf(F1,"%d %d %e %e %d %d %d  %d %d %d  %e %e\n",i,j,ell_Cluster[nl1],ell_Cluster[nl2], nzc1, nN1,nzs1, nzc2, nN2,nzs2,c_g, c_ng);
            }
          }
        }
      }
    }
  }
  fclose(F1);
}

void run_cov_shear_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start)
{
  int nz1,nz2, nN2, nl1, nzc1, i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nz1 = Z1(N1);
  nz2 = Z2(N1);
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      cov = 0.;
      i = like.Ncl*N1+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;

      if (ell[nl1] < like.lmax_shear){cov =cov_shear_N(ell[nl1],nz1,nz2, nzc2, nN2);}
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., nz1, nz2, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_shear_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN1, nN2, nzs1, nzs2,nl2, nzc2, nzs3,i,j;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzs1 = Z1(N1);
  nzs2 = Z2(N1);
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*N1+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;
        c_g = 0.;
        c_ng = 0.;
        if (ell[nl1] < like.lmax_shear){
          c_ng = cov_NG_shear_cgl(ell[nl1],ell_Cluster[nl2],nzs1, nzs2, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.001){
            c_g =cov_G_shear_cgl(ell[nl1],dell_Cluster[nl2],nzs1,nzs2, nzc2, nN2,nzs3);
          }
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], nzs1, nzs2, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_ggl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start)
{
  int zl,zs, nN2, nl1, nzc1, i,j;
  double cov,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(N1);
  zs = ZS(N1);
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+N1)+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;
      cov = 0.;
      weight = test_kmax(ell[nl1],zl);
      if (weight){
        cov =cov_ggl_N(ell[nl1],zl,zs, nzc2, nN2);
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., zl, zs, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_ggl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN2, zl, zs, nzs1, nl2, nzc2, nzs3,i,j;
  double c_g, c_ng,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(N1);
  zs = ZS(N1);
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*(tomo.shear_Npowerspectra+N1)+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

        c_g = 0; c_ng = 0.;
        weight = test_kmax(ell[nl1],zl);
        if (weight){
          c_ng = cov_NG_ggl_cgl(ell[nl1],ell_Cluster[nl2],zl,zs, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.1){c_g =cov_G_ggl_cgl(ell[nl1],dell_Cluster[nl2],zl,zs, nzc2, nN2,nzs3);}
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], zl, zs, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_cl_N (char *OUTFILE, char *PATH, double *ell, double *dell,int N1, int nzc2, int start)
{
  int zl,zs, nN2, nl1, nzc1, i,j;
  double cov,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+N1)+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;
      cov = 0.;
      weight = test_kmax(ell[nl1],N1);
      if (weight){
        cov =cov_cl_N(ell[nl1],N1,N1,nzc2,nN2);
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., N1, N1, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_cl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN2,nzc2, nzs3,i,j,nl2;
  double c_g, c_ng,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+N1)+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

        c_g = 0; c_ng = 0.;
        weight = test_kmax(ell[nl1],N1);
        if (weight){
          c_ng = cov_NG_cl_cgl(ell[nl1],ell_Cluster[nl2],N1,N1, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.1){
            c_g =cov_G_cl_cgl(ell[nl1],dell_Cluster[nl2],N1,N1, nzc2, nN2,nzs3);
          }
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], N1,N1, nzc2, nzs3,c_g, c_ng);
        //printf("%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], N1,N1, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_ggl_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl,zs,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(n1); zs = ZS(n1);
  printf("\nN_ggl = %d (%d, %d)\n", n1,zl,zs);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(zl,z3)*test_zoverlap(zl,z4)){
          c_ng = cov_NG_gl_shear_tomo(ell[nl1],ell[nl2],zl,zs,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_gl_shear_tomo(ell[nl1],dell[nl1],zl,zs,z3,z4);
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],zl,zs,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl = %d \n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(z1,z3)*test_zoverlap(z1,z4)){
          c_ng = cov_NG_cl_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_cl_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        }
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,zl,zs,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  zl = ZL(n2); zs = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl,zs);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (z1 == zl){
        weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],zl);
        if (weight){
          c_ng = cov_NG_cl_gl_tomo(ell[nl1],ell[nl2],z1,z2,zl,zs);
          if (nl1 == nl2){
            c_g =  cov_G_cl_gl_tomo(ell[nl1],dell[nl1],z1,z2,zl,zs);
          }
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2, ell[nl1],ell[nl2],z1,z2,zl,zs,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  z3 = n2; z4 = n2;
  printf("N_cl_2 = %d\n", n2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (z1 == z3){
        weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],z3);
        if (weight) {
          c_ng = cov_NG_cl_cl_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_cl_cl_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra + n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl1,zl2,zs1,zs2,nl1,nl2, weight;
  double c_ng, c_g;
  double fsky = survey.area/41253.0;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");

  zl1 = ZL(n1); zs1 = ZS(n1);
  printf("\nN_tomo_1 = %d (%d, %d)\n", n1,zl1,zs1);
  zl2 = ZL(n2); zs2 = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl2,zs2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl1)*test_kmax(ell[nl2],zl2);
      if (weight && zl1 == zl2) {
        c_ng = cov_NG_gl_gl_tomo(ell[nl1],ell[nl2],zl1,zs1,zl2,zs2);
      }
      if (nl1 == nl2){
        c_g =  cov_G_gl_gl_tomo(ell[nl1],dell[nl1],zl1,zs1,zl2,zs2);
      }
      if (weight ==0 && n2 != n1){
        c_g = 0;
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2, ell[nl1],ell[nl2],zl1,zs1,zl2,zs2,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_shear_shear(char *OUTFILE, char *PATH, double *ell, double *dell,int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  z1 = Z1(n1); z2 = Z2(n1);
  printf("N_shear = %d\n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n",n2,z3,z4);
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (ell[nl1] < like.lmax_shear && ell[nl2] < like.lmax_shear){
        c_ng = cov_NG_shear_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
      }
      if (nl1 == nl2){
        c_g =  cov_G_shear_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        if (ell[nl1] > like.lmax_shear && n1!=n2){c_g = 0.;}
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",like.Ncl*n1+nl1,like.Ncl*(n2)+nl2,ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
      //printf("%d %d %e %e %d %d %d %d %e %e\n", like.Ncl*n1+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

int main(int argc, char** argv)
{

  int i,l,m,n,o,s,p,nl1,t,k;
  char OUTFILE[400],filename[400],arg1[400],arg2[400];

  int N_scenarios=36;
  double area_table[36]={7623.22,14786.3,9931.47,8585.43,17681.8,15126.9,9747.99,8335.08,9533.42,18331.3,12867.8,17418.9,19783.1,12538.8,15260.0,16540.7,19636.8,11112.7,10385.5,16140.2,18920.1,17976.2,11352.0,9214.77,16910.7,11995.6,16199.8,14395.1,8133.86,13510.5,19122.3,15684.5,12014.8,14059.7,10919.3,13212.7};
  double nsource_table[36]={13.991,33.3975,17.069,28.4875,35.3643,10.3802,11.1105,29.5904,31.4608,15.7463,13.3066,9.07218,9.58719,16.3234,25.8327,10.7415,38.0412,32.8821,19.8893,27.0369,15.1711,14.2418,19.1851,26.926,22.0012,12.6553,18.6304,11.9787,36.8728,22.5265,17.4381,12.3424,10.0095,23.5979,20.8771,24.8956};

  double nlens_table[36]={22.9726 ,61.5948 ,28.7809 ,51.4354 ,65.7226 ,16.3777 ,17.6899 ,53.6984 ,57.562 ,26.266 ,21.703 ,14.0587 ,14.9667 ,27.36 ,46.0364 ,17.0253 ,71.39 ,60.5184 ,34.2283 ,48.4767 ,25.1812 ,23.44 ,32.8578 ,48.2513 ,38.3766 ,20.5029 ,31.783 ,19.2647 ,68.9094 ,39.4168 ,29.4874 ,19.9292 ,15.7162 ,41.5487 ,36.1616 ,44.148};

  char source_zfile[36][400]={"wl_redshift_model0_WLz01.880307e-01_WLalpha8.485694e-01.txt", "wl_redshift_model1_WLz01.731166e-01_WLalpha7.662434e-01.txt", "wl_redshift_model2_WLz01.846221e-01_WLalpha8.297540e-01.txt", "wl_redshift_model3_WLz01.758423e-01_WLalpha7.812893e-01.txt", "wl_redshift_model4_WLz01.721357e-01_WLalpha7.608290e-01.txt", "wl_redshift_model5_WLz01.931476e-01_WLalpha8.768147e-01.txt", "wl_redshift_model6_WLz01.919821e-01_WLalpha8.703814e-01.txt", "wl_redshift_model7_WLz01.751912e-01_WLalpha7.776952e-01.txt", "wl_redshift_model8_WLz01.741405e-01_WLalpha7.718958e-01.txt", "wl_redshift_model9_WLz01.860048e-01_WLalpha8.373865e-01.txt", "wl_redshift_model10_WLz01.888904e-01_WLalpha8.533151e-01.txt", "wl_redshift_model11_WLz01.954564e-01_WLalpha8.895592e-01.txt", "wl_redshift_model12_WLz01.945099e-01_WLalpha8.843347e-01.txt", "wl_redshift_model13_WLz01.853877e-01_WLalpha8.339802e-01.txt", "wl_redshift_model14_WLz01.775192e-01_WLalpha7.905458e-01.txt", "wl_redshift_model15_WLz01.925611e-01_WLalpha8.735775e-01.txt", "wl_redshift_model16_WLz01.708849e-01_WLalpha7.539247e-01.txt", "wl_redshift_model17_WLz01.733832e-01_WLalpha7.677150e-01.txt", "wl_redshift_model18_WLz01.820009e-01_WLalpha8.152851e-01.txt", "wl_redshift_model19_WLz01.767381e-01_WLalpha7.862345e-01.txt", "wl_redshift_model20_WLz01.866426e-01_WLalpha8.409073e-01.txt", "wl_redshift_model21_WLz01.877261e-01_WLalpha8.468882e-01.txt", "wl_redshift_model22_WLz01.826188e-01_WLalpha8.186960e-01.txt", "wl_redshift_model23_WLz01.768086e-01_WLalpha7.866234e-01.txt", "wl_redshift_model24_WLz01.802711e-01_WLalpha8.057362e-01.txt", "wl_redshift_model25_WLz01.897506e-01_WLalpha8.580632e-01.txt", "wl_redshift_model26_WLz01.831217e-01_WLalpha8.214720e-01.txt", "wl_redshift_model27_WLz01.906926e-01_WLalpha8.632630e-01.txt", "wl_redshift_model28_WLz01.714197e-01_WLalpha7.568766e-01.txt", "wl_redshift_model29_WLz01.798667e-01_WLalpha8.035040e-01.txt", "wl_redshift_model30_WLz01.842554e-01_WLalpha8.277299e-01.txt", "wl_redshift_model31_WLz01.901797e-01_WLalpha8.604322e-01.txt", "wl_redshift_model32_WLz01.937710e-01_WLalpha8.802561e-01.txt", "wl_redshift_model33_WLz01.790701e-01_WLalpha7.991072e-01.txt", "wl_redshift_model34_WLz01.811701e-01_WLalpha8.106987e-01.txt", "wl_redshift_model35_WLz01.781525e-01_WLalpha7.940419e-01.txt"};

  char lens_zfile[36][400]={"LSS_redshift_model0_LSSz02.629496e-01_LSSalpha9.285983e-01.txt","LSS_redshift_model1_LSSz02.852923e-01_LSSalpha8.985943e-01.txt","LSS_redshift_model2_LSSz02.664822e-01_LSSalpha9.186036e-01.txt","LSS_redshift_model3_LSSz02.798758e-01_LSSalpha9.014201e-01.txt","LSS_redshift_model4_LSSz02.873873e-01_LSSalpha8.978683e-01.txt","LSS_redshift_model5_LSSz02.593970e-01_LSSalpha9.470921e-01.txt","LSS_redshift_model6_LSSz02.600213e-01_LSSalpha9.425114e-01.txt","LSS_redshift_model7_LSSz02.811155e-01_LSSalpha9.006370e-01.txt","LSS_redshift_model8_LSSz02.831875e-01_LSSalpha8.995165e-01.txt","LSS_redshift_model9_LSSz02.649368e-01_LSSalpha9.224338e-01.txt","LSS_redshift_model10_LSSz02.622058e-01_LSSalpha9.314128e-01.txt","LSS_redshift_model11_LSSz02.584820e-01_LSSalpha9.568082e-01.txt","LSS_redshift_model12_LSSz02.588053e-01_LSSalpha9.527220e-01.txt","LSS_redshift_model13_LSSz02.656075e-01_LSSalpha9.206866e-01.txt","LSS_redshift_model14_LSSz02.768397e-01_LSSalpha9.037492e-01.txt","LSS_redshift_model15_LSSz02.596975e-01_LSSalpha9.447600e-01.txt","LSS_redshift_model16_LSSz02.901709e-01_LSSalpha8.971658e-01.txt","LSS_redshift_model17_LSSz02.847362e-01_LSSalpha8.988183e-01.txt","LSS_redshift_model18_LSSz02.698330e-01_LSSalpha9.121821e-01.txt","LSS_redshift_model19_LSSz02.782257e-01_LSSalpha9.026084e-01.txt","LSS_redshift_model20_LSSz02.642756e-01_LSSalpha9.243038e-01.txt","LSS_redshift_model21_LSSz02.632273e-01_LSSalpha9.276296e-01.txt","LSS_redshift_model22_LSSz02.689934e-01_LSSalpha9.135968e-01.txt","LSS_redshift_model23_LSSz02.780987e-01_LSSalpha9.027073e-01.txt","LSS_redshift_model24_LSSz02.723464e-01_LSSalpha9.085463e-01.txt","LSS_redshift_model25_LSSz02.615210e-01_LSSalpha9.343470e-01.txt","LSS_redshift_model26_LSSz02.683327e-01_LSSalpha9.147934e-01.txt","LSS_redshift_model27_LSSz02.608392e-01_LSSalpha9.376962e-01.txt","LSS_redshift_model28_LSSz02.889655e-01_LSSalpha8.974355e-01.txt","LSS_redshift_model29_LSSz02.729686e-01_LSSalpha9.077654e-01.txt","LSS_redshift_model30_LSSz02.669178e-01_LSSalpha9.176391e-01.txt","LSS_redshift_model31_LSSz02.612016e-01_LSSalpha9.358553e-01.txt","LSS_redshift_model32_LSSz02.591077e-01_LSSalpha9.496317e-01.txt","LSS_redshift_model33_LSSz02.742326e-01_LSSalpha9.063038e-01.txt","LSS_redshift_model34_LSSz02.710103e-01_LSSalpha9.103760e-01.txt","LSS_redshift_model35_LSSz02.757517e-01_LSSalpha9.047459e-01.txt"};

  char survey_designation[1][200]={"LSST"};
  char tomo_binning_source[1][200]={"source_std"};
  char tomo_binning_lens[1][200]={"LSST_gold"};
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // use this remapping only if files fail !!!!!!!!!
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // int fail[5]={1134, 259, 497, 623, 718};
    // hit=fail[hit-1];
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  int hit=atoi(argv[1]);
  Ntable.N_a=20;
  k=1;

  for(t=0;t<1;t++){

    //RUN MODE setup
    init_cosmo_runmode("emu");
    init_binning_fourier(15,20.0,3000.0,3000.0,21.0,10,10);
    init_survey(survey_designation[0],nsource_table[t],nlens_table[t],area_table[t]);
    sprintf(arg1,"zdistris/%s",source_zfile[t]);
    sprintf(arg2,"zdistris/%s",lens_zfile[t]);
    init_galaxies(arg1,arg2,"gaussian","gaussian",tomo_binning_source[0],tomo_binning_lens[0]);
    init_IA("none", "GAMA");
    init_probes("3x2pt");

    //set l-bins for shear, ggl, clustering, clusterWL
    double logdl=(log(like.lmax)-log(like.lmin))/like.Ncl;
    double *ell, *dell, *ell_Cluster, *dell_Cluster;
    ell=create_double_vector(0,like.Ncl-1);
    dell=create_double_vector(0,like.Ncl-1);
    ell_Cluster=create_double_vector(0,Cluster.lbin-1);
    dell_Cluster=create_double_vector(0,Cluster.lbin-1);
    int j=0;
    for(i=0;i<like.Ncl;i++){
      ell[i]=exp(log(like.lmin)+(i+0.5)*logdl);
      dell[i]=exp(log(like.lmin)+(i+1)*logdl)-exp(log(like.lmin)+(i*logdl));
      if(ell[i]<like.lmax_shear) printf("%le\n",ell[i]);
      if(ell[i]>like.lmax_shear){
        ell_Cluster[j]=ell[i];
        dell_Cluster[j]=dell[i];
        printf("%le %le\n",ell[i],ell_Cluster[j]);
        j++;
      }
    }

    printf("----------------------------------\n");
    survey.area=area_table[t];
    survey.n_gal=nsource_table[t];
    survey.n_lens=nlens_table[t];

    sprintf(survey.name,"%s_area%le_ng%le_nl%le",survey_designation[0],survey.area,survey.n_gal,survey.n_lens);
    printf("area: %le n_source: %le n_lens: %le\n",survey.area,survey.n_gal,survey.n_lens);

//    sprintf(covparams.outdir,"/home/u17/timeifler/covparallel/");
    sprintf(covparams.outdir,"/aurora_nobackup/cosmos/teifler/covparallel/");

    printf("----------------------------------\n");
    sprintf(OUTFILE,"%s_ssss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_shear_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }

    sprintf(OUTFILE,"%s_lsls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=l;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_ggl(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
       //printf("%d\n",k);
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_llll_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){ //auto bins only for now!
      for (m=l;m<tomo.clustering_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d %d %d\n",l,m,k);
      }
    }
    sprintf(OUTFILE,"%s_llss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }
    sprintf(OUTFILE,"%s_llls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering_ggl(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }
    sprintf(OUTFILE,"%s_lsss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_ggl_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }

    //******************************
    //******cluster covariance******
    //******************************

   //  sprintf(OUTFILE,"%s_nn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.cluster_Nbin; l++){
   //    for (m=0;m<tomo.cluster_Nbin; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //          run_cov_N_N (OUTFILE,covparams.outdir,l,m,k);
   //        }
   //      }
   //      k=k+1;
   //      //printf("%d\n",k);
   //    }
   //  }
   //  sprintf(OUTFILE,"%s_cscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.cgl_Npowerspectra; l++){
   //    for (m=0;m<tomo.cgl_Npowerspectra; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //          run_cov_cgl_cgl (OUTFILE,covparams.outdir,ell_Cluster,dell_Cluster,l,m,k);
   //        }
   //      }
   //      k=k+1;
   //    }
   //  }
   //  sprintf(OUTFILE,"%s_csn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.cgl_Npowerspectra; l++){
   //    for (m=0;m<tomo.cluster_Nbin; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //          run_cov_cgl_N (OUTFILE,covparams.outdir,ell_Cluster,dell_Cluster,l,m,k);
   //        }
   //      }
   //      k=k+1;
   //    }
   //  }
   //  //shear X cluster
   //  sprintf(OUTFILE,"%s_ssn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.shear_Npowerspectra; l++){
   //    for (m=0;m<tomo.cluster_Nbin; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //         run_cov_shear_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
   //       }
   //     }
   //     k=k+1;
   //   }
   // }
   // sprintf(OUTFILE,"%s_sscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   // for (l=0;l<tomo.shear_Npowerspectra; l++){
   //    for (m=0;m<tomo.cgl_Npowerspectra; m++){
   //      //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
   //        if(k==hit){
   //          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //          if (fopen(filename, "r") != NULL){exit(1);}
   //          else {
   //            run_cov_shear_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
   //          }
   //        }
   //        k=k+1;
   //      //}
   //    }
   //  }
   //  // ggl X cluster
   //  sprintf(OUTFILE,"%s_lsn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.ggl_Npowerspectra; l++){
   //    for (m=0;m<tomo.cluster_Nbin; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //         run_cov_ggl_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
   //        }
   //      }
   //      k=k+1;
   //    }
   //  }
   //  sprintf(OUTFILE,"%s_lscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.ggl_Npowerspectra; l++){
   //    for (m=0;m<tomo.cgl_Npowerspectra; m++){
   //      //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
   //        if(k==hit){
   //          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //          if (fopen(filename, "r") != NULL){exit(1);}
   //          else {
   //            run_cov_ggl_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
   //          }
   //        }
   //        k=k+1;
   //      //}
   //    }
   //  }
   //  // clustering X cluster
   //  sprintf(OUTFILE,"%s_lln_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.clustering_Npowerspectra; l++){
   //    for (m=0;m<tomo.cluster_Nbin; m++){
   //      if(k==hit){
   //        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //        if (fopen(filename, "r") != NULL){exit(1);}
   //        else {
   //         run_cov_cl_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
   //        }
   //      }
   //      //printf("%d\n",k);
   //      k=k+1;
   //    }
   //  }
   //  sprintf(OUTFILE,"%s_llcs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   //  for (l=0;l<tomo.clustering_Npowerspectra; l++){
   //    for (m=0;m<tomo.cgl_Npowerspectra; m++){
   //      //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
   //        if(k==hit){
   //          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
   //          if (fopen(filename, "r") != NULL){exit(1);}
   //          else {
   //            run_cov_cl_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
   //          }
   //        }
   //        k=k+1;
   //      //}
   //    }
   //  }
  }
  printf("number of cov blocks for parallelization: %d\n",k-1);
  printf("-----------------\n");
  printf("PROGRAM EXECUTED\n");
  printf("-----------------\n");
  return 0;
}


text:
#ifndef petscutils_h
#define petscutils_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>

CeedMemType MemTypeP2C(PetscMemType mtype);
PetscErrorCode ProjectToUnitSphere(DM dm);
PetscErrorCode Kershaw(DM dm_orig, PetscScalar eps);
typedef PetscErrorCode (*BCFunction)(PetscInt dim, PetscReal time,
                                     const PetscReal x[],
                                     PetscInt num_comp_u, PetscScalar *u, void *ctx);
PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt num_comp_u,
                               PetscInt topo_dim,
                               bool enforce_bc,  BCFunction bc_func);
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P,
    CeedInt topo_dim, CeedInt height, DMLabel domain_label, CeedInt value,
    CeedElemRestriction *elem_restr);

#endif // petscutils_h

ext: h

max_stars_repo_name: LeilaGhaffari/libCEED

max_issues_repo_name: LeilaGhaffari/libCEED

max_forks_repo_name: LeilaGhaffari/libCEED


text:
/* filter/test_gaussian.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.
 */

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* compute Gaussian filter by explicitely constructing window and computing weighted sum */
int
slow_gaussian(const gsl_filter_end_t etype, const double alpha, const size_t order, const gsl_vector * x,
              gsl_vector * y, const size_t K)
{
  const size_t n = x->size;
  const size_t H = K / 2;
  double *window = malloc(K * sizeof(double));
  double *kernel = malloc(K * sizeof(double));
  gsl_vector_view k = gsl_vector_view_array(kernel, K);
  size_t i;

  gsl_filter_gaussian_kernel(alpha, order, 1, &k.vector);

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, H, window);
      double sum = 0.0;
      size_t j;

      for (j = 0; j < wsize; ++j)
        sum += window[j] * kernel[wsize - j - 1];

      gsl_vector_set(y, i, sum);
    }

  free(window);
  free(kernel);

  return GSL_SUCCESS;
}

static void
fdiff(const gsl_vector * x, gsl_vector * dx)
{
  const size_t N = x->size;
  size_t i;

  for (i = 1; i < N - 1; ++i)
    {
      double xm1 = gsl_vector_get(x, i - 1);
      double xp1 = gsl_vector_get(x, i + 1);
      gsl_vector_set(dx, i, 0.5 * (xp1 - xm1));
    }

  gsl_vector_set(dx, 0, gsl_vector_get(x, 1) - gsl_vector_get(x, 0));
  gsl_vector_set(dx, N - 1, gsl_vector_get(x, N - 1) - gsl_vector_get(x, N - 2));
}

static void
test_gaussian_kernel(const double alpha, const size_t K)
{
  const size_t max_order = 3;
  gsl_vector * kernel = gsl_vector_alloc(K);
  gsl_vector * deriv = gsl_vector_alloc(K);
  gsl_vector * deriv_fd = gsl_vector_alloc(K);
  char buf[2048];
  size_t order;

  gsl_filter_gaussian_kernel(alpha, 0, 0, kernel);

  for (order = 1; order <= max_order; ++order)
    {
      gsl_filter_gaussian_kernel(alpha, order, 0, deriv);
      fdiff(kernel, deriv_fd);

      sprintf(buf, "gaussian kernel order=%zu alpha=%g K=%zu", order, alpha, K);
      compare_vectors(1.0e-2, deriv_fd, deriv, buf);

      gsl_vector_memcpy(kernel, deriv);
    }

  gsl_vector_free(kernel);
  gsl_vector_free(deriv);
  gsl_vector_free(deriv_fd);
}

static void
test_gaussian_proc(const double tol, const double alpha, const size_t order, const size_t n, const size_t K,
                   const gsl_filter_end_t etype, gsl_rng * rng_p)
{
  gsl_filter_gaussian_workspace * w = gsl_filter_gaussian_alloc(K);
  gsl_vector * x = gsl_vector_alloc(n);
  gsl_vector * y = gsl_vector_alloc(n);
  gsl_vector * z = gsl_vector_alloc(n);
  char buf[2048];

  random_vector(x, rng_p);

  /* y = filter(x) with slow brute force method */
  slow_gaussian(etype, alpha, order, x, y, K);

  /* y = filter(x) with fast method */
  gsl_filter_gaussian(etype, alpha, order, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu K=%zu endtype=%u alpha=%g order=%zu gaussian random", n, K, etype, alpha, order);
  compare_vectors(tol, z, y, buf);

  /* z = filter(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_filter_gaussian(etype, alpha, order, z, z, w);

  sprintf(buf, "n=%zu K=%zu endtype=%u alpha=%g order=%zu gaussian random in-place", n, K, etype, alpha, order);
  compare_vectors(tol, z, y, buf);

  gsl_filter_gaussian_free(w);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
}

static void
test_gaussian_deriv(const double alpha, const size_t n, const size_t K)
{
#if 0
  const double f_low = 1.0;
  const double f_high = 50.0;
  const double gamma = 2.0 * M_PI / (n - 1.0);
  const double dt = 1.0 / (n - 1.0);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *dx = gsl_vector_alloc(n);
  gsl_vector *y1 = gsl_vector_alloc(n);
  gsl_vector *y2 = gsl_vector_alloc(n);
  gsl_filter_gaussian_workspace *w = gsl_filter_gaussian_alloc(K);
  size_t i;

  /* make input signal composed of two sine waves at different frequencies */
  for (i = 0; i < n; ++i)
    {
      double xi = sin(gamma * f_low * i) + sin(gamma * f_high * i);
      double dxi = gamma * f_low * cos(gamma * f_low * i) +
                   gamma * f_high * cos(gamma * f_high * i);

      gsl_vector_set(x, i, xi);
      gsl_vector_set(dx, i, dxi);
    }

  /* compute y1 = G * dx(t)/dt */
  gsl_filter_gaussian(alpha, 0, dx, y1, w);

  /* compute y2 = dG/dt * x(t) */
  gsl_filter_gaussian(alpha, 1, x, y2, w);

  for (i = 0; i < n; ++i)
    {
      printf("%zu %.12e %.12e %.12e %.12e\n",
             i,
             gsl_vector_get(x, i),
             gsl_vector_get(dx, i),
             gsl_vector_get(y1, i),
             gsl_vector_get(y2, i));
    }

  gsl_vector_free(x);
  gsl_vector_free(dx);
  gsl_vector_free(y1);
  gsl_vector_free(y2);
  gsl_filter_gaussian_free(w);
#endif
}

static void
test_gaussian(gsl_rng * r)
{
  const double tol = 1.0e-10;
  size_t order;

  test_gaussian_kernel(3.0, 2001);

  for (order = 0; order <= 3; ++order)
    {
      test_gaussian_proc(tol, 2.5, order, 1000, 21, GSL_FILTER_END_PADZERO, r);
      test_gaussian_proc(tol, 3.0, order, 500, 11, GSL_FILTER_END_PADZERO, r);
      test_gaussian_proc(tol, 1.0, order, 50, 101, GSL_FILTER_END_PADZERO, r);
      test_gaussian_proc(tol, 2.0, order, 50, 11, GSL_FILTER_END_PADZERO, r);

      test_gaussian_proc(tol, 2.5, order, 1000, 21, GSL_FILTER_END_PADVALUE, r);
      test_gaussian_proc(tol, 3.0, order, 500, 11, GSL_FILTER_END_PADVALUE, r);
      test_gaussian_proc(tol, 1.0, order, 50, 101, GSL_FILTER_END_PADVALUE, r);
      test_gaussian_proc(tol, 2.0, order, 50, 11, GSL_FILTER_END_PADVALUE, r);

      test_gaussian_proc(tol, 2.5, order, 1000, 21, GSL_FILTER_END_TRUNCATE, r);
      test_gaussian_proc(tol, 3.0, order, 500, 11, GSL_FILTER_END_TRUNCATE, r);
      test_gaussian_proc(tol, 1.0, order, 50, 101, GSL_FILTER_END_TRUNCATE, r);
      test_gaussian_proc(tol, 2.0, order, 50, 11, GSL_FILTER_END_TRUNCATE, r);
    }
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>

#include "common-driver-utils.h"
#include "pc_GtKG.h"
#include "pc_ScaledGtKG.h"

PetscErrorCode BSSCR_PetscExtStokesSolversInitialize( void )
{
    Stg_PCRegister( "gtkg", "Solvers/KSPSolvers/src/BSSCR", "BSSCR_PCCreate_GtKG", BSSCR_PCCreate_GtKG );
	//Stg_PCRegister( "bfbt", "Solvers/KSPSolvers/src/BSSCR", "BSSCR_PCCreate_GtKG", BSSCR_PCCreate_GtKG );
	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PetscExtStokesSolversFinalize( void )
{

#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=4))
//  PCFinalizePackage();
#else
//  PCRegisterDestroy();
#endif

	PetscFunctionReturn(0);
}


text:
//
//  Encoder.h
//  PretendToWork
//
//  Created by tqtifnypmb on 08/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include "../types.h"

#include <string>
#include <gsl/gsl>

namespace brick
{

struct ASCIIConverter {
    typedef std::string result_type;

    static detail::CodePointList encode(gsl::span<const char> bytes);
    static result_type decode(const detail::CodePointList& cplist);
};

//struct UTF8Converter {
//    typedef std::string result_type;
//
//    static detail::CodePointList encode(gsl::span<const char> bytes);
//    static result_type decode(const detail::CodePointList& cplist);
//    static size_t numOfLine(const detail::CodePointList& cplist);
//};
//
}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
#include <cblas.h>
#include <clapack.h>
main(int nargs, char **args)
{
   double A[1]={1.0}, b[1]={1.0};
   int ipiv[1];
   clapack_dgesv(CblasColMajor, 1, 1, A, 1, ipiv, b, 1);
}

max_stars_repo_name: kevleyski/math-atlas

max_issues_repo_name: kevleyski/math-atlas

max_forks_repo_name: kevleyski/math-atlas


text:
/* Copyright 2016, Alistair Boyle, 3-clause BSD License */
#include "config.h"
#include <stdio.h>
#include <stdlib.h> /* for exit */
#include <math.h>  /* for fabs */
#include <cblas.h> /* for cblas_dgemm */

int
main (void)
{
    double const a[] = { 0.11, 0.12, 0.13,
                         0.21, 0.22, 0.23
                       };
    double const b[] = { 1011.0, 1012.0,
                         1021.0, 1022.0,
                         1031.0, 1032.0
                       };
    double c[] = { 0.00, 0.00,
                   0.00, 0.00
                 };
    /* Compute C = a A B + c C */
    cblas_dgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans,
                 3, 2, 3, /* m, n, k, */
                 1.0, a, 3, b, 2, /* a, A, m, B, n, */
                 0.0, c, 2); /* c, C, k */
    /* A is m x k, B is k x n, C is m x n */
    /* matlab: expect = a*b */
    double const expect[] = { 367.7600, 368.1200,
                              674.0600, 674.7200
                            };
    printf("\tactual\t\t\texpect\n");
    printf ("[ %g, %g    ==  [ %g, %g\n", c[0], c[1], expect[0], expect[1]);
    printf ("  %g, %g ]        %g, %g ]\n", c[2], c[3], expect[2], expect[3]);
    int i;
    for(i = 0; i < 4; i++) {
        if ( fabs(expect[i] - c[i]) > 1e-12 ) {
            printf("FAIL\n");
            return 1;
        }
    }
    printf("PASS\n");
    exit(EXIT_SUCCESS);
}

max_stars_repo_name: boyle/zedhat

max_issues_repo_name: boyle/zedhat

max_forks_repo_name: boyle/zedhat


text:
//
// Created by Harold on 2020/9/16.
//

#ifndef M_MATH_M_CURVEFIT_H
#define M_MATH_M_CURVEFIT_H

#include <map>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <numeric>
#include <vector>
#include <cmath>

namespace M_MATH {
    class CurveFit{
    public:
        enum robust_type {
            DEFAULT,
            BISQUARE,
            CAUCHY,
            FAIR,
            HUBER,
            OLS,
            WELSCH
        };
    public:
        bool is_valid() const { return !m_coefficient.empty(); }
        bool linearfit(double const* x, double const* y, size_t n);
        bool polyfit(double const* x, double const* y, size_t length, unsigned poly_n);
        bool robustfit(double const* x, double const* y, size_t length, robust_type type, unsigned poly_n);

        double getY(double x) const;

        template<typename IT_X, typename IT_Y>
        void getYs(IT_X x_begin, IT_X x_end, IT_Y y_begin) const;

        double get_coefficient(size_t n) const;
        void set_coefficient(size_t n, double c) { m_coefficient[n] = c; };
        size_t get_coefficient_size() const { return m_coefficient.size(); };

        double get_slope() { return m_coefficient[1]; };
        double get_intercept() { return m_coefficient[0]; };
        double get_SSR() const { return m_ssr; };
        double get_SSE() const { return m_sse; };
        double get_SST() const { return m_sst; };
        double get_RMSE() const { return m_rmse; };
        double get_RSquare() const { return m_RSquare; };
        double get_Goodness() const { return m_goodness; };

        void clear() {
            m_coefficient.clear();
            m_err.clear();
        }

    private:
        std::map<size_t, double> m_coefficient;  // <index, fitting coefficient>, [0]: slope, [1] intercept
        std::map<size_t, double> m_err;
        double m_cov;  // covariance
        double m_ssr;  // sum of squares due to regression
        double m_sse;  // sum of squares error
        double m_sst;  // m_sst = m_ssr + m_sse
        double m_rmse;  // RMS error
        double m_RSquare;  // R-Square

        double m_sumsq;  // sum of squares
        double m_goodness;  // fit goodness

        static void get_determinate_parameters(double const* y, double const* yf,
                                               size_t length,
                                               double &ssr,
                                               double &sse,
                                               double &sst,
                                               double &rmse,
                                               double &RSquare) {
            double y_mean = std::accumulate(y, y+length, 0.0);
            ssr = 0.0;
            for (auto i = 0; i < length; ++i) {
                ssr += (yf[i] - y_mean) * (yf[i] - y_mean);
                sse += (y[i] - yf[i]) * (y[i] - yf[i]);
            }
            sst = ssr + sse;
            rmse = std::sqrt(sse/(double)length);
            RSquare = 1.0 - ssr/sst;
        }

        static gsl_multifit_robust_type const* get_robust_type(robust_type type) {
            switch (type) {
                case DEFAULT:
                    return gsl_multifit_robust_default;
                case BISQUARE:
                    return gsl_multifit_robust_bisquare;
                case CAUCHY:
                    return gsl_multifit_robust_cauchy;
                case FAIR:
                    return gsl_multifit_robust_fair;
                case HUBER:
                    return gsl_multifit_robust_huber;
                case OLS:
                    return gsl_multifit_robust_ols;
                case WELSCH:
                    return gsl_multifit_robust_welsch;
            }
            return nullptr;
        }
    };

    inline bool CurveFit::linearfit(const double *x, const double *y, size_t n) {
        clear();
        m_coefficient[0] = 0;
        m_coefficient[1] = 1;
        m_err[0] = 0;
        m_err[0] = 0;
        // double const* x, size_t const x_stride,
        // double const* y, size_t const y_stride,
        // size_t n,
        // double &intercept,
        // double &slope,
        // double &intercept_err,
        // double &slope_err,
        // double &cov,
        // double &wssr
        int ret = gsl_fit_linear(x, 1, y, 1, n,
                                 &m_coefficient[0], &m_coefficient[1],
                                 &m_err[0], &m_err[1],
                                 &m_cov, &m_sumsq);
        if (ret != 0)
            return false;
        // gammaQ function
        m_goodness = gsl_cdf_chisq_Q(m_sumsq / 2.0, double(n - 2) / 2.0);
        {
            std::vector<double> yf(n, 0);
            getYs(x, x+n, yf.begin());
            get_determinate_parameters(y, &yf[0], n, m_ssr, m_sse, m_sst, m_rmse, m_RSquare);
        }
        return true;
    }

    inline bool CurveFit::polyfit(const double *x, const double *y, size_t length, unsigned poly_n) {
        gsl_matrix *Xf = gsl_matrix_alloc(length, poly_n+1);
        gsl_vector *Yf = gsl_vector_alloc(length);
        gsl_vector *c = gsl_vector_alloc(poly_n+1);
        gsl_matrix *cov = gsl_matrix_alloc(poly_n+1, poly_n+1);

        for (auto i = 0; i < length; i++) {
            gsl_matrix_set(Xf, i, 0, 1.0);
            gsl_vector_set(Yf, i, y[i]);
            for (auto j = 1; j < poly_n + 1; j++)
                gsl_matrix_set(Xf, i, j, std::pow(x[i], int(j)));
        }

        auto *workspace = gsl_multifit_linear_alloc(length, poly_n+1);
        double chisq;
        // temporary vec to hold computed coefficients
        std::vector<double> coe(poly_n+1, 0);

        int ret = gsl_multifit_linear(Xf, Yf, c, cov, &chisq, workspace);
        gsl_multifit_linear_free(workspace);
        for (auto i = 0; i < c->size; i++)
            coe[i] = gsl_vector_get(c, i);

        gsl_matrix_free(Xf);
        gsl_vector_free(Yf);
        gsl_vector_free(c);
        gsl_matrix_free(cov);

        // if not successfully computed, not change former results
        if (ret != 0)
            return false;

        // update results
        m_goodness = gsl_cdf_chisq_Q(chisq/2.0, double(length-2)/2.0);
        clear();
        for (auto i = 0; i < poly_n+1; i++)
            m_coefficient[i] = coe[i];
        {
            std::vector<double> yf(length, 0);
            getYs(x, x+length, yf.begin());
            get_determinate_parameters(y, &yf[0], length, m_ssr, m_sse, m_sst, m_rmse, m_RSquare);
        }
        return true;
    }

    inline double CurveFit::getY(double x) const {
        double ret = 0;
        for (auto const& it : m_coefficient)
            ret += (it.second) * std::pow(x, it.first);
        return ret;
    }

    template<typename IT_X, typename IT_Y>
    inline void CurveFit::getYs(IT_X x_begin, IT_X x_end, IT_Y y_begin) const {
        for (; x_begin != x_end; ++x_begin, ++y_begin)
            *y_begin = getY(*x_begin);
    }

    inline double CurveFit::get_coefficient(size_t n) const {
        auto it = m_coefficient.find(n);
        if (it != m_coefficient.end())
            return it->second;
        return 0.0;
    }

    // poly_n >= 2
    // poly_n = 2: linear fit
    bool CurveFit::robustfit(const double *x, const double *y, size_t length, robust_type type, unsigned int poly_n) {
        gsl_matrix *Xf = gsl_matrix_alloc(length, poly_n);
        gsl_vector *Yf = gsl_vector_alloc(length);
        gsl_vector *c = gsl_vector_alloc(poly_n);
        gsl_matrix *cov = gsl_matrix_alloc(poly_n, poly_n);

        for (auto i = 0; i < length; i++) {
            gsl_matrix_set(Xf, i, 0, 1.0);
            gsl_vector_set(Yf, i, y[i]);
            for (auto j = 1; j < poly_n; j++)
                gsl_matrix_set(Xf, i, j, std::pow(x[i], int(j)));
        }

        auto *workspace = gsl_multifit_robust_alloc(get_robust_type(type), length, poly_n);
        // temporary vec to hold computed coefficients
        std::vector<double> coe(poly_n, 0);

        int ret = gsl_multifit_robust(Xf, Yf, c, cov, workspace);
        gsl_multifit_robust_free(workspace);
        for (auto i = 0; i < c->size; i++)
            coe[i] = gsl_vector_get(c, i);

        gsl_matrix_free(Xf);
        gsl_vector_free(Yf);
        gsl_vector_free(c);
        gsl_matrix_free(cov);

        // if not successfully computed, not change former results
        if (ret != 0)
            return false;

        // update results
        clear();
        for (auto i = 0; i < poly_n; i++)
            m_coefficient[i] = coe[i];
        {
            std::vector<double> yf(length, 0);
            getYs(x, x+length, yf.begin());
            get_determinate_parameters(y, &yf[0], length, m_ssr, m_sse, m_sst, m_rmse, m_RSquare);
        }
        return true;
    }
}

#endif //M_MATH_M_CURVEFIT_H


text:
#pragma once

#include <gsl/gsl>

#include <memory>
#include <list>
#include <chrono>

#include <EASTL/fixed_vector.h>
#include <EASTL/map.h>

#include "../platform/d3d.h"

#include "buffers.h"
#include "shaders.h"
#include "camera.h"
#include "textures.h"

struct ID3D11Texture2D;

namespace gfx {

	class RenderingDevice;
	class Material;
	struct RasterizerSpec;
	class RasterizerState;
	struct BlendSpec;
	class BlendState;
	struct DepthStencilSpec;
	class DepthStencilState;
	struct SamplerSpec;
	class SamplerState;
	struct MaterialSamplerSpec;
	class BufferBinding;
	class TextEngine;

	using SamplerStatePtr = std::shared_ptr<SamplerState>;
	using DepthStencilStatePtr = std::shared_ptr<DepthStencilState>;
	using BlendStatePtr = std::shared_ptr<BlendState>;
	using RasterizerStatePtr = std::shared_ptr<RasterizerState>;

	using ResizeListener = std::function<void(int w, int h)>;

	// RAII style listener registration
	class ResizeListenerRegistration {
		friend class RenderingDevice;
	public:
		~ResizeListenerRegistration();

		NO_COPY(ResizeListenerRegistration);

		ResizeListenerRegistration(ResizeListenerRegistration &&o) : mDevice(o.mDevice), mKey(o.mKey) {
			o.mKey = 0;
		}
		ResizeListenerRegistration &operator =(ResizeListenerRegistration &&o) = delete;
	private:
		ResizeListenerRegistration(gfx::RenderingDevice &device, uint32_t key) : mDevice(device), mKey(key) {}

		gfx::RenderingDevice &mDevice;
		uint32_t mKey;
	};

	class ResourceListener {
	public:
		virtual ~ResourceListener();

		virtual void CreateResources(RenderingDevice&) = 0;
		virtual void FreeResources(RenderingDevice&) = 0;
	};

	enum class StandardSlotSemantic {
		ViewProjMatrix,
		UiProjMatrix
	};

	// RAII class for managing resource listener registrations
	class ResourceListenerRegistration {
	public:
		explicit ResourceListenerRegistration(RenderingDevice& device, ResourceListener* listener);
		~ResourceListenerRegistration();

		NO_COPY_OR_MOVE(ResourceListenerRegistration);

	private:
		RenderingDevice& mDevice;
		ResourceListener* mListener;
	};

	using DynamicTexturePtr = std::shared_ptr<class DynamicTexture>;
	using RenderTargetTexturePtr = std::shared_ptr<class RenderTargetTexture>;
	using RenderTargetDepthStencilPtr = std::shared_ptr<class RenderTargetDepthStencil>;

	// An output of a display device (think: monitor)
	struct DisplayDeviceOutput {
		// Technical id for use in a configuration or log file
		std::string id;
		// Name to display to the end user
		std::string name;
	};

	// A display device that can be used to render the game (think: GPU)
	struct DisplayDevice {
		// Technical id for use in a configuration or log file
		size_t id;
		// Name to display to the end user
		std::string name;
		// Outputs associated with this display device
		eastl::fixed_vector<DisplayDeviceOutput, 4> outputs;
	};

	enum class PrimitiveType {
		TriangleList,
		TriangleStrip,
		LineStrip,
		LineList,
		PointList,
	};

	enum class BufferFormat {
		A8,
		A8R8G8B8,
		X8R8G8B8
	};

	enum class MapMode {
		Read,
		Discard,
		NoOverwrite
	};

	template<typename TBuffer, typename TElement>
	class MappedBuffer {
	public:

		using View = gsl::span<TElement>;
		using Iterator = gsl::contiguous_span_iterator<View>;

		MappedBuffer(TBuffer &buffer,
			RenderingDevice &device,
			View data,
			uint32_t rowPitch) : mBuffer(buffer), mDevice(device), mData(data), mRowPitch(rowPitch) {}

		MappedBuffer(MappedBuffer&& o) : mBuffer(o.mBuffer), mDevice(o.mDevice), mData(o.mData), mRowPitch(o.rowPitch) {
			o.mMoved = true;
		}

		~MappedBuffer();

		Iterator begin() {
			return mData.begin();
		}

		Iterator end() {
			return mData.end();
		}

		TElement &operator[](size_t idx) {
			assert(!mMoved);
			return mData[idx];
		}

		size_t size() const {
			return mData.size();
		}

		TElement *GetData() const {
			return mData.data();
		}

		// Only relevant for mapped texture buffers
		size_t GetRowPitch() const {
			return mRowPitch;
		}

		void Unmap();

	private:
		TBuffer& mBuffer;
		RenderingDevice &mDevice;
		View mData;
		size_t mRowPitch;
		bool mMoved = false;
	};

	template<typename TElement>
	using MappedVertexBuffer = MappedBuffer<VertexBuffer, TElement>;
	using MappedIndexBuffer = MappedBuffer<IndexBuffer, uint16_t>;
	using MappedTexture = MappedBuffer<DynamicTexture, uint8_t>;

	class RenderingDevice {
	template<typename T>
	friend class Shader;
	friend class BufferBinding;
	friend class TextureLoader;
	friend class DebugUI;
	public:
		RenderingDevice(HWND mWindowHandle, uint32_t adapterIdx = 0, bool debugDevice = false);
		~RenderingDevice();

		bool BeginFrame();
		bool Present();
		void PresentForce();
		void Flush();

		void ClearCurrentColorTarget(XMCOLOR color);

		void ClearCurrentDepthTarget(bool clearDepth = true,
			bool clearStencil = true,
			float depthValue = 1.0f,
			uint8_t stencilValue = 0);

		using Clock = std::chrono::high_resolution_clock;
		using TimePoint = std::chrono::time_point<Clock>;

		TimePoint GetLastFrameStart() const {
			return mLastFrameStart;
		}
		TimePoint GetDeviceCreated() const {
			return mDeviceCreated;
		}

		const eastl::vector<DisplayDevice> &GetDisplayDevices();

		// Resize the back buffer
		void ResizeBuffers(int w, int h);

		Material CreateMaterial(
			const BlendSpec &blendSpec,
			const DepthStencilSpec &depthStencilSpec,
			const RasterizerSpec &rasterizerSpec,
			const std::vector<MaterialSamplerSpec> &samplerSpecs,
			const VertexShaderPtr &vs,
			const PixelShaderPtr &ps
		);

		BlendStatePtr CreateBlendState(const BlendSpec &spec);
		DepthStencilStatePtr CreateDepthStencilState(const DepthStencilSpec &spec);
		RasterizerStatePtr CreateRasterizerState(const RasterizerSpec &spec);
		SamplerStatePtr CreateSamplerState(const SamplerSpec &spec);

		// Changes the current scissor rect to the given rectangle
		void SetScissorRect(int x, int y, int width, int height);

		// Resets the scissor rect to the current render target's size
		void ResetScissorRect();

		std::shared_ptr<class IndexBuffer> CreateEmptyIndexBuffer(size_t count);
		std::shared_ptr<class VertexBuffer> CreateEmptyVertexBuffer(size_t count, bool forPoints = false);
		DynamicTexturePtr CreateDynamicTexture(gfx::BufferFormat format, int width, int height);
		DynamicTexturePtr CreateDynamicStagingTexture(gfx::BufferFormat format, int width, int height);
		void CopyRenderTarget(gfx::RenderTargetTexture &renderTarget, gfx::DynamicTexture &stagingTexture);
		RenderTargetTexturePtr CreateRenderTargetTexture(gfx::BufferFormat format, int width, int height, bool multiSampled = false);
		RenderTargetTexturePtr CreateRenderTargetForNativeSurface(ID3D11Texture2D *surface);
		RenderTargetTexturePtr CreateRenderTargetForSharedSurface(IUnknown *surface);
		RenderTargetDepthStencilPtr CreateRenderTargetDepthStencil(int width, int height, bool multiSampled = false);

		template<typename T>
		VertexBufferPtr CreateVertexBuffer(gsl::span<T> data, bool immutable = true);
		VertexBufferPtr CreateVertexBufferRaw(gsl::span<const uint8_t> data, bool immutable = true);
		IndexBufferPtr CreateIndexBuffer(gsl::span<const uint16_t> data, bool immutable = true);

		void SetMaterial(const Material &material);
		void SetVertexShaderConstant(uint32_t startRegister, StandardSlotSemantic semantic);
		void SetPixelShaderConstant(uint32_t startRegister, StandardSlotSemantic semantic);

		void SetRasterizerState(const RasterizerState &state);
		void SetBlendState(const BlendState &state);
		void SetDepthStencilState(const DepthStencilState &state);
		void SetSamplerState(int samplerIdx, const SamplerState &state);

		void SetTexture(uint32_t slot, gfx::Texture &texture);

		void SetIndexBuffer(const gfx::IndexBuffer &indexBuffer);

		void Draw(PrimitiveType type, uint32_t vertexCount, uint32_t startVertex = 0);

		void DrawIndexed(PrimitiveType type, uint32_t vertexCount, uint32_t indexCount, uint32_t startVertex = 0, uint32_t vertexBase = 0);

		/*
		  Changes the currently used cursor to the given surface.
		 */
		void SetCursor(int hotspotX, int hotspotY, const gfx::TextureRef &texture);
		void ShowCursor();
		void HideCursor();

		/*
			Take a screenshot with the given size. The image will be stretched
			to the given size.
		*/
		void TakeScaledScreenshot(const std::string& filename, int width, int height, int quality = 90);

		// Creates a buffer binding for a MDF material that
		// is preinitialized with the correct shader
		BufferBinding CreateMdfBufferBinding();

		Shaders& GetShaders() {
			return mShaders;
		}

		Textures& GetTextures() {
			return mTextures;
		}

		WorldCamera& GetCamera() {
			return mCamera;
		}

		void SetAntiAliasing(bool enable, uint32_t samples, uint32_t quality);

		template<typename T>
		void UpdateBuffer(VertexBuffer &buffer, gsl::span<T> data) {
			UpdateBuffer(buffer, data.data(), data.size_bytes());
		}

		void UpdateBuffer(VertexBuffer &buffer, const void *data, size_t size);

		void UpdateBuffer(IndexBuffer &buffer, gsl::span<uint16_t> data);

		template<typename TElement>
		MappedVertexBuffer<TElement> Map(VertexBuffer &buffer, gfx::MapMode mode = gfx::MapMode::Discard) {
			auto data = MapVertexBufferRaw(buffer, mode);
			auto castData = gsl::span<TElement>((TElement*)data.data(), data.size() / sizeof(TElement));
			return MappedVertexBuffer<TElement>(buffer, *this, castData, 0);
		}
		void Unmap(VertexBuffer &buffer);

		// Index buffer memory mapping techniques
		MappedIndexBuffer Map(IndexBuffer &buffer, gfx::MapMode mode = gfx::MapMode::Discard);
		void Unmap(IndexBuffer &buffer);

		MappedTexture Map(DynamicTexture &texture, gfx::MapMode mode = gfx::MapMode::Discard);
		void Unmap(DynamicTexture &texture);

		static constexpr uint32_t MaxVsConstantBufferSize = 2048;

		template<typename T>
		void SetVertexShaderConstants(uint32_t slot, const T &buffer) {
			static_assert(sizeof(T) <= MaxVsConstantBufferSize, "Constant buffer exceeds maximum size");
			UpdateResource(mVsConstantBuffer, &buffer, sizeof(T));
			VSSetConstantBuffer(slot, mVsConstantBuffer);
		}

		static constexpr uint32_t MaxPsConstantBufferSize = 512;

		template<typename T>
		void SetPixelShaderConstants(uint32_t slot, const T &buffer) {
			static_assert(sizeof(T) <= MaxPsConstantBufferSize, "Constant buffer exceeds maximum size");
			UpdateResource(mPsConstantBuffer, &buffer, sizeof(T));
			PSSetConstantBuffer(slot, mPsConstantBuffer);
		}

		const gfx::Size &GetBackBufferSize() const;

		// Pushes the back buffer and it's depth buffer as the current render target
		void PushBackBufferRenderTarget();
		void PushRenderTarget(
			const gfx::RenderTargetTexturePtr &colorBuffer,
			const gfx::RenderTargetDepthStencilPtr &depthStencilBuffer
		);
		void PopRenderTarget();
		const gfx::RenderTargetTexturePtr &GetCurrentRederTargetColorBuffer() const {
			return mRenderTargetStack.back().colorBuffer;
		}
		const gfx::RenderTargetDepthStencilPtr &GetCurrentRenderTargetDepthStencilBuffer() const {
			return mRenderTargetStack.back().depthStencilBuffer;
		}

		ResizeListenerRegistration AddResizeListener(ResizeListener listener);

		bool IsDebugDevice() const;

		/**
		 * Emits the start of a rendering call group if the debug device is being used.
		 * This information can be used in the graphic debugger.
		 */
		template<typename... T>
		void BeginPerfGroup(const char *format, const T &... args) const {
			if (IsDebugDevice()) {
				BeginPerfGroupInternal(fmt::format(format, args...).c_str());
			}
		}
		/**
		 * Ends a previously started performance group.
		 */
		void EndPerfGroup() const;

		TextEngine& GetTextEngine() const;

	private:
		friend class ResourceListenerRegistration;
		friend class ResizeListenerRegistration;

		void BeginPerfGroupInternal(const char *msg) const;

		void RemoveResizeListener(uint32_t key);

		void AddResourceListener(ResourceListener* resourceListener);
		void RemoveResourceListener(ResourceListener* resourceListener);

		void UpdateResource(ID3D11Resource *resource, const void *data, size_t size);
		CComPtr<ID3D11Buffer> CreateConstantBuffer(const void *initialData, size_t initialDataSize);
		void VSSetConstantBuffer(uint32_t slot, ID3D11Buffer *buffer);
		void PSSetConstantBuffer(uint32_t slot, ID3D11Buffer *buffer);

		gsl::span<uint8_t> MapVertexBufferRaw(VertexBuffer &buffer, MapMode mode);

		CComPtr<IDXGIAdapter1> GetAdapter(size_t index);

		int mBeginSceneDepth = 0;

		HWND mWindowHandle;

		CComPtr<IDXGIFactory1> mDxgiFactory;

		// The DXGI adapter we use
		CComPtr<IDXGIAdapter1> mAdapter;

		// D3D11 device and related
		CComPtr<ID3D11Device> mD3d11Device;
		CComPtr<ID3D11Device1> mD3d11Device1;
		DXGI_SWAP_CHAIN_DESC mSwapChainDesc;
		CComPtr<IDXGISwapChain> mSwapChain;
		CComPtr<ID3D11DeviceContext> mContext;
		gfx::RenderTargetTexturePtr mBackBufferNew;
		gfx::RenderTargetDepthStencilPtr mBackBufferDepthStencil;

		struct RenderTarget {
			gfx::RenderTargetTexturePtr colorBuffer;
			gfx::RenderTargetDepthStencilPtr depthStencilBuffer;
		};
		eastl::fixed_vector<RenderTarget, 16> mRenderTargetStack;

		D3D_FEATURE_LEVEL mFeatureLevel = D3D_FEATURE_LEVEL_9_1;

		eastl::vector<DisplayDevice> mDisplayDevices;

		CComPtr<ID3D11Buffer> mVsConstantBuffer;
		CComPtr<ID3D11Buffer> mPsConstantBuffer;

		eastl::map<uint32_t, ResizeListener> mResizeListeners;
		uint32_t mResizeListenersKey = 0;

		std::list<ResourceListener*> mResourcesListeners;
		bool mResourcesCreated = false;

		TimePoint mLastFrameStart = Clock::now();
		TimePoint mDeviceCreated = Clock::now();

		size_t mUsedSamplers = 0;

		Shaders mShaders;
		Textures mTextures;
		WorldCamera mCamera;

		struct Impl;
		std::unique_ptr<Impl> mImpl;
	};

	template <typename T>
	VertexBufferPtr RenderingDevice::CreateVertexBuffer(gsl::span<T> data, bool immutable) {
		return CreateVertexBufferRaw(gsl::span(reinterpret_cast<const uint8_t*>(&data[0]), data.size_bytes()), immutable);
	}

	extern RenderingDevice *renderingDevice;

	template<typename TBuffer, typename TElement>
	inline MappedBuffer<TBuffer, TElement>::~MappedBuffer()
	{
		if (!mMoved) {
			mDevice.Unmap(mBuffer);
		}
	}

	template<typename TBuffer, typename TElement>
	inline void MappedBuffer<TBuffer, TElement>::Unmap()
	{
		assert(!mMoved);
		mMoved = true;
		mDevice.Unmap(mBuffer);
	}

	inline ResizeListenerRegistration::~ResizeListenerRegistration() {
		mDevice.RemoveResizeListener(mKey);
	}

	// RAII style demarcation of render call groups for performance debugging
	class PerfGroup {
	public:
		template<typename... T>
		PerfGroup(RenderingDevice &device, const char *format, T... args) : mDevice(device) {
			device.BeginPerfGroup(format, args...);
		}
		~PerfGroup() {
			mDevice.EndPerfGroup();
		}
	private:
		const RenderingDevice &mDevice;
	};

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_rot (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int N = 1;
   float c = 0.0f;
   float s = 0.0f;
   float X[] = { -0.314f };
   int incX = 1;
   float Y[] = { -0.406f };
   int incY = -1;
   float x_expected[] = { 0.0f };
   float y_expected[] = { 0.0f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 558)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 559)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.866025403784f;
   float s = 0.5f;
   float X[] = { -0.314f };
   int incX = 1;
   float Y[] = { -0.406f };
   int incY = -1;
   float x_expected[] = { -0.474932f };
   float y_expected[] = { -0.194606f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 560)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 561)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.0f;
   float s = -1.0f;
   float X[] = { -0.314f };
   int incX = 1;
   float Y[] = { -0.406f };
   int incY = -1;
   float x_expected[] = { 0.406f };
   float y_expected[] = { -0.314f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 562)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 563)");
     }
   };
  };

  {
   int N = 1;
   float c = -1.0f;
   float s = 0.0f;
   float X[] = { -0.314f };
   int incX = 1;
   float Y[] = { -0.406f };
   int incY = -1;
   float x_expected[] = { 0.314f };
   float y_expected[] = { 0.406f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 564)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 565)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = 0;
   double X[] = { -0.493 };
   int incX = 1;
   double Y[] = { -0.014 };
   int incY = -1;
   double x_expected[] = { 0.0 };
   double y_expected[] = { 0.0 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 566)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 567)");
     }
   };
  };

  {
   int N = 1;
   double c = 0.866025403784;
   double s = 0.5;
   double X[] = { -0.493 };
   int incX = 1;
   double Y[] = { -0.014 };
   int incY = -1;
   double x_expected[] = { -0.433950524066 };
   double y_expected[] = { 0.234375644347 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 568)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 569)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = -1;
   double X[] = { -0.493 };
   int incX = 1;
   double Y[] = { -0.014 };
   int incY = -1;
   double x_expected[] = { 0.014 };
   double y_expected[] = { -0.493 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 570)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 571)");
     }
   };
  };

  {
   int N = 1;
   double c = -1;
   double s = 0;
   double X[] = { -0.493 };
   int incX = 1;
   double Y[] = { -0.014 };
   int incY = -1;
   double x_expected[] = { 0.493 };
   double y_expected[] = { 0.014 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 572)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 573)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.0f;
   float s = 0.0f;
   float X[] = { -0.808f };
   int incX = -1;
   float Y[] = { -0.511f };
   int incY = 1;
   float x_expected[] = { 0.0f };
   float y_expected[] = { 0.0f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 574)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 575)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.866025403784f;
   float s = 0.5f;
   float X[] = { -0.808f };
   int incX = -1;
   float Y[] = { -0.511f };
   int incY = 1;
   float x_expected[] = { -0.955249f };
   float y_expected[] = { -0.038539f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 576)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 577)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.0f;
   float s = -1.0f;
   float X[] = { -0.808f };
   int incX = -1;
   float Y[] = { -0.511f };
   int incY = 1;
   float x_expected[] = { 0.511f };
   float y_expected[] = { -0.808f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 578)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 579)");
     }
   };
  };

  {
   int N = 1;
   float c = -1.0f;
   float s = 0.0f;
   float X[] = { -0.808f };
   int incX = -1;
   float Y[] = { -0.511f };
   int incY = 1;
   float x_expected[] = { 0.808f };
   float y_expected[] = { 0.511f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 580)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 581)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = 0;
   double X[] = { -0.176 };
   int incX = -1;
   double Y[] = { -0.165 };
   int incY = 1;
   double x_expected[] = { 0.0 };
   double y_expected[] = { 0.0 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 582)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 583)");
     }
   };
  };

  {
   int N = 1;
   double c = 0.866025403784;
   double s = 0.5;
   double X[] = { -0.176 };
   int incX = -1;
   double Y[] = { -0.165 };
   int incY = 1;
   double x_expected[] = { -0.234920471066 };
   double y_expected[] = { -0.0548941916244 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 584)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 585)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = -1;
   double X[] = { -0.176 };
   int incX = -1;
   double Y[] = { -0.165 };
   int incY = 1;
   double x_expected[] = { 0.165 };
   double y_expected[] = { -0.176 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 586)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 587)");
     }
   };
  };

  {
   int N = 1;
   double c = -1;
   double s = 0;
   double X[] = { -0.176 };
   int incX = -1;
   double Y[] = { -0.165 };
   int incY = 1;
   double x_expected[] = { 0.176 };
   double y_expected[] = { 0.165 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 588)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 589)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.0f;
   float s = 0.0f;
   float X[] = { -0.201f };
   int incX = -1;
   float Y[] = { 0.087f };
   int incY = -1;
   float x_expected[] = { 0.0f };
   float y_expected[] = { 0.0f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 590)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 591)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.866025403784f;
   float s = 0.5f;
   float X[] = { -0.201f };
   int incX = -1;
   float Y[] = { 0.087f };
   int incY = -1;
   float x_expected[] = { -0.130571f };
   float y_expected[] = { 0.175844f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 592)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 593)");
     }
   };
  };

  {
   int N = 1;
   float c = 0.0f;
   float s = -1.0f;
   float X[] = { -0.201f };
   int incX = -1;
   float Y[] = { 0.087f };
   int incY = -1;
   float x_expected[] = { -0.087f };
   float y_expected[] = { -0.201f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 594)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 595)");
     }
   };
  };

  {
   int N = 1;
   float c = -1.0f;
   float s = 0.0f;
   float X[] = { -0.201f };
   int incX = -1;
   float Y[] = { 0.087f };
   int incY = -1;
   float x_expected[] = { 0.201f };
   float y_expected[] = { -0.087f };
   cblas_srot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], flteps, "srot(case 596)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "srot(case 597)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = 0;
   double X[] = { -0.464 };
   int incX = -1;
   double Y[] = { 0.7 };
   int incY = -1;
   double x_expected[] = { 0.0 };
   double y_expected[] = { 0.0 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 598)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 599)");
     }
   };
  };

  {
   int N = 1;
   double c = 0.866025403784;
   double s = 0.5;
   double X[] = { -0.464 };
   int incX = -1;
   double Y[] = { 0.7 };
   int incY = -1;
   double x_expected[] = { -0.051835787356 };
   double y_expected[] = { 0.838217782649 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 600)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 601)");
     }
   };
  };

  {
   int N = 1;
   double c = 0;
   double s = -1;
   double X[] = { -0.464 };
   int incX = -1;
   double Y[] = { 0.7 };
   int incY = -1;
   double x_expected[] = { -0.7 };
   double y_expected[] = { -0.464 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 602)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 603)");
     }
   };
  };

  {
   int N = 1;
   double c = -1;
   double s = 0;
   double X[] = { -0.464 };
   int incX = -1;
   double Y[] = { 0.7 };
   int incY = -1;
   double x_expected[] = { 0.464 };
   double y_expected[] = { -0.7 };
   cblas_drot(N, X, incX, Y, incY, c, s);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], x_expected[i], dbleps, "drot(case 604)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "drot(case 605)");
     }
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#include <ceed.h>
#include <petsc.h>
#include "../problems/problems.h"

PetscErrorCode RegisterProblems(ProblemFunctions problem_functions) {
  PetscErrorCode ierr;

  PetscFunctionBegin;

  SOLIDS_PROBLEM_REGISTER(problem_functions, "Linear", ElasLinear, NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "SS-NH", ElasSSNH, NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "FSCurrent-NH1", ElasFSCurrentNH1,
                          NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "FSCurrent-NH2", ElasFSCurrentNH2,
                          NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "FSInitial-NH1", ElasFSInitialNH1,
                          NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "FSInitial-NH2", ElasFSInitialNH2,
                          NH);
  SOLIDS_PROBLEM_REGISTER(problem_functions, "FSInitial-MR1", ElasFSInitialMR1,
                          MR);

  PetscFunctionReturn(0);
};

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED