text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

#ifndef setupsphere_h
#define setupsphere_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include <stdbool.h>
#include <string.h>
#include "qfunctions/bps/bp1sphere.h"
#include "qfunctions/bps/bp2sphere.h"
#include "qfunctions/bps/bp3sphere.h"
#include "qfunctions/bps/bp4sphere.h"
#include "qfunctions/bps/common.h"

#if PETSC_VERSION_LT(3,14,0)
#  define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i)
#  define DMPlexRestoreClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexRestoreClosureIndices(a,b,c,d,f,g,i)
#endif

#if PETSC_VERSION_LT(3,14,0)
#  define DMPlexCreateSphereMesh(a,b,c,d,e) DMPlexCreateSphereMesh(a,b,c,e)
#endif

// -----------------------------------------------------------------------------
// PETSc Operator Structs
// -----------------------------------------------------------------------------

// Data for PETSc
typedef struct UserO_ *UserO;
struct UserO_ {
  MPI_Comm comm;
  DM dm;
  Vec Xloc, Yloc, diag;
  CeedVector xceed, yceed;
  CeedOperator op;
  Ceed ceed;
};

// Data for PETSc Interp/Restrict operators
typedef struct UserIR_ *UserIR;
struct UserIR_ {
  MPI_Comm comm;
  DM dmc, dmf;
  Vec Xloc, Yloc, mult;
  CeedVector ceedvecc, ceedvecf;
  CeedOperator op;
  Ceed ceed;
};

// -----------------------------------------------------------------------------
// libCEED Data Struct
// -----------------------------------------------------------------------------

// libCEED data struct for level
typedef struct CeedData_ *CeedData;
struct CeedData_ {
  Ceed ceed;
  CeedBasis basisx, basisu, basisctof;
  CeedElemRestriction Erestrictx, Erestrictu, Erestrictui, Erestrictqdi;
  CeedQFunction qf_apply;
  CeedOperator op_apply, op_restrict, op_interp;
  CeedVector qdata, xceed, yceed;
};

// -----------------------------------------------------------------------------
// Command Line Options
// -----------------------------------------------------------------------------

// Coarsening options
typedef enum {
  COARSEN_UNIFORM = 0, COARSEN_LOGARITHMIC = 1
} coarsenType;
static const char *const coarsenTypes [] = {"uniform","logarithmic",
                                            "coarsenType","COARSEN",0
                                           };

// -----------------------------------------------------------------------------
// BP Option Data
// -----------------------------------------------------------------------------

// BP options
typedef enum {
  CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2,
  CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5
} bpType;
static const char *const bpTypes[] = {"bp1","bp2","bp3","bp4","bp5","bp6",
                                      "bpType","CEED_BP",0
                                     };

// BP specific data
typedef struct {
  CeedInt ncompu, qdatasize, qextra;
  CeedQFunctionUser setupgeo, setuprhs, apply, error;
  const char *setupgeofname, *setuprhsfname, *applyfname, *errorfname;
  CeedEvalMode inmode, outmode;
  CeedQuadMode qmode;
  PetscBool enforce_bc;
  PetscErrorCode (*bcs_func)(PetscInt, PetscReal, const PetscReal *,
                             PetscInt, PetscScalar *, void *);
} bpData;

static bpData bpOptions[6] = {
  [CEED_BP1] = {
    .ncompu = 1,
    .qdatasize = 1,
    .qextra = 1,
    .setupgeo = SetupMassGeo,
    .setuprhs = SetupMassRhs,
    .apply = Mass,
    .error = Error,
    .setupgeofname = SetupMassGeo_loc,
    .setuprhsfname = SetupMassRhs_loc,
    .applyfname = Mass_loc,
    .errorfname = Error_loc,
    .inmode = CEED_EVAL_INTERP,
    .outmode = CEED_EVAL_INTERP,
    .qmode = CEED_GAUSS
  },
  [CEED_BP2] = {
    .ncompu = 3,
    .qdatasize = 1,
    .qextra = 1,
    .setupgeo = SetupMassGeo,
    .setuprhs = SetupMassRhs3,
    .apply = Mass3,
    .error = Error3,
    .setupgeofname = SetupMassGeo_loc,
    .setuprhsfname = SetupMassRhs3_loc,
    .applyfname = Mass3_loc,
    .errorfname = Error3_loc,
    .inmode = CEED_EVAL_INTERP,
    .outmode = CEED_EVAL_INTERP,
    .qmode = CEED_GAUSS
  },
  [CEED_BP3] = {
    .ncompu = 1,
    .qdatasize = 4,
    .qextra = 1,
    .setupgeo = SetupDiffGeo,
    .setuprhs = SetupDiffRhs,
    .apply = Diff,
    .error = Error,
    .setupgeofname = SetupDiffGeo_loc,
    .setuprhsfname = SetupDiffRhs_loc,
    .applyfname = Diff_loc,
    .errorfname = Error_loc,
    .inmode = CEED_EVAL_GRAD,
    .outmode = CEED_EVAL_GRAD,
    .qmode = CEED_GAUSS
  },
  [CEED_BP4] = {
    .ncompu = 3,
    .qdatasize = 4,
    .qextra = 1,
    .setupgeo = SetupDiffGeo,
    .setuprhs = SetupDiffRhs3,
    .apply = Diff3,
    .error = Error3,
    .setupgeofname = SetupDiffGeo_loc,
    .setuprhsfname = SetupDiffRhs3_loc,
    .applyfname = Diff_loc,
    .errorfname = Error3_loc,
    .inmode = CEED_EVAL_GRAD,
    .outmode = CEED_EVAL_GRAD,
    .qmode = CEED_GAUSS
  },
  [CEED_BP5] = {
    .ncompu = 1,
    .qdatasize = 4,
    .qextra = 0,
    .setupgeo = SetupDiffGeo,
    .setuprhs = SetupDiffRhs,
    .apply = Diff,
    .error = Error,
    .setupgeofname = SetupDiffGeo_loc,
    .setuprhsfname = SetupDiffRhs_loc,
    .applyfname = Diff_loc,
    .errorfname = Error_loc,
    .inmode = CEED_EVAL_GRAD,
    .outmode = CEED_EVAL_GRAD,
    .qmode = CEED_GAUSS_LOBATTO
  },
  [CEED_BP6] = {
    .ncompu = 3,
    .qdatasize = 4,
    .qextra = 0,
    .setupgeo = SetupDiffGeo,
    .setuprhs = SetupDiffRhs3,
    .apply = Diff3,
    .error = Error3,
    .setupgeofname = SetupDiffGeo_loc,
    .setuprhsfname = SetupDiffRhs3_loc,
    .applyfname = Diff_loc,
    .errorfname = Error3_loc,
    .inmode = CEED_EVAL_GRAD,
    .outmode = CEED_EVAL_GRAD,
    .qmode = CEED_GAUSS_LOBATTO
  }
};

// -----------------------------------------------------------------------------
// PETSc sphere auxiliary functions
// -----------------------------------------------------------------------------

// Utility function taken from petsc/src/dm/impls/plex/examples/tutorials/ex7.c
static PetscErrorCode ProjectToUnitSphere(DM dm) {
  Vec            coordinates;
  PetscScalar   *coords;
  PetscInt       Nv, v, dim, d;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr);
  ierr = VecGetLocalSize(coordinates, &Nv); CHKERRQ(ierr);
  ierr = VecGetBlockSize(coordinates, &dim); CHKERRQ(ierr);
  Nv  /= dim;
  ierr = VecGetArray(coordinates, &coords); CHKERRQ(ierr);
  for (v = 0; v < Nv; ++v) {
    PetscReal r = 0.0;

    for (d = 0; d < dim; ++d) r += PetscSqr(PetscRealPart(coords[v*dim+d]));
    r = PetscSqrtReal(r);
    for (d = 0; d < dim; ++d) coords[v*dim+d] /= r;
  }
  ierr = VecRestoreArray(coordinates, &coords); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// PETSc FE Boilerplate
// -----------------------------------------------------------------------------

// Create FE by degree
static int PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc,
                                 PetscBool isSimplex, const char prefix[],
                                 PetscInt order, PetscFE *fem) {
  PetscQuadrature q, fq;
  DM              K;
  PetscSpace      P;
  PetscDualSpace  Q;
  PetscInt        quadPointsPerEdge;
  PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
  PetscErrorCode  ierr;

  PetscFunctionBeginUser;
  /* Create space */
  ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix); CHKERRQ(ierr);
  ierr = PetscSpacePolynomialSetTensor(P, tensor); CHKERRQ(ierr);
  ierr = PetscSpaceSetFromOptions(P); CHKERRQ(ierr);
  ierr = PetscSpaceSetNumComponents(P, Nc); CHKERRQ(ierr);
  ierr = PetscSpaceSetNumVariables(P, dim); CHKERRQ(ierr);
  ierr = PetscSpaceSetDegree(P, order, order); CHKERRQ(ierr);
  ierr = PetscSpaceSetUp(P); CHKERRQ(ierr);
  ierr = PetscSpacePolynomialGetTensor(P, &tensor); CHKERRQ(ierr);
  /* Create dual space */
  ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
  CHKERRQ(ierr);
  ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); CHKERRQ(ierr);
  ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetDM(Q, K); CHKERRQ(ierr);
  ierr = DMDestroy(&K); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetNumComponents(Q, Nc); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetOrder(Q, order); CHKERRQ(ierr);
  ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetFromOptions(Q); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetUp(Q); CHKERRQ(ierr);
  /* Create element */
  ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); CHKERRQ(ierr);
  ierr = PetscFESetFromOptions(*fem); CHKERRQ(ierr);
  ierr = PetscFESetBasisSpace(*fem, P); CHKERRQ(ierr);
  ierr = PetscFESetDualSpace(*fem, Q); CHKERRQ(ierr);
  ierr = PetscFESetNumComponents(*fem, Nc); CHKERRQ(ierr);
  ierr = PetscFESetUp(*fem); CHKERRQ(ierr);
  ierr = PetscSpaceDestroy(&P); CHKERRQ(ierr);
  ierr = PetscDualSpaceDestroy(&Q); CHKERRQ(ierr);
  /* Create quadrature */
  quadPointsPerEdge = PetscMax(order + 1,1);
  if (isSimplex) {
    ierr = PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0,
                                          &q); CHKERRQ(ierr);
    ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0,
                                          &fq); CHKERRQ(ierr);
  } else {
    ierr = PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0,
                                        &q); CHKERRQ(ierr);
    ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0,
                                        &fq); CHKERRQ(ierr);
  }
  ierr = PetscFESetQuadrature(*fem, q); CHKERRQ(ierr);
  ierr = PetscFESetFaceQuadrature(*fem, fq); CHKERRQ(ierr);
  ierr = PetscQuadratureDestroy(&q); CHKERRQ(ierr);
  ierr = PetscQuadratureDestroy(&fq); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// PETSc Setup for Level
// -----------------------------------------------------------------------------

// This function sets up a DM for a given degree
static int SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompu,
                           PetscInt dim) {
  PetscInt ierr;
  PetscFE fe;

  PetscFunctionBeginUser;

  // Setup FE
  ierr = PetscFECreateByDegree(dm, dim, ncompu, PETSC_FALSE, NULL, degree, &fe);
  CHKERRQ(ierr);
  ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr);

  // Setup DM
  ierr = DMCreateDS(dm); CHKERRQ(ierr);
  ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);
  CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// libCEED Setup for Level
// -----------------------------------------------------------------------------

// Destroy libCEED operator objects
static PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) {
  PetscInt ierr;

  CeedVectorDestroy(&data->qdata);
  CeedVectorDestroy(&data->xceed);
  CeedVectorDestroy(&data->yceed);
  CeedBasisDestroy(&data->basisx);
  CeedBasisDestroy(&data->basisu);
  CeedElemRestrictionDestroy(&data->Erestrictu);
  CeedElemRestrictionDestroy(&data->Erestrictx);
  CeedElemRestrictionDestroy(&data->Erestrictui);
  CeedElemRestrictionDestroy(&data->Erestrictqdi);
  CeedQFunctionDestroy(&data->qf_apply);
  CeedOperatorDestroy(&data->op_apply);
  if (i > 0) {
    CeedOperatorDestroy(&data->op_interp);
    CeedBasisDestroy(&data->basisctof);
    CeedOperatorDestroy(&data->op_restrict);
  }
  ierr = PetscFree(data); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Get CEED restriction data from DMPlex
static int CreateRestrictionPlex(Ceed ceed, CeedInt P, CeedInt ncomp,
                                 CeedElemRestriction *Erestrict, DM dm) {
  PetscInt ierr;
  PetscInt c, cStart, cEnd, nelem, nnodes, *erestrict, eoffset;
  PetscSection section;
  Vec Uloc;

  PetscFunctionBegin;

  // Get Nelem
  ierr = DMGetSection(dm, &section); CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd); CHKERRQ(ierr);
  nelem = cEnd - cStart;

  // Get indices
  ierr = PetscMalloc1(nelem*P*P, &erestrict); CHKERRQ(ierr);
  for (c=cStart, eoffset = 0; c<cEnd; c++) {
    PetscInt numindices, *indices, i;
    ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE,
                                   &numindices, &indices, NULL, NULL);
    CHKERRQ(ierr);
    for (i=0; i<numindices; i+=ncomp) {
      for (PetscInt j=0; j<ncomp; j++) {
        if (indices[i+j] != indices[i] + copysign(j, indices[i]))
          SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
                   "Cell %D closure indices not interlaced", c);
      }
      // NO BC on closed surfaces
      PetscInt loc = indices[i];
      erestrict[eoffset++] = loc;
    }
    ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE,
                                       &numindices, &indices, NULL, NULL);
    CHKERRQ(ierr);
  }

  // Setup CEED restriction
  ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr);
  ierr = VecGetLocalSize(Uloc, &nnodes); CHKERRQ(ierr);

  ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr);
  CeedElemRestrictionCreate(ceed, nelem, P*P, ncomp, 1, nnodes, CEED_MEM_HOST,
                            CEED_COPY_VALUES, erestrict, Erestrict);
  ierr = PetscFree(erestrict); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Set up libCEED for a given degree
static int SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
                                CeedInt topodim, CeedInt qextra,
                                PetscInt ncompx, PetscInt ncompu,
                                PetscInt gsize, PetscInt xlsize,
                                bpType bpChoice, CeedData data,
                                PetscBool setup_rhs, CeedVector rhsceed,
                                CeedVector *target) {
  int ierr;
  DM dmcoord;
  Vec coords;
  const PetscScalar *coordArray;
  CeedBasis basisx, basisu;
  CeedElemRestriction Erestrictx, Erestrictu, Erestrictui, Erestrictqdi;
  CeedQFunction qf_setupgeo, qf_apply;
  CeedOperator op_setupgeo, op_apply;
  CeedVector xcoord, qdata, xceed, yceed;
  CeedInt P, Q, cStart, cEnd, nelem, qdatasize = bpOptions[bpChoice].qdatasize;
  CeedScalar R = 1,                      // radius of the sphere
             l = 1.0/PetscSqrtReal(3.0); // half edge of the inscribed cube

  // CEED bases
  P = degree + 1;
  Q = P + qextra;
  CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompu, P, Q,
                                  bpOptions[bpChoice].qmode, &basisu);
  CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompx, 2, Q,
                                  bpOptions[bpChoice].qmode, &basisx);

  // CEED restrictions
  ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr);
  ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL);
  CHKERRQ(ierr);

  ierr = CreateRestrictionPlex(ceed, 2, ncompx, &Erestrictx, dmcoord);
  CHKERRQ(ierr);
  ierr = CreateRestrictionPlex(ceed, P, ncompu, &Erestrictu, dm); CHKERRQ(ierr);

  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
  nelem = cEnd - cStart;

  CeedElemRestrictionCreateStrided(ceed, nelem, Q*Q, ncompu, ncompu*nelem*Q*Q,
                                   CEED_STRIDES_BACKEND, &Erestrictui);
  CeedElemRestrictionCreateStrided(ceed, nelem, Q*Q, qdatasize,
                                   qdatasize*nelem*Q*Q,
                                   CEED_STRIDES_BACKEND, &Erestrictqdi);

  // Element coordinates
  ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr);
  ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr);

  CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);
  CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES,
                     (PetscScalar *)coordArray);
  ierr = VecRestoreArrayRead(coords, &coordArray);

  // Create the persistent vectors that will be needed in setup and apply
  CeedInt nqpts;
  CeedBasisGetNumQuadraturePoints(basisu, &nqpts);
  CeedVectorCreate(ceed, qdatasize*nelem*nqpts, &qdata);
  CeedVectorCreate(ceed, xlsize, &xceed);
  CeedVectorCreate(ceed, xlsize, &yceed);

  // Create the Q-function that builds the operator (i.e. computes its
  // quadrature data) and set its context data
  CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].setupgeo,
                              bpOptions[bpChoice].setupgeofname, &qf_setupgeo);
  CeedQFunctionAddInput(qf_setupgeo, "x", ncompx, CEED_EVAL_INTERP);
  CeedQFunctionAddInput(qf_setupgeo, "dx", ncompx*topodim, CEED_EVAL_GRAD);
  CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT);
  CeedQFunctionAddOutput(qf_setupgeo, "qdata", qdatasize, CEED_EVAL_NONE);

  // Set up PDE operator
  CeedInt inscale = bpOptions[bpChoice].inmode==CEED_EVAL_GRAD ? topodim : 1;
  CeedInt outscale = bpOptions[bpChoice].outmode==CEED_EVAL_GRAD ? topodim : 1;
  CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].apply,
                              bpOptions[bpChoice].applyfname, &qf_apply);
  CeedQFunctionAddInput(qf_apply, "u", ncompu*inscale,
                        bpOptions[bpChoice].inmode);
  CeedQFunctionAddInput(qf_apply, "qdata", qdatasize, CEED_EVAL_NONE);
  CeedQFunctionAddOutput(qf_apply, "v", ncompu*outscale,
                         bpOptions[bpChoice].outmode);

  // Create the operator that builds the quadrature data for the operator
  CeedOperatorCreate(ceed, qf_setupgeo, NULL, NULL, &op_setupgeo);
  CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx,
                       CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx,
                       CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx,
                       CEED_VECTOR_NONE);
  CeedOperatorSetField(op_setupgeo, "qdata", Erestrictqdi,
                       CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

  // Create the mass or diff operator
  CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply);
  CeedOperatorSetField(op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_apply, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED,
                       qdata);
  CeedOperatorSetField(op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE);

  // Setup qdata
  CeedOperatorApply(op_setupgeo, xcoord, qdata, CEED_REQUEST_IMMEDIATE);

  // Set up RHS if needed
  if (setup_rhs) {
    CeedQFunction qf_setuprhs;
    CeedOperator op_setuprhs;
    CeedVectorCreate(ceed, nelem*nqpts*ncompu, target);

    // Create the q-function that sets up the RHS and true solution
    CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].setuprhs,
                                bpOptions[bpChoice].setuprhsfname, &qf_setuprhs);
    CeedQFunctionAddInput(qf_setuprhs, "x", ncompx, CEED_EVAL_INTERP);
    CeedQFunctionAddInput(qf_setuprhs, "qdata", qdatasize, CEED_EVAL_NONE);
    CeedQFunctionAddOutput(qf_setuprhs, "true_soln", ncompu, CEED_EVAL_NONE);
    CeedQFunctionAddOutput(qf_setuprhs, "rhs", ncompu, CEED_EVAL_INTERP);

    // Create the operator that builds the RHS and true solution
    CeedOperatorCreate(ceed, qf_setuprhs, NULL, NULL, &op_setuprhs);
    CeedOperatorSetField(op_setuprhs, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE);
    CeedOperatorSetField(op_setuprhs, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED,
                         qdata);
    CeedOperatorSetField(op_setuprhs, "true_soln", Erestrictui,
                         CEED_BASIS_COLLOCATED,
                         *target);
    CeedOperatorSetField(op_setuprhs, "rhs", Erestrictu, basisu,
                         CEED_VECTOR_ACTIVE);

    // Set up the libCEED context
    CeedQFunctionContext rhsSetup;
    CeedQFunctionContextCreate(ceed, &rhsSetup);
    CeedScalar rhsSetupData[2] = {R, l};
    CeedQFunctionContextSetData(rhsSetup, CEED_MEM_HOST, CEED_COPY_VALUES,
                                sizeof rhsSetupData, &rhsSetupData);
    CeedQFunctionSetContext(qf_setuprhs, rhsSetup);
    CeedQFunctionContextDestroy(&rhsSetup);

    // Setup RHS and target
    CeedOperatorApply(op_setuprhs, xcoord, rhsceed, CEED_REQUEST_IMMEDIATE);
    CeedVectorTakeArray(rhsceed, CEED_MEM_HOST, NULL);

    // Cleanup
    CeedQFunctionDestroy(&qf_setuprhs);
    CeedOperatorDestroy(&op_setuprhs);
  }

  // Cleanup
  CeedQFunctionDestroy(&qf_setupgeo);
  CeedOperatorDestroy(&op_setupgeo);
  CeedVectorDestroy(&xcoord);

  // Save libCEED data required for level
  data->basisx = basisx; data->basisu = basisu;
  data->Erestrictx = Erestrictx;
  data->Erestrictu = Erestrictu;
  data->Erestrictui = Erestrictui;
  data->Erestrictqdi = Erestrictqdi;
  data->qf_apply = qf_apply;
  data->op_apply = op_apply;
  data->qdata = qdata;
  data->xceed = xceed;
  data->yceed = yceed;

  PetscFunctionReturn(0);
}

#ifdef multigrid
// Setup libCEED level transfer operator objects
static PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt numlevels,
    CeedInt ncompu, bpType bpChoice, CeedData *data, CeedInt *leveldegrees,
    CeedQFunction qf_restrict, CeedQFunction qf_prolong) {
  // Return early if numlevels=1
  if (numlevels==1)
    PetscFunctionReturn(0);

  // Set up each level
  for (CeedInt i=1; i<numlevels; i++) {
    // P coarse and P fine
    CeedInt Pc = leveldegrees[i-1] + 1;
    CeedInt Pf = leveldegrees[i] + 1;

    // Restriction - Fine to corse
    CeedBasis basisctof;
    CeedOperator op_restrict;

    // Basis
    CeedBasisCreateTensorH1Lagrange(ceed, 3, ncompu, Pc, Pf,
                                    CEED_GAUSS_LOBATTO, &basisctof);

    // Create the restriction operator
    CeedOperatorCreate(ceed, qf_restrict, NULL, NULL, &op_restrict);
    CeedOperatorSetField(op_restrict, "input", data[i]->Erestrictu,
                         CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
    CeedOperatorSetField(op_restrict, "output", data[i-1]->Erestrictu,
                         basisctof, CEED_VECTOR_ACTIVE);

    // Save libCEED data required for level
    data[i]->basisctof = basisctof;
    data[i]->op_restrict = op_restrict;

    // Interpolation - Corse to fine
    CeedOperator op_interp;

    // Create the prolongation operator
    CeedOperatorCreate(ceed, qf_prolong, NULL, NULL, &op_interp);
    CeedOperatorSetField(op_interp, "uin", data[i-1]->Erestrictu,
                         basisctof, CEED_VECTOR_ACTIVE);
    CeedOperatorSetField(op_interp, "uout", data[i]->Erestrictu,
                         CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

    // Save libCEED data required for level
    data[i]->op_interp = op_interp;
  }

  PetscFunctionReturn(0);
}
#endif

// -----------------------------------------------------------------------------
// Mat Shell Functions
// -----------------------------------------------------------------------------

#ifdef multigrid
// This function returns the computed diagonal of the operator
static PetscErrorCode MatGetDiag(Mat A, Vec D) {
  PetscErrorCode ierr;
  UserO user;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);

  ierr = VecCopy(user->diag, D); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}
#endif

// This function uses libCEED to compute the action of the Laplacian with
// Dirichlet boundary conditions
static PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) {
  PetscErrorCode ierr;
  UserO user;
  PetscScalar *x, *y;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);

  // Global-to-local
  ierr = DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr);

  // Setup CEED vectors
  ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr);
  ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr);
  CeedVectorSetArray(user->xceed, CEED_MEM_HOST, CEED_USE_POINTER, x);
  CeedVectorSetArray(user->yceed, CEED_MEM_HOST, CEED_USE_POINTER, y);

  // Apply CEED operator
  CeedOperatorApply(user->op, user->xceed, user->yceed, CEED_REQUEST_IMMEDIATE);

  // Restore PETSc vectors
  CeedVectorTakeArray(user->xceed, CEED_MEM_HOST, NULL);
  CeedVectorTakeArray(user->yceed, CEED_MEM_HOST, NULL);
  ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
  CHKERRQ(ierr);
  ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr);

  // Local-to-global
  ierr = VecZeroEntries(Y); CHKERRQ(ierr);
  ierr = DMLocalToGlobalBegin(user->dm, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);
  ierr = DMLocalToGlobalEnd(user->dm, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#ifdef multigrid
// This function uses libCEED to compute the action of the interp operator
static PetscErrorCode MatMult_Interp(Mat A, Vec X, Vec Y) {
  PetscErrorCode ierr;
  UserIR user;
  PetscScalar *x, *y;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);

  // Global-to-local
  ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(user->dmc, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(user->dmc, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr);

  // Setup CEED vectors
  ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr);
  ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr);
  CeedVectorSetArray(user->ceedvecc, CEED_MEM_HOST, CEED_USE_POINTER, x);
  CeedVectorSetArray(user->ceedvecf, CEED_MEM_HOST, CEED_USE_POINTER, y);

  // Apply CEED operator
  CeedOperatorApply(user->op, user->ceedvecc, user->ceedvecf,
                    CEED_REQUEST_IMMEDIATE);

  // Restore PETSc vectors
  CeedVectorTakeArray(user->ceedvecc, CEED_MEM_HOST, NULL);
  CeedVectorTakeArray(user->ceedvecf, CEED_MEM_HOST, NULL);
  ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
  CHKERRQ(ierr);
  ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr);

  // Multiplicity
  ierr = VecPointwiseMult(user->Yloc, user->Yloc, user->mult);

  // Local-to-global
  ierr = VecZeroEntries(Y); CHKERRQ(ierr);
  ierr = DMLocalToGlobalBegin(user->dmf, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);
  ierr = DMLocalToGlobalEnd(user->dmf, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// This function uses libCEED to compute the action of the restriction operator
static PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) {
  PetscErrorCode ierr;
  UserIR user;
  PetscScalar *x, *y;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);

  // Global-to-local
  ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(user->dmf, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(user->dmf, X, INSERT_VALUES, user->Xloc);
  CHKERRQ(ierr);
  ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr);

  // Multiplicity
  ierr = VecPointwiseMult(user->Xloc, user->Xloc, user->mult); CHKERRQ(ierr);

  // Setup CEED vectors
  ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr);
  ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr);
  CeedVectorSetArray(user->ceedvecf, CEED_MEM_HOST, CEED_USE_POINTER, x);
  CeedVectorSetArray(user->ceedvecc, CEED_MEM_HOST, CEED_USE_POINTER, y);

  // Apply CEED operator
  CeedOperatorApply(user->op, user->ceedvecf, user->ceedvecc,
                    CEED_REQUEST_IMMEDIATE);

  // Restore PETSc vectors
  CeedVectorTakeArray(user->ceedvecf, CEED_MEM_HOST, NULL);
  CeedVectorTakeArray(user->ceedvecc, CEED_MEM_HOST, NULL);
  ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
  CHKERRQ(ierr);
  ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr);

  // Local-to-global
  ierr = VecZeroEntries(Y); CHKERRQ(ierr);
  ierr = DMLocalToGlobalBegin(user->dmc, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);
  ierr = DMLocalToGlobalEnd(user->dmc, user->Yloc, ADD_VALUES, Y);
  CHKERRQ(ierr);

  PetscFunctionReturn(0);
}
#endif

// This function calculates the error in the final solution
static PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error,
                                      Vec X, CeedVector target,
                                      PetscReal *maxerror) {
  PetscErrorCode ierr;
  PetscScalar *x;
  CeedVector collocated_error;
  CeedInt length;

  PetscFunctionBeginUser;
  CeedVectorGetLength(target, &length);
  CeedVectorCreate(user->ceed, length, &collocated_error);

  // Global-to-local
  ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);

  // Setup CEED vector
  ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr);
  CeedVectorSetArray(user->xceed, CEED_MEM_HOST, CEED_USE_POINTER, x);

  // Apply CEED operator
  CeedOperatorApply(op_error, user->xceed, collocated_error,
                    CEED_REQUEST_IMMEDIATE);

  // Restore PETSc vector
  VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr);

  // Reduce max error
  *maxerror = 0;
  const CeedScalar *e;
  CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e);
  for (CeedInt i=0; i<length; i++) {
    *maxerror = PetscMax(*maxerror, PetscAbsScalar(e[i]));
  }
  CeedVectorRestoreArrayRead(collocated_error, &e);
  ierr = MPI_Allreduce(MPI_IN_PLACE, maxerror,
                       1, MPIU_REAL, MPIU_MAX, user->comm); CHKERRQ(ierr);

  // Cleanup
  CeedVectorDestroy(&collocated_error);

  PetscFunctionReturn(0);
}

#endif // setupsphere_h


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_rotg (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   float a = -1.5f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -2.12132034356f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 166)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 167)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 168)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 169)");
  };

  {
   float a = -1.5f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.80277563773f;
   float z_expected = 0.554700196225f;
   float c_expected = 0.832050294338f;
   float s_expected = 0.554700196225f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 170)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 171)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 172)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 173)");
  };

  {
   float a = -1.5f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = -1.50332963784f;
   float z_expected = 0.0665190105238f;
   float c_expected = 0.997785157857f;
   float s_expected = 0.0665190105238f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 174)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 175)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 176)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 177)");
  };

  {
   float a = -1.5f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = -1.5f;
   float z_expected = -0.0f;
   float c_expected = 1.0f;
   float s_expected = -0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 178)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 179)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 180)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 181)");
  };

  {
   float a = -1.5f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = -1.50332963784f;
   float z_expected = -0.0665190105238f;
   float c_expected = 0.997785157857f;
   float s_expected = -0.0665190105238f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 182)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 183)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 184)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 185)");
  };

  {
   float a = -1.5f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = -1.80277563773f;
   float z_expected = -0.554700196225f;
   float c_expected = 0.832050294338f;
   float s_expected = -0.554700196225f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 186)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 187)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 188)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 189)");
  };

  {
   float a = -1.5f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 2.12132034356f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 190)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 191)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 192)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 193)");
  };

  {
   float a = -1.0f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -1.80277563773f;
   float z_expected = 1.80277563773f;
   float c_expected = 0.554700196225f;
   float s_expected = 0.832050294338f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 194)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 195)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 196)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 197)");
  };

  {
   float a = -1.0f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.41421356237f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 198)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 199)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 200)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 201)");
  };

  {
   float a = -1.0f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = -1.00498756211f;
   float z_expected = 0.099503719021f;
   float c_expected = 0.99503719021f;
   float s_expected = 0.099503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 202)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 203)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 204)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 205)");
  };

  {
   float a = -1.0f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = -1.0f;
   float z_expected = -0.0f;
   float c_expected = 1.0f;
   float s_expected = -0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 206)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 207)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 208)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 209)");
  };

  {
   float a = -1.0f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = -1.00498756211f;
   float z_expected = -0.099503719021f;
   float c_expected = 0.99503719021f;
   float s_expected = -0.099503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 210)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 211)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 212)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 213)");
  };

  {
   float a = -1.0f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.41421356237f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 214)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 215)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 216)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 217)");
  };

  {
   float a = -1.0f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 1.80277563773f;
   float z_expected = -1.80277563773f;
   float c_expected = -0.554700196225f;
   float s_expected = 0.832050294338f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 218)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 219)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 220)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 221)");
  };

  {
   float a = -0.1f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -1.50332963784f;
   float z_expected = 15.0332963784f;
   float c_expected = 0.0665190105238f;
   float s_expected = 0.997785157857f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 222)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 223)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 224)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 225)");
  };

  {
   float a = -0.1f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.00498756211f;
   float z_expected = 10.0498756211f;
   float c_expected = 0.099503719021f;
   float s_expected = 0.99503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 226)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 227)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 228)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 229)");
  };

  {
   float a = -0.1f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = -0.141421356237f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 230)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 231)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 232)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 233)");
  };

  {
   float a = -0.1f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = -0.1f;
   float z_expected = -0.0f;
   float c_expected = 1.0f;
   float s_expected = -0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 234)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 235)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 236)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 237)");
  };

  {
   float a = -0.1f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = 0.141421356237f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 238)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 239)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 240)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 241)");
  };

  {
   float a = -0.1f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.00498756211f;
   float z_expected = -10.0498756211f;
   float c_expected = -0.099503719021f;
   float s_expected = 0.99503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 242)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 243)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 244)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 245)");
  };

  {
   float a = -0.1f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 1.50332963784f;
   float z_expected = -15.0332963784f;
   float c_expected = -0.0665190105238f;
   float s_expected = 0.997785157857f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 246)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 247)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 248)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 249)");
  };

  {
   float a = 0.0f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -1.5f;
   float z_expected = 1.0f;
   float c_expected = -0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 250)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 251)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 252)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 253)");
  };

  {
   float a = 0.0f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.0f;
   float z_expected = 1.0f;
   float c_expected = -0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 254)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 255)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 256)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 257)");
  };

  {
   float a = 0.0f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = -0.1f;
   float z_expected = 1.0f;
   float c_expected = -0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 258)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 259)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 260)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 261)");
  };

  {
   float a = 0.0f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = 0.0f;
   float z_expected = 0.0f;
   float c_expected = 1.0f;
   float s_expected = 0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 262)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 263)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 264)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 265)");
  };

  {
   float a = 0.0f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = 0.1f;
   float z_expected = 1.0f;
   float c_expected = 0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 266)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 267)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 268)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 269)");
  };

  {
   float a = 0.0f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.0f;
   float z_expected = 1.0f;
   float c_expected = 0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 270)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 271)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 272)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 273)");
  };

  {
   float a = 0.0f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 1.5f;
   float z_expected = 1.0f;
   float c_expected = 0.0f;
   float s_expected = 1.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 274)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 275)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 276)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 277)");
  };

  {
   float a = 0.1f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -1.50332963784f;
   float z_expected = -15.0332963784f;
   float c_expected = -0.0665190105238f;
   float s_expected = 0.997785157857f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 278)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 279)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 280)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 281)");
  };

  {
   float a = 0.1f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.00498756211f;
   float z_expected = -10.0498756211f;
   float c_expected = -0.099503719021f;
   float s_expected = 0.99503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 282)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 283)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 284)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 285)");
  };

  {
   float a = 0.1f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = -0.141421356237f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 286)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 287)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 288)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 289)");
  };

  {
   float a = 0.1f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = 0.1f;
   float z_expected = 0.0f;
   float c_expected = 1.0f;
   float s_expected = 0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 290)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 291)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 292)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 293)");
  };

  {
   float a = 0.1f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = 0.141421356237f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 294)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 295)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 296)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 297)");
  };

  {
   float a = 0.1f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.00498756211f;
   float z_expected = 10.0498756211f;
   float c_expected = 0.099503719021f;
   float s_expected = 0.99503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 298)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 299)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 300)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 301)");
  };

  {
   float a = 0.1f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 1.50332963784f;
   float z_expected = 15.0332963784f;
   float c_expected = 0.0665190105238f;
   float s_expected = 0.997785157857f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 302)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 303)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 304)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 305)");
  };

  {
   float a = 1.0f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -1.80277563773f;
   float z_expected = -1.80277563773f;
   float c_expected = -0.554700196225f;
   float s_expected = 0.832050294338f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 306)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 307)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 308)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 309)");
  };

  {
   float a = 1.0f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = -1.41421356237f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 310)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 311)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 312)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 313)");
  };

  {
   float a = 1.0f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = 1.00498756211f;
   float z_expected = -0.099503719021f;
   float c_expected = 0.99503719021f;
   float s_expected = -0.099503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 314)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 315)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 316)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 317)");
  };

  {
   float a = 1.0f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = 1.0f;
   float z_expected = 0.0f;
   float c_expected = 1.0f;
   float s_expected = 0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 318)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 319)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 320)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 321)");
  };

  {
   float a = 1.0f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = 1.00498756211f;
   float z_expected = 0.099503719021f;
   float c_expected = 0.99503719021f;
   float s_expected = 0.099503719021f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 322)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 323)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 324)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 325)");
  };

  {
   float a = 1.0f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.41421356237f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 326)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 327)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 328)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 329)");
  };

  {
   float a = 1.0f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 1.80277563773f;
   float z_expected = 1.80277563773f;
   float c_expected = 0.554700196225f;
   float s_expected = 0.832050294338f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 330)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 331)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 332)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 333)");
  };

  {
   float a = 1.5f;
   float b = -1.5f;
   float c;
   float s;
   float r_expected = -2.12132034356f;
   float z_expected = -1.41421356237f;
   float c_expected = -0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 334)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 335)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 336)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 337)");
  };

  {
   float a = 1.5f;
   float b = -1.0f;
   float c;
   float s;
   float r_expected = 1.80277563773f;
   float z_expected = -0.554700196225f;
   float c_expected = 0.832050294338f;
   float s_expected = -0.554700196225f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 338)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 339)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 340)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 341)");
  };

  {
   float a = 1.5f;
   float b = -0.1f;
   float c;
   float s;
   float r_expected = 1.50332963784f;
   float z_expected = -0.0665190105238f;
   float c_expected = 0.997785157857f;
   float s_expected = -0.0665190105238f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 342)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 343)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 344)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 345)");
  };

  {
   float a = 1.5f;
   float b = 0.0f;
   float c;
   float s;
   float r_expected = 1.5f;
   float z_expected = 0.0f;
   float c_expected = 1.0f;
   float s_expected = 0.0f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 346)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 347)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 348)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 349)");
  };

  {
   float a = 1.5f;
   float b = 0.1f;
   float c;
   float s;
   float r_expected = 1.50332963784f;
   float z_expected = 0.0665190105238f;
   float c_expected = 0.997785157857f;
   float s_expected = 0.0665190105238f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 350)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 351)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 352)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 353)");
  };

  {
   float a = 1.5f;
   float b = 1.0f;
   float c;
   float s;
   float r_expected = 1.80277563773f;
   float z_expected = 0.554700196225f;
   float c_expected = 0.832050294338f;
   float s_expected = 0.554700196225f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 354)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 355)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 356)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 357)");
  };

  {
   float a = 1.5f;
   float b = 1.5f;
   float c;
   float s;
   float r_expected = 2.12132034356f;
   float z_expected = 1.41421356237f;
   float c_expected = 0.707106781187f;
   float s_expected = 0.707106781187f;
   cblas_srotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, flteps, "srotg(case 358)");
   gsl_test_rel(b, z_expected, flteps, "srotg(case 359)");
   gsl_test_rel(c, c_expected, flteps, "srotg(case 360)");
   gsl_test_rel(s, s_expected, flteps, "srotg(case 361)");
  };

  {
   double a = -1.5;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -2.12132034356;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 362)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 363)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 364)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 365)");
  };

  {
   double a = -1.5;
   double b = -1;
   double c;
   double s;
   double r_expected = -1.80277563773;
   double z_expected = 0.554700196225;
   double c_expected = 0.832050294338;
   double s_expected = 0.554700196225;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 366)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 367)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 368)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 369)");
  };

  {
   double a = -1.5;
   double b = -0.1;
   double c;
   double s;
   double r_expected = -1.50332963784;
   double z_expected = 0.0665190105238;
   double c_expected = 0.997785157857;
   double s_expected = 0.0665190105238;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 370)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 371)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 372)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 373)");
  };

  {
   double a = -1.5;
   double b = 0;
   double c;
   double s;
   double r_expected = -1.5;
   double z_expected = -0;
   double c_expected = 1;
   double s_expected = -0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 374)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 375)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 376)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 377)");
  };

  {
   double a = -1.5;
   double b = 0.1;
   double c;
   double s;
   double r_expected = -1.50332963784;
   double z_expected = -0.0665190105238;
   double c_expected = 0.997785157857;
   double s_expected = -0.0665190105238;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 378)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 379)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 380)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 381)");
  };

  {
   double a = -1.5;
   double b = 1;
   double c;
   double s;
   double r_expected = -1.80277563773;
   double z_expected = -0.554700196225;
   double c_expected = 0.832050294338;
   double s_expected = -0.554700196225;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 382)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 383)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 384)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 385)");
  };

  {
   double a = -1.5;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 2.12132034356;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 386)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 387)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 388)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 389)");
  };

  {
   double a = -1;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -1.80277563773;
   double z_expected = 1.80277563773;
   double c_expected = 0.554700196225;
   double s_expected = 0.832050294338;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 390)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 391)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 392)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 393)");
  };

  {
   double a = -1;
   double b = -1;
   double c;
   double s;
   double r_expected = -1.41421356237;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 394)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 395)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 396)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 397)");
  };

  {
   double a = -1;
   double b = -0.1;
   double c;
   double s;
   double r_expected = -1.00498756211;
   double z_expected = 0.099503719021;
   double c_expected = 0.99503719021;
   double s_expected = 0.099503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 398)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 399)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 400)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 401)");
  };

  {
   double a = -1;
   double b = 0;
   double c;
   double s;
   double r_expected = -1;
   double z_expected = -0;
   double c_expected = 1;
   double s_expected = -0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 402)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 403)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 404)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 405)");
  };

  {
   double a = -1;
   double b = 0.1;
   double c;
   double s;
   double r_expected = -1.00498756211;
   double z_expected = -0.099503719021;
   double c_expected = 0.99503719021;
   double s_expected = -0.099503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 406)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 407)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 408)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 409)");
  };

  {
   double a = -1;
   double b = 1;
   double c;
   double s;
   double r_expected = 1.41421356237;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 410)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 411)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 412)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 413)");
  };

  {
   double a = -1;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 1.80277563773;
   double z_expected = -1.80277563773;
   double c_expected = -0.554700196225;
   double s_expected = 0.832050294338;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 414)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 415)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 416)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 417)");
  };

  {
   double a = -0.1;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -1.50332963784;
   double z_expected = 15.0332963784;
   double c_expected = 0.0665190105238;
   double s_expected = 0.997785157857;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 418)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 419)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 420)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 421)");
  };

  {
   double a = -0.1;
   double b = -1;
   double c;
   double s;
   double r_expected = -1.00498756211;
   double z_expected = 10.0498756211;
   double c_expected = 0.099503719021;
   double s_expected = 0.99503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 422)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 423)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 424)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 425)");
  };

  {
   double a = -0.1;
   double b = -0.1;
   double c;
   double s;
   double r_expected = -0.141421356237;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 426)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 427)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 428)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 429)");
  };

  {
   double a = -0.1;
   double b = 0;
   double c;
   double s;
   double r_expected = -0.1;
   double z_expected = -0;
   double c_expected = 1;
   double s_expected = -0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 430)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 431)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 432)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 433)");
  };

  {
   double a = -0.1;
   double b = 0.1;
   double c;
   double s;
   double r_expected = 0.141421356237;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 434)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 435)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 436)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 437)");
  };

  {
   double a = -0.1;
   double b = 1;
   double c;
   double s;
   double r_expected = 1.00498756211;
   double z_expected = -10.0498756211;
   double c_expected = -0.099503719021;
   double s_expected = 0.99503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 438)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 439)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 440)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 441)");
  };

  {
   double a = -0.1;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 1.50332963784;
   double z_expected = -15.0332963784;
   double c_expected = -0.0665190105238;
   double s_expected = 0.997785157857;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 442)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 443)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 444)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 445)");
  };

  {
   double a = 0;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -1.5;
   double z_expected = 1;
   double c_expected = -0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 446)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 447)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 448)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 449)");
  };

  {
   double a = 0;
   double b = -1;
   double c;
   double s;
   double r_expected = -1;
   double z_expected = 1;
   double c_expected = -0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 450)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 451)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 452)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 453)");
  };

  {
   double a = 0;
   double b = -0.1;
   double c;
   double s;
   double r_expected = -0.1;
   double z_expected = 1;
   double c_expected = -0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 454)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 455)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 456)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 457)");
  };

  {
   double a = 0;
   double b = 0;
   double c;
   double s;
   double r_expected = 0;
   double z_expected = 0;
   double c_expected = 1;
   double s_expected = 0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 458)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 459)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 460)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 461)");
  };

  {
   double a = 0;
   double b = 0.1;
   double c;
   double s;
   double r_expected = 0.1;
   double z_expected = 1;
   double c_expected = 0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 462)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 463)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 464)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 465)");
  };

  {
   double a = 0;
   double b = 1;
   double c;
   double s;
   double r_expected = 1;
   double z_expected = 1;
   double c_expected = 0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 466)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 467)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 468)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 469)");
  };

  {
   double a = 0;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 1.5;
   double z_expected = 1;
   double c_expected = 0;
   double s_expected = 1;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 470)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 471)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 472)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 473)");
  };

  {
   double a = 0.1;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -1.50332963784;
   double z_expected = -15.0332963784;
   double c_expected = -0.0665190105238;
   double s_expected = 0.997785157857;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 474)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 475)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 476)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 477)");
  };

  {
   double a = 0.1;
   double b = -1;
   double c;
   double s;
   double r_expected = -1.00498756211;
   double z_expected = -10.0498756211;
   double c_expected = -0.099503719021;
   double s_expected = 0.99503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 478)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 479)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 480)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 481)");
  };

  {
   double a = 0.1;
   double b = -0.1;
   double c;
   double s;
   double r_expected = -0.141421356237;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 482)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 483)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 484)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 485)");
  };

  {
   double a = 0.1;
   double b = 0;
   double c;
   double s;
   double r_expected = 0.1;
   double z_expected = 0;
   double c_expected = 1;
   double s_expected = 0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 486)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 487)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 488)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 489)");
  };

  {
   double a = 0.1;
   double b = 0.1;
   double c;
   double s;
   double r_expected = 0.141421356237;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 490)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 491)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 492)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 493)");
  };

  {
   double a = 0.1;
   double b = 1;
   double c;
   double s;
   double r_expected = 1.00498756211;
   double z_expected = 10.0498756211;
   double c_expected = 0.099503719021;
   double s_expected = 0.99503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 494)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 495)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 496)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 497)");
  };

  {
   double a = 0.1;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 1.50332963784;
   double z_expected = 15.0332963784;
   double c_expected = 0.0665190105238;
   double s_expected = 0.997785157857;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 498)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 499)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 500)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 501)");
  };

  {
   double a = 1;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -1.80277563773;
   double z_expected = -1.80277563773;
   double c_expected = -0.554700196225;
   double s_expected = 0.832050294338;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 502)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 503)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 504)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 505)");
  };

  {
   double a = 1;
   double b = -1;
   double c;
   double s;
   double r_expected = -1.41421356237;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 506)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 507)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 508)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 509)");
  };

  {
   double a = 1;
   double b = -0.1;
   double c;
   double s;
   double r_expected = 1.00498756211;
   double z_expected = -0.099503719021;
   double c_expected = 0.99503719021;
   double s_expected = -0.099503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 510)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 511)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 512)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 513)");
  };

  {
   double a = 1;
   double b = 0;
   double c;
   double s;
   double r_expected = 1;
   double z_expected = 0;
   double c_expected = 1;
   double s_expected = 0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 514)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 515)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 516)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 517)");
  };

  {
   double a = 1;
   double b = 0.1;
   double c;
   double s;
   double r_expected = 1.00498756211;
   double z_expected = 0.099503719021;
   double c_expected = 0.99503719021;
   double s_expected = 0.099503719021;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 518)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 519)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 520)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 521)");
  };

  {
   double a = 1;
   double b = 1;
   double c;
   double s;
   double r_expected = 1.41421356237;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 522)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 523)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 524)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 525)");
  };

  {
   double a = 1;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 1.80277563773;
   double z_expected = 1.80277563773;
   double c_expected = 0.554700196225;
   double s_expected = 0.832050294338;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 526)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 527)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 528)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 529)");
  };

  {
   double a = 1.5;
   double b = -1.5;
   double c;
   double s;
   double r_expected = -2.12132034356;
   double z_expected = -1.41421356237;
   double c_expected = -0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 530)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 531)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 532)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 533)");
  };

  {
   double a = 1.5;
   double b = -1;
   double c;
   double s;
   double r_expected = 1.80277563773;
   double z_expected = -0.554700196225;
   double c_expected = 0.832050294338;
   double s_expected = -0.554700196225;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 534)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 535)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 536)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 537)");
  };

  {
   double a = 1.5;
   double b = -0.1;
   double c;
   double s;
   double r_expected = 1.50332963784;
   double z_expected = -0.0665190105238;
   double c_expected = 0.997785157857;
   double s_expected = -0.0665190105238;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 538)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 539)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 540)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 541)");
  };

  {
   double a = 1.5;
   double b = 0;
   double c;
   double s;
   double r_expected = 1.5;
   double z_expected = 0;
   double c_expected = 1;
   double s_expected = 0;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 542)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 543)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 544)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 545)");
  };

  {
   double a = 1.5;
   double b = 0.1;
   double c;
   double s;
   double r_expected = 1.50332963784;
   double z_expected = 0.0665190105238;
   double c_expected = 0.997785157857;
   double s_expected = 0.0665190105238;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 546)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 547)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 548)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 549)");
  };

  {
   double a = 1.5;
   double b = 1;
   double c;
   double s;
   double r_expected = 1.80277563773;
   double z_expected = 0.554700196225;
   double c_expected = 0.832050294338;
   double s_expected = 0.554700196225;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 550)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 551)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 552)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 553)");
  };

  {
   double a = 1.5;
   double b = 1.5;
   double c;
   double s;
   double r_expected = 2.12132034356;
   double z_expected = 1.41421356237;
   double c_expected = 0.707106781187;
   double s_expected = 0.707106781187;
   cblas_drotg(&a, &b, &c, &s);
   gsl_test_rel(a, r_expected, dbleps, "drotg(case 554)");
   gsl_test_rel(b, z_expected, dbleps, "drotg(case 555)");
   gsl_test_rel(c, c_expected, dbleps, "drotg(case 556)");
   gsl_test_rel(s, s_expected, dbleps, "drotg(case 557)");
  };

}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/* TwoAdvSelfSims.c

Simulation calculating fixation probability of second beneficial allele, given existing
ben allele at initial frequency p. For use in the study "Linkage and the limits to natural
selection in partially selfing populations".

Simulation uses routines found with the GNU Scientific Library (GSL)
(http://www.gnu.org/software/gsl/)
Since GSL is distributed under the GNU General Public License
(http://www.gnu.org/copyleft/gpl.html), you must download it
separately from this file.

This program can be compiled with e.g. GCC using a command like:
gcc TwoAdvSelfSims -lm -lgsl -lgslcblas -I/usr/local/include -L/usr/local/lib TwoAdvSelfSims.c

Then run by executing:
./TwoAdvSelfSims N self rec ha sa hb sb p reps
Where:
- N is the population size
- self is the rate of self-fertilisation
- rec is recombination rate
- ha, hb is dominance at the original, introduced beneficial allele
- sa, sb is selection coefficient of the original, introduced beneficial allele
- p is the initial frequency of the first beneficial allele (when the second is introduced)
- reps is how many times the second allele should FIX before simulation stops
(the number of actual runs is greater due to stochastic loss of second allele)

Note that haplotypes are defined as:
x1 = ab
x2 = Ab
x3 = aB
x4 = AB

Genotypes defined as:
g11 = g1 = ab/ab
g12 = g2 = Ab/ab
g13 = g3 = aB/ab
g14 = g4 = AB/ab
g22 = g5 = Ab/Ab
g23 = g6 = Ab/aB
g24 = g7 = Ab/AB
g33 = g8 = aB/aB
g34 = g9 = aB/AB
g44 = g10 = AB/AB

Output files are the parameters;
followed by number of times each haplotype fixed;
followed by average total generations elapsed in each case;
Then total number of simulations ran;
Then fixation prob of allele, both unscaled and scaled to unlinked case,
along with 95% CI intervals for the latter case;
then number of allele fixations.

Note that 'fixation' DIFFERS depending on the inputs of sa, sb.
If sa >= sb (interference case) then 'fixation' counts as fixation of second allele on any genetic background.
If sa < sb (replacement case) then 'fixation' only considers fixation of second allele with neutral haplotype
(I.e. where the 'less fit' neutral allele at locus A fixes, instead of selected allele).

*/

/* Preprocessor statements */
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

/* Function prototypes */
void geninit(double *geninit, double FIS, const gsl_rng *r);
void selection(double *geninit);
void reproduction(double *geninit);
unsigned int hcheck(double *geninit, double *haps, unsigned int *hf, unsigned int stype);

/* Global variable declaration */
unsigned int N = 0;		/* Pop size */
double rec = 0;			/* Recombination rate */
double self = 0;		/* Rate of self-fertilisation */
double ha = 0;			/* Dominance of site A */
double sa = 0;			/* Fitness of site A */
double hb = 0;			/* Dominance of site B */
double sb = 0;			/* Fitness of site B */
double pee = 0;			/* Freq of initial sweep */

/* Main program */
int main(int argc, char *argv[]){
	unsigned int i; 				/* A counter */
	unsigned int reps;				/* Length of simulation (no. of introductions of neutral site) */
	unsigned int stype = 0;			/* What type of sim (replacement or hitch-hiking)? */
	unsigned int nsfix = 0;			/* sims where target type fixed */
	unsigned int nstot = 0;			/* total sims ran */
	unsigned int gens = 0;			/* Number gens elapsed */
	unsigned int isfin = 0;			/* Is sim finished? */
	unsigned int hf = 0;			/* The hap that fixed */
	double pf = 0;					/* Overall fix prob */
	double FIS = 0;					/* Wright's FIS */
	double StdFix = 0;				/* Standard Fixation prob if unlinked */
	double citop = 0;
	double cibot = 0;
	double nsfix2 = 0;
	double nstot2 = 0;
	char selfchar[10];
	char recchar[15];
	char hchar[10];
	char pchar[10];
	char fname[64];
	FILE *ofp_tr;

	/* GSL random number definitions */
	const gsl_rng_type * T;
	gsl_rng * r;

	/* This reads in data from command line. */
	if(argc != 10){
		fprintf(stderr,"Invalid number of input values.\n");
		exit(1);
	}
	N = strtod(argv[1],NULL);
	self = strtod(argv[2],NULL);
	rec = strtod(argv[3],NULL);
	ha = strtod(argv[4],NULL);
	sa = strtod(argv[5],NULL);
	hb = strtod(argv[6],NULL);
	sb = strtod(argv[7],NULL);
	pee = strtod(argv[8],NULL);
	reps = strtod(argv[9],NULL);

	if(sa >= sb){
		stype = 0;
	}else if(sa < sb){
		stype = 1;
	}
	FIS = self/(2.0-self);

	/* Arrays definition and memory assignment */
	double *genotype = calloc(10,sizeof(double));				/* Genotype frequencies */
	unsigned int *gensamp = calloc(10,sizeof(unsigned int));	/* New population samples */
	double *haps = calloc(4,sizeof(double));					/* Haplotypes */
	unsigned int *pfix = calloc(4,sizeof(unsigned int));		/* Haplotypes that fix */
	unsigned int *tfix = calloc(4,sizeof(unsigned int));		/* time that haps fix */

	/* create a generator chosen by the
    environment variable GSL_RNG_TYPE */

	gsl_rng_env_setup();
	if (!getenv("GSL_RNG_SEED")) gsl_rng_default_seed = time(0);
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);

	nsfix = 0;
	nstot = 0;

    while(nsfix < reps){
    	nstot++;

		/* Initialising genotypes */
		geninit(genotype,FIS,r);
		gens = 0;
		isfin = 0;
		hf = 0;

		while(isfin == 0){
    		/* Selection routine */
    		selection(genotype);

	    	/* Reproduction routine */
    		reproduction(genotype);

    		/* Sampling based on new frequencies */
	       	gsl_ran_multinomial(r,10,N,genotype,gensamp);
       		for(i = 0; i < 10; i++){
    	   		*(genotype + i) = (*(gensamp + i))/(1.0*N);
	       	}

	       	gens++;
	       	isfin = hcheck(genotype, haps, &hf, stype);
       	}

       	if(isfin == 1){
       		(*(pfix + hf))++;
    	  	(*(tfix + hf)) += gens;
	       	if(stype == 0){
       			nsfix += (*(haps + 2) + *(haps + 3));
       		}else if(stype == 1){
    	   		nsfix += (*(haps + 2));
	       	}
       	}

	}	/* End of simulation */

	pf = nsfix/(1.0*nstot);
	StdFix = 2*sb*((hb+FIS-hb*FIS)/(1+FIS));		/* Fix prob of new allele if unlinked */

	nstot2 = nstot + 3.84;
	nsfix2 = (1.0/(nstot2))*(nsfix + (3.84/2.0));
	citop = nsfix2 + 1.96*sqrt((1.0/(nstot2))*nsfix2*(1.0-nsfix2));
	cibot = nsfix2 - 1.96*sqrt((1.0/(nstot2))*nsfix2*(1.0-nsfix2));
	citop = citop/(1.0*StdFix);
	cibot = cibot/(1.0*StdFix);

	/* Printing solutions to file */
	/* First, converting values to strings */
    sprintf(selfchar, "%0.2lf",self);
    sprintf(recchar, "%0.7lf",rec);
    sprintf(hchar, "%0.3lf",ha);
    sprintf(pchar, "%0.5lf",pee);

	strcpy(fname,"sim_self");
	strcat(fname,selfchar);
	strcat(fname,"rec");
	strcat(fname,recchar);
	strcat(fname,"h");
	strcat(fname,hchar);
	strcat(fname,"p");
	strcat(fname,pchar);
	strcat(fname,".sim");

	ofp_tr = fopen(fname,"a+");
	fprintf(ofp_tr,"%d %lf %lf %lf %lf %lf %lf %lf ",N,self,rec,ha,sa,hb,sb,pee);
	for(i = 0; i < 4; i++){
		fprintf(ofp_tr,"%d ",*(pfix + i));
	}
	for(i = 0; i < 4; i++){
		fprintf(ofp_tr,"%lf ",((*(tfix + i)))/(1.0*(*(pfix + i))));
	}
	fprintf(ofp_tr,"%d %lf %lf %lf %lf %d\n",nstot,pf,pf/(1.0*StdFix),citop,cibot,reps);
	fclose(ofp_tr);

	/* Freeing memory and wrapping up */
 	gsl_rng_free(r);
 	free(tfix);
 	free(pfix);
 	free(haps);
 	free(gensamp);
	free(genotype);
	return 0;
}

/* Initialising genotypes */
void geninit(double *geninit, double FIS, const gsl_rng *r){

	unsigned int htype = 0;		/* Type of het assignment */

	double *ptype = calloc(3,sizeof(double));
	unsigned int *ctype = calloc(3,sizeof(double));

	/* Routine to determine initial background of second mutant, given first is at frequency p. */
	*(ptype + 0) = (1-pee)*(1-pee) + FIS*pee*(1-pee);
	*(ptype + 1) = 2*pee*(1-pee)*(1-FIS);
	*(ptype + 2) = pee*pee + FIS*pee*(1-pee);

	/* First initialise baseline freqs */
	*(geninit + 0) = *(ptype + 0);
	*(geninit + 1) = *(ptype + 1);
	*(geninit + 2) = 0;
	*(geninit + 3) = 0;
	*(geninit + 4) = *(ptype + 2);
	*(geninit + 5) = 0;
	*(geninit + 6) = 0;
	*(geninit + 7) = 0;
	*(geninit + 8) = 0;
	*(geninit + 9) = 0;

	/* Then decide where to add new mutant */
	gsl_ran_multinomial(r,3,1,ptype,ctype);
	if(*(ctype + 0) == 1){
		*(geninit + 0) -= 1/(1.0*N);
		*(geninit + 2) += 1/(1.0*N);
	}else if(*(ctype + 1) == 1){
		*(geninit + 1) -= 1/(1.0*N);
	  	htype = gsl_ran_bernoulli(r,0.5);
	  	if(htype == 0){
	  		*(geninit + 5) += 1/(1.0*N);
	  	}else if(htype == 1){
	  		*(geninit + 3) += 1/(1.0*N);
	  	}
	}else if(*(ctype + 2) == 1){
		*(geninit + 4) -= 1/(1.0*N);
		*(geninit + 6) += 1/(1.0*N);
	}

	free(ctype);
	free(ptype);

}	/* End of gen initiation routine */

/* Selection routine */
void selection(double *geninit){
	/* Fitness of each genotype */
	double W11, W12, W13, W14, W22, W23, W24, W33, W34, W44;
	double Wmean;				/* Mean fitness */

	W11 = 1;
	W12 = 1 + ha*sa;
	W13 = 1 + hb*sb;
	W14 = 1 + ha*sa + hb*sb;
	W22 = 1 + sa;
	W23 = 1 + ha*sa + hb*sb;
	W24 = 1 + sa + hb*sb;
	W33 = 1 + sb;
	W34 = 1 + ha*sa + sb;
	W44 = 1 + sa + sb;

	/* Mean fitness calculation */
	Wmean = ((*(geninit + 0))*W11) + ((*(geninit + 1))*W12) + ((*(geninit + 2))*W13) + ((*(geninit + 3))*W14) + ((*(geninit + 4))*W22) + ((*(geninit + 5))*W23) + ((*(geninit + 6))*W24) + ((*(geninit + 7))*W33) + ((*(geninit + 8))*W34) + ((*(geninit + 9))*W44);

	/* Changing frequencies by selection */
	*(geninit + 0) = ((*(geninit + 0))*W11)/Wmean;
	*(geninit + 1) = ((*(geninit + 1))*W12)/Wmean;
	*(geninit + 2) = ((*(geninit + 2))*W13)/Wmean;
	*(geninit + 3) = ((*(geninit + 3))*W14)/Wmean;
	*(geninit + 4) = ((*(geninit + 4))*W22)/Wmean;
	*(geninit + 5) = ((*(geninit + 5))*W23)/Wmean;
	*(geninit + 6) = ((*(geninit + 6))*W24)/Wmean;
	*(geninit + 7) = ((*(geninit + 7))*W33)/Wmean;
	*(geninit + 8) = ((*(geninit + 8))*W34)/Wmean;
	*(geninit + 9) = ((*(geninit + 9))*W44)/Wmean;

}	/* End of selection routine */

/* Reproduction routine */
void reproduction(double *geninit){
	/* Fed-in genotype frequencies (for ease of programming) */
	double g11s, g12s, g13s, g14s, g22s, g23s, g24s, g33s, g34s, g44s;
	/* Haplotypes */
	double x1, x2, x3, x4;

	/* Initial definition of genotypes */
	g11s = *(geninit + 0);
	g12s = *(geninit + 1);
	g13s = *(geninit + 2);
	g14s = *(geninit + 3);
	g22s = *(geninit + 4);
	g23s = *(geninit + 5);
	g24s = *(geninit + 6);
	g33s = *(geninit + 7);
	g34s = *(geninit + 8);
	g44s = *(geninit + 9);

	/* Baseline change in haplotype frequencies */
	x1 = g11s + (g12s + g13s + g14s)/2.0 - ((g14s - g23s)*rec)/2.0;
	x2 = g22s + (g12s + g23s + g24s)/2.0 + ((g14s - g23s)*rec)/2.0;
	x3 = g33s + (g13s + g23s + g34s)/2.0 + ((g14s - g23s)*rec)/2.0;
	x4 = g44s + (g14s + g24s + g34s)/2.0 - ((g14s - g23s)*rec)/2.0;

	/* Change in SEXUAL frequencies (both outcrossing and selfing) */
	*(geninit + 0) = (g11s + (g12s + g13s + g14s*pow((1 - rec),2) + g23s*pow(rec,2))/4.0)*self + (1 - self)*pow(x1,2);
	*(geninit + 4) = (g22s + (g12s + g24s + g23s*pow((1 - rec),2) + g14s*pow(rec,2))/4.0)*self + (1 - self)*pow(x2,2);
	*(geninit + 7) = (g33s + (g13s + g34s + g23s*pow((1 - rec),2) + g14s*pow(rec,2))/4.0)*self + (1 - self)*pow(x3,2);
	*(geninit + 9) = (g44s + (g24s + g34s + g14s*pow((1 - rec),2) + g23s*pow(rec,2))/4.0)*self + (1 - self)*pow(x4,2);
	*(geninit + 1) = ((g12s + (g14s + g23s)*(1 - rec)*rec)*self)/2.0 + 2.0*(1 - self)*x1*x2;
	*(geninit + 2) = ((g13s + (g14s + g23s)*(1 - rec)*rec)*self)/2.0 + 2.0*(1 - self)*x1*x3;
	*(geninit + 3) = ((g14s*pow((1 - rec),2) + g23s*pow(rec,2))*self)/2.0 + 2.0*(1 - self)*x1*x4;
	*(geninit + 5) = ((g23s*pow((1 - rec),2) + g14s*pow(rec,2))*self)/2.0 + 2.0*(1 - self)*x2*x3;
	*(geninit + 6) = ((g24s + (g14s + g23s)*(1 - rec)*rec)*self)/2.0 + 2.0*(1 - self)*x2*x4;
	*(geninit + 8) = ((g34s + (g14s + g23s)*(1 - rec)*rec)*self)/2.0 + 2.0*(1 - self)*x3*x4;

}	/* End of reproduction routine */

/* Has any allele fixed or not? */
unsigned int hcheck(double *geninit, double *haps, unsigned int *hf, unsigned int stype){
	/* Fed-in genotype frequencies (for ease of programming) */
	double g11s, g12s, g13s, g14s, g22s, g23s, g24s, g33s, g34s, g44s;
	unsigned int retval = 0;

	/* Initial definition of genotypes */
	g11s = *(geninit + 0);
	g12s = *(geninit + 1);
	g13s = *(geninit + 2);
	g14s = *(geninit + 3);
	g22s = *(geninit + 4);
	g23s = *(geninit + 5);
	g24s = *(geninit + 6);
	g33s = *(geninit + 7);
	g34s = *(geninit + 8);
	g44s = *(geninit + 9);

	/* Calculation of haplotypes */
	*(haps + 0) = g11s + (g12s + g13s + g14s)/2.0;
	*(haps + 1) = g22s + (g12s + g23s + g24s)/2.0;
	*(haps + 2) = g33s + (g13s + g23s + g34s)/2.0;
	*(haps + 3) = g44s + (g14s + g24s + g34s)/2.0;

/* 	printf("Haps are %lf %lf %lf %lf\n",*(haps + 0),*(haps + 1),*(haps + 2),*(haps + 3));*/

	if(*(haps + 0) == 1){
		retval = 1;
		*hf = 0;
	}
	else if(*(haps + 1) == 1){
		retval = 1;
		*hf = 1;
	}
	else if(*(haps + 2) == 1){
		retval = 1;
		*hf = 2;
	}
	else if(*(haps + 3) == 1){
		retval = 1;
		*hf = 3;
	}else if(stype == 1){
		if( (*(haps + 2) + *(haps + 3) ) == 0 ){
			retval = 2;
			*hf = 0;
		}
	}

	return retval;

}	/* End of hap check routine */

/* End of program */

max_stars_repo_name: MattHartfield/TwoAdvSelfSims

max_issues_repo_name: MattHartfield/TwoAdvSelfSims

max_forks_repo_name: MattHartfield/TwoAdvSelfSims


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_chebyshev.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_randist.h>

double funcExp(double x, void *p)
{
    (void)(p); //avoid unused parameter warning
    return exp(pow(x, 2));
}

double funcAbs(double x, void *p)
{
    (void)(p); //avoid unused parameter warning
    return fabs(x + pow(x, 3));
}

double funcSign(double x, void *p)
{
    (void)(p); //avoid unused parameter warning
    if (x == 0) return 0.0;
    else if (x < 0) return -1.0;

    return 1.0;
}

void calculateExp()
{
    const int stepsCount = 1000;
    FILE* outputR1 = fopen("cheb_exp_r1.txt", "w");
    FILE* outputR2 = fopen("cheb_exp_r2.txt", "w");
    FILE* outputR4 = fopen("cheb_exp_r4.txt", "w");
    FILE* outputR8 = fopen("cheb_exp_r8.txt", "w");
    FILE* outputR16 = fopen("cheb_exp_r16.txt", "w");
    FILE* outputR25 = fopen("cheb_exp_r25.txt", "w");
    FILE* outputR40 = fopen("cheb_exp_r40.txt", "w");

    gsl_cheb_series* cs = gsl_cheb_alloc(40);

    gsl_function function;

    function.function = funcExp;
    function.params = 0;

    gsl_cheb_init(cs, &function, -1.0, 1.0);
    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double r1 = GSL_FN_EVAL(&function, x);
        double r2 = gsl_cheb_eval_n(cs, 2, x);
        double r4 = gsl_cheb_eval_n(cs, 4, x);
        double r8 = gsl_cheb_eval_n(cs, 8, x);
        double r16 = gsl_cheb_eval_n(cs, 16, x);
        double r25 = gsl_cheb_eval_n(cs, 25, x);
        double r40 = gsl_cheb_eval_n(cs, 40, x);

        fprintf (outputR1,"%g %g\n", x, r1);
        fprintf (outputR2,"%g %g\n", x, r2);
        fprintf (outputR4,"%g %g\n", x, r4);
        fprintf (outputR8,"%g %g\n", x, r8);
        fprintf (outputR16,"%g %g\n", x, r16);
        fprintf (outputR25,"%g %g\n", x, r25);
        fprintf (outputR40,"%g %g\n", x, r40);
    }

    gsl_cheb_free(cs);
}

void calculateAbs()
{
    const int stepsCount = 1000;
    FILE* outputR1 = fopen("cheb_abs_r1.txt", "w");
    FILE* outputR2 = fopen("cheb_abs_r2.txt", "w");
    FILE* outputR4 = fopen("cheb_abs_r4.txt", "w");
    FILE* outputR8 = fopen("cheb_abs_r8.txt", "w");
    FILE* outputR16 = fopen("cheb_abs_r16.txt", "w");
    FILE* outputR25 = fopen("cheb_abs_r25.txt", "w");
    FILE* outputR40 = fopen("cheb_abs_r40.txt", "w");

    gsl_cheb_series* cs = gsl_cheb_alloc(40);

    gsl_function function;

    function.function = funcAbs;
    function.params = 0;

    gsl_cheb_init(cs, &function, -1.0, 1.0);
    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double r1 = GSL_FN_EVAL(&function, x);
        double r2 = gsl_cheb_eval_n(cs, 2, x);
        double r4 = gsl_cheb_eval_n(cs, 4, x);
        double r8 = gsl_cheb_eval_n(cs, 8, x);
        double r16 = gsl_cheb_eval_n(cs, 16, x);
        double r25 = gsl_cheb_eval_n(cs, 25, x);
        double r40 = gsl_cheb_eval_n(cs, 40, x);

        fprintf (outputR1,"%g %g\n", x, r1);
        fprintf (outputR2,"%g %g\n", x, r2);
        fprintf (outputR4,"%g %g\n", x, r4);
        fprintf (outputR8,"%g %g\n", x, r8);
        fprintf (outputR16,"%g %g\n", x, r16);
        fprintf (outputR25,"%g %g\n", x, r25);
        fprintf (outputR40,"%g %g\n", x, r40);
    }

    gsl_cheb_free(cs);
}

void calculateSign()
{
    const int stepsCount = 1000;
    FILE* outputR1 = fopen("cheb_sign_r1.txt", "w");
    FILE* outputR2 = fopen("cheb_sign_r2.txt", "w");
    FILE* outputR4 = fopen("cheb_sign_r4.txt", "w");
    FILE* outputR8 = fopen("cheb_sign_r8.txt", "w");
    FILE* outputR16 = fopen("cheb_sign_r16.txt", "w");
    FILE* outputR25 = fopen("cheb_sign_r25.txt", "w");
    FILE* outputR40 = fopen("cheb_sign_r40.txt", "w");

    gsl_cheb_series* cs = gsl_cheb_alloc(40);

    gsl_function function;

    function.function = funcSign;
    function.params = 0;

    gsl_cheb_init(cs, &function, -1.0, 1.0);
    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double r1 = GSL_FN_EVAL(&function, x);
        double r2 = gsl_cheb_eval_n(cs, 2, x);
        double r4 = gsl_cheb_eval_n(cs, 4, x);
        double r8 = gsl_cheb_eval_n(cs, 8, x);
        double r16 = gsl_cheb_eval_n(cs, 16, x);
        double r25 = gsl_cheb_eval_n(cs, 25, x);
        double r40 = gsl_cheb_eval_n(cs, 40, x);

        fprintf (outputR1,"%g %g\n", x, r1);
        fprintf (outputR2,"%g %g\n", x, r2);
        fprintf (outputR4,"%g %g\n", x, r4);
        fprintf (outputR8,"%g %g\n", x, r8);
        fprintf (outputR16,"%g %g\n", x, r16);
        fprintf (outputR25,"%g %g\n", x, r25);
        fprintf (outputR40,"%g %g\n", x, r40);
    }

    gsl_cheb_free(cs);
}

void fitLinear()
{
    const int stepsCount = 1000;
    double xArray[stepsCount * 2];
    double yArray[stepsCount * 2];

    FILE* outputOriginal = fopen("linear_original.txt", "w");
    FILE* outputFit = fopen("linear_fit.txt", "w");
    FILE* outputCheb = fopen("linear_cheb.txt", "w");

    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double y = pow(0.5, pow(x, 2) + 2 * x);

        xArray[i + stepsCount] = x;
        yArray[i + stepsCount] = y;

        fprintf (outputOriginal,"%g %g\n", x, y);
    }

    double c0 = 0;
    double c1 = 0;
    double cov00 = 0;
    double cov01 = 0;
    double cov11 = 0;
    double sumsq = 0;
    gsl_fit_linear(xArray, 1, yArray, 1, stepsCount * 2, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

//------
    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double y = c0 + c1 * x;

        xArray[i + stepsCount] = x;
        yArray[i + stepsCount] = y;

        fprintf (outputFit,"%g %g\n", x, y);
    }

    printf("c0: %f --- c1: %f\n", c0, c1);
}

double funcCheb(double x, void *p)
{
    (void)(p); //avoid unused parameter warning
    return pow(0.5, pow(x, 2) + 2 * x);
}

void fitCheb()
{
    const int stepsCount = 1000;
    FILE* outputR4 = fopen("cheb_fit_r4.txt", "w");
    FILE* outputR40 = fopen("cheb_fit_r40.txt", "w");

    gsl_cheb_series* cs = gsl_cheb_alloc(40);
    gsl_function function;
    function.function = funcCheb;
    function.params = 0;

    gsl_cheb_init(cs, &function, -1.0, 1.0);
    for (int i = -stepsCount; i <= stepsCount; i++)
    {
        double x = (double)i / (double)stepsCount;
        double r4 = gsl_cheb_eval_n(cs, 4, x);
        double r40 = gsl_cheb_eval_n(cs, 40, x);
        fprintf (outputR4,"%g %g\n", x, r4);
        fprintf (outputR40,"%g %g\n", x, r40);
    }
    gsl_cheb_free(cs);
}

int main (int argc, char* argv[])
{
    calculateExp();
    calculateSign();
    calculateAbs();
    fitLinear();
    fitCheb();
    return 0;
}

/*

plot "cheb_exp_r1.txt" with lines, \
    "cheb_exp_r2.txt" with lines, \
    "cheb_exp_r4.txt" with lines, \
    "cheb_exp_r8.txt" with lines, \
    "cheb_exp_r16.txt" with lines, \
    "cheb_exp_r25.txt" with lines, \
    "cheb_exp_r40.txt" with lines

plot "cheb_abs_r1.txt" with lines, \
    "cheb_abs_r2.txt" with lines, \
    "cheb_abs_r4.txt" with lines, \
    "cheb_abs_r8.txt" with lines, \
    "cheb_abs_r16.txt" with lines, \
    "cheb_abs_r25.txt" with lines, \
    "cheb_abs_r40.txt" with lines

plot "cheb_sign_r1.txt" with lines, \
    "cheb_sign_r2.txt" with lines, \
    "cheb_sign_r4.txt" with lines, \
    "cheb_sign_r8.txt" with lines, \
    "cheb_sign_r16.txt" with lines, \
    "cheb_sign_r25.txt" with lines, \
    "cheb_sign_r40.txt" with lines

================

plot "linear_original.txt" with lines title 'Original', \
    "linear_fit.txt" with lines title 'Linear least sqr fit', \
    "cheb_fit_r4.txt" with lines title 'Chebyshev approx R4'

*/