text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "ellipsoid/ellipsoid.h"
#include "ellipsoid/ellSolv.h"
#include "sphere/sphere.h"
#include "sphere/sphSolv.h"
#include "constants.h"
#include "testProblem.h"

//compares ellipsoidal and spherical
//has ellipsoidal distribution of charges just inside ellipsoidal boundary
//computes the solution at two point, one inside smallest enclosing sphere where spherical diverges
//other solution point is outside of sphere and ellipse to compare convergence
#undef __FUNCT__
#define __FUNCT__ "EllVsSphConvergence"
PetscErrorCode EllVsSphConvergence()
{
  PetscErrorCode ierr;
  //semi-axes of ellipsoid
  PetscReal a = 3.0;
  PetscReal b = 2.5;
  PetscReal c = 1.2;
  //sphere radius
  PetscReal sphRad = 3;
  //dielecric permitivities both set to 1
  const PetscReal eps1 = 1.0;
  const PetscReal eps2 = 1.0;
  //charges placed on ellipsoid with sem-axes smaller by epsilon
  PetscReal epsilon = .5;
  //interior charge XYZ and magnitude
  const PetscInt SQRT_NUM_CHARGES = 3;
  Vec chargeXYZ;
  Vec chargeMag;
  //maximum expansion order for test
  const PetscInt MAX_N = 20;
  //solutions vector
  const PetscInt NUM_SOL_PTS = 2;
  Vec ellSolutions[MAX_N], sphSolutions[MAX_N];
  Vec ellError[MAX_N], sphError[MAX_N];
  Vec exactSolution;
  Vec solXYZ;
  //solution point inside of sphere
  PetscReal sol1X = .1;
  PetscReal sol1Y = .1;
  PetscReal sol1Z = 1.5;
  //solution point outside of sphere
  PetscReal sol2X = 2.1;
  PetscReal sol2Y = 2.6;
  PetscReal sol2Z = 1;
  //pointer for accessing petsc vectors
  PetscScalar *vecPt;
  PetscFunctionBegin;

  //create, set charge magnitudes to 1
  ierr = VecCreateSeq(PETSC_COMM_SELF, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, &chargeMag);CHKERRQ(ierr);
  ierr = VecSet(chargeMag, 1.0);CHKERRQ(ierr);

  //create, set charge XYZ values to be on a-eps,b-eps,c-eps ellipsoid
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, &chargeXYZ);CHKERRQ(ierr);
  ierr = VecGetArray(chargeXYZ, &vecPt);CHKERRQ(ierr);
  PetscReal pu, pv;
  for(PetscInt i=0; i<SQRT_NUM_CHARGES; ++i) {
    for(PetscInt j=0; j<SQRT_NUM_CHARGES; ++j) {
      pu = i*2.0*PETSC_PI/(SQRT_NUM_CHARGES-1);
      pv = j*PETSC_PI/(SQRT_NUM_CHARGES-1);
      PetscInt ind = i*SQRT_NUM_CHARGES + j;
      //there is a slight issue with on-axes coordinate transforms which is what the +.01 is for
      vecPt[3*ind + 0] = (a-epsilon)*PetscCosReal(pu)*PetscSinReal(pv) +.01;
      vecPt[3*ind + 1] = (b-epsilon)*PetscSinReal(pu)*PetscSinReal(pv) + .01;
      vecPt[3*ind + 2] = (c-epsilon)*PetscCosReal(pv) + .01;
      printf("pu: %3.3f pv: %3.3f\n", pu, pv);
    }
  }
  ierr = VecRestoreArray(chargeXYZ, &vecPt);CHKERRQ(ierr);

  //set solution points vector
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*NUM_SOL_PTS, &solXYZ);CHKERRQ(ierr);
  ierr = VecGetArray(solXYZ, &vecPt);CHKERRQ(ierr);
  vecPt[0] = sol1X; vecPt[1] = sol1Y; vecPt[2] = sol1Z;
  vecPt[3] = sol2X; vecPt[4] = sol2Y; vecPt[5] = sol2Z;
  ierr = VecRestoreArray(solXYZ, &vecPt);CHKERRQ(ierr);

  //caulculate exact solution
  ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, &exactSolution);CHKERRQ(ierr);
  ierr = CoulombExact(eps1, chargeXYZ, chargeMag, solXYZ, exactSolution);CHKERRQ(ierr);

  //initialize solutions and error vectors
  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, ellSolutions+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, sphSolutions+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, ellError+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, sphError+k);CHKERRQ(ierr);
  }

  //compute solutions for expansions up to order MAX_N
  for(PetscInt k=1; k<=MAX_N; ++k) {
    printf("N = %d\n", k);
    printf("calculating ellipsoidal solution...\n");
    ierr = CalcEllipsoidCoulombPotential(a, b, c, eps1, eps2, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, chargeXYZ, chargeMag, NUM_SOL_PTS, solXYZ, k, ellSolutions[k-1]);CHKERRQ(ierr);
    printf("calculating spherical solution...\n");
    ierr = CalcSphericalCoulombPotential(1.0, eps1, eps2, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, chargeXYZ, chargeMag, NUM_SOL_PTS, solXYZ, k, sphSolutions[k-1]);CHKERRQ(ierr);
  }

  //compute errors
  printf("Computing errors...\n");
  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecCopy(exactSolution, ellError[k]);CHKERRQ(ierr);
    ierr = VecCopy(exactSolution, sphError[k]);CHKERRQ(ierr);
    ierr = VecAXPY(ellError[k], -1.0, ellSolutions[k]);CHKERRQ(ierr);
    ierr = VecAXPY(sphError[k], -1.0, sphSolutions[k]);CHKERRQ(ierr);
    ierr = VecAbs(ellError[k]);CHKERRQ(ierr);
    ierr = VecAbs(sphError[k]);CHKERRQ(ierr);
  }

  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecGetArrayRead(ellError[k], &vecPt);
    printf("Error[%d] = %4.4e\n", k, vecPt[0]);
    ierr = VecRestoreArrayRead(ellError[k], &vecPt);
  }
  printf("Exact: \n");
  ierr = VecView(exactSolution, PETSC_VIEWER_STDOUT_SELF);
  printf("Ell: \n");
  ierr = VecView(ellSolutions[MAX_N-1], PETSC_VIEWER_STDOUT_SELF);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "runTest1"
void runTest1() {

  const int nCharges = 1;

  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0, 64);

  //grid dimensions
  int nx = 10;
  int ny = 10;
  int nz = 2;
  int nPoints = nx*ny*nz;
  double xa = -5;
  double xb = 5;
  double ya = -5;
  double yb = 5;
  double za = -3;
  double zb = 3;
  double xh = (xb - xa)/(nx-1);
  double yh = (yb - ya)/(ny-1);
  double zh = (zb - za)/(nz-1);

  double charges[nCharges];
  for(int i=0; i<nCharges; ++i)
    charges[i] = 1.0;

  //source charges
  SPoint sPoints[nCharges];
  Point  ePoints[nCharges];
  sPoints[0].x1 = 2.02;
  sPoints[0].x2 = -0.4;
  sPoints[0].x3 = -.463;
  ePoints[0].x1 = sPoints[0].x1;
  ePoints[0].x2 = sPoints[0].x2;
  ePoints[0].x3 = sPoints[0].x3;

  //permitivity constants
  double e1 = 1.0;
  double e2 = 1.0;

  printf("entering\n");
  //solution vector
  SPoint *rs = (SPoint*) malloc(sizeof(SPoint)*nPoints);
  Point  *re = (Point*)  malloc(sizeof(Point) *nPoints);
  SPoint *ro = (SPoint*) malloc(sizeof(SPoint)*nPoints);
  double *exact = (double*) malloc(sizeof(double)*nPoints);
  for(int i=0; i<nx; ++i) {
    for(int j=0; j<ny; ++j) {
      for(int k=0; k<nz; ++k) {
	int index = k*ny*nx + j*nx + i;
	rs[index].x1 = xa + i*xh;
	rs[index].x2 = ya + j*yh;
	rs[index].x3 = za + k*zh;
	printf("rs[3] = %15.15f\n", rs[index].x3);
	re[index].x1 = rs[index].x1;
	re[index].x2 = rs[index].x2;
	re[index].x3 = rs[index].x3;
	ro[index].x1 = re[index].x1;
	ro[index].x2 = re[index].x2;
	ro[index].x3 = re[index].x3;

	exact[index] = 0.0;
	for(int num=0; num<nCharges; ++num)
	  exact[index] += charges[num]/(e1*sqrt((sPoints[num].x1 - ro[index].x1)*(sPoints[num].x1 - ro[index].x1) +
					    (sPoints[num].x2 - ro[index].x2)*(sPoints[num].x2 - ro[index].x2) +
					     (sPoints[num].x3 - ro[index].x3)*(sPoints[num].x3 - ro[index].x3)));
      }
    }
  }
  printf("exiting\n");

  //convert to spherical/ellipsoidal coordinates
  for(int i=0; i<nCharges; ++i) {
    cartesianToSpherical(sPoints+i);
    cartesianToEllipsoidal(&e, ePoints+i);
  }

  //convert to spherical/ellipsoidal coordinates
  for(int i=0; i<nPoints; ++i) {
    cartesianToSpherical(rs + i);
    cartesianToEllipsoidal(&e, re + i);
  }

  //set up problem
  SProblem sProblem = {.positions = sPoints, .charges = charges, .nCharges = nCharges, .e1 = e1, .e2 = e2, .b = 3.0};
  Problem  eProblem = {.positions = ePoints, .charges = charges, .nCharges = nCharges, .e1 = e1, .e2 = e2, .e = &e};

  double *ssolution = (double*) malloc(sizeof(double)*nPoints);
  double *esolution = (double*) malloc(sizeof(double)*nPoints);
  calcCoulombSphericalGrid(&sProblem, 14, rs, nPoints, ssolution);
  calcCoulombEllipsoidalGrid(&eProblem, 14, re, nPoints, esolution);

  //add coulomb inside for plotting. Add later to calculation
  for(int i=0; i<nPoints; ++i) {
    for(int k=0; k<nCharges; ++k) {
  //if inside sphere
  //
      //if(
    }
  }

  printf("exact[0] = %15.15f\n", exact[0]);
  //compare errors
  double *sError = (double*) malloc(sizeof(double)*nPoints);
  double *eError = (double*) malloc(sizeof(double)*nPoints);
  for(int i=0; i<nx; ++i) {
    for(int j=0; j<ny; ++j) {
      for(int k=0; k<nz; ++k) {
	int index = k*ny*nx + j*nx + i;
	sError[index] = fabs(ssolution[index] - exact[index]);
	eError[index] = fabs(esolution[index] - exact[index]);
	if(i==1)
	  printf("%d - s: %4.4e e: %4.4e\n", i, sError[index], eError[index]);
      }
    }
  }

}

#undef __FUNCT__
#define __FUNCT__ "RunArg"
PetscErrorCode RunArg()
{
  PetscErrorCode ierr;
  char buff[64];
  int npts;
  double a, b, c;
  double *xyz;
  double *chargeValues;
  double *solution;
  double freeEnergy;

  PetscFunctionBegin;

  //constants
  const double q     = ELECTRON_CHARGE;
  const double Na    = AVOGADRO_NUMBER;
  const double JperC = 4.184; /* Jouled/Calorie */
  const double cf    = Na * (q*q/EPSILON_0)/JperC * (1e10/1000) * 1/4/PETSC_PI; /* kcal ang/mol */

  FILE *fp = fopen("ellArgData.txt", "r");
  // READ NUMBER OF POINTS
  fscanf(fp, "%s", buff); //skip "size: "
  fscanf(fp, "%d", &npts);

  //initialize local charges/coordinates arrays
  xyz     = (double*) malloc(sizeof(double)*3*npts);
  chargeValues = (double*) malloc(sizeof(double)*npts);
  // READ A,B,C ELLIPSOID VALUES
  fscanf(fp, "%s", buff); //skip "a: "
  fscanf(fp, "%lf", &a); //read a value
  fscanf(fp, "%s", buff); //skip "b: "
  fscanf(fp, "%lf", &b); //read b value
  fscanf(fp, "%s", buff); //skip "c: "
  fscanf(fp, "%lf", &c);
  // READ CHARGE X,Y,Z AND CHARGE DATA
  for(int i=0; i<npts; ++i) {
    fscanf(fp, "%lf", xyz+3*i+0); //read x coordinate
    fscanf(fp, "%lf", xyz+3*i+1); //read y coordinate
    fscanf(fp, "%lf", xyz+3*i+2); //read z coordinate
    fscanf(fp, "%lf", chargeValues+i); //read charge
    printf("%lf %lf %lf %lf\n", xyz[3*i+0], xyz[3*i+1], xyz[3*i+2], chargeValues[i]);
  }
  fclose(fp);

  //initialize ellipsoidal system
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, a, b, c, 64);

  //set charge positions and calculation points to be the same
  Point *solPoints = (Point*) malloc(sizeof(Point)*npts);
  Point *chargePoints = (Point*) malloc(sizeof(Point)*npts);
  for(int i=0; i<npts; ++i) {
    solPoints[i].x1 = xyz[3*i+0];
    solPoints[i].x2 = xyz[3*i+1];
    solPoints[i].x3 = xyz[3*i+2];
    solPoints[i].type = 'c';
    chargePoints[i] = solPoints[i];
  }

  //convert to ellipsoidal coordinates
  for(int i=0; i<npts; ++i) {
    cartesianToEllipsoidal(&e, solPoints+i);
    cartesianToEllipsoidal(&e, chargePoints+i);
  }

  //initialize ellipsoidal problem context
  Problem prob;
  prob.e = &e;
  prob.positions = chargePoints;
  prob.charges = chargeValues;
  prob.nCharges = npts;
  prob.e1 = 4.0;
  prob.e2 = 80.0;

  //solution vector
  solution = (double*) malloc(sizeof(double)*npts);

  //calculate potential on solPoints
  calcCoulombEllipsoidalGrid(&prob, 10, solPoints, npts, solution);
  printf("wow!\n");
  for(int i=0; i<npts; ++i)
    printf("solution[%d] = %15.15f\n", i, solution[i]);
  for(int i=0; i<npts; ++i)
    printf("chargeValues[%d] = %15.15f\n", i, chargeValues[i]);

  //calculate free energy
  freeEnergy = 0;
  for(int i=0; i<npts; ++i) {
    freeEnergy += solution[i]*chargeValues[i];
    printf("freeEnergy[%d] = %15.15f\n", i, freeEnergy*.5*cf);
  }

  printf("\n");
  printf("Free energy: %15.15f\n", freeEnergy*.5*cf);

  free(xyz);
  free(chargeValues);
  free(solPoints);
  free(solution);
  free(chargePoints);
  PetscFunctionReturn(0);
}

/* DOESN'T DO ANYTHING YET */
#undef __FUNCT__
#define __FUNCT__ "RunArgWorkPrec"
PetscErrorCode RunArgWorkPrec()
{
  PetscErrorCode ierr;
  char buff[64];
  int npts;
  double a, b, c;
  double *xyz;
  double *chargeValues;
  double *solution;
  double freeEnergy;

  PetscFunctionBegin;

  //constants
  const double q     = ELECTRON_CHARGE;
  const double Na    = AVOGADRO_NUMBER;
  const double JperC = 4.184; /* Jouled/Calorie */
  const double cf    = Na * (q*q/EPSILON_0)/JperC * (1e10/1000) * 1/4/PETSC_PI; /* kcal ang/mol */

  FILE *fp = fopen("ellArgData.txt", "r");
  // READ NUMBER OF POINTS
  fscanf(fp, "%s", buff); //skip "size: "
  fscanf(fp, "%d", &npts);

  //initialize local charges/coordinates arrays
  xyz     = (double*) malloc(sizeof(double)*3*npts);
  chargeValues = (double*) malloc(sizeof(double)*npts);
  // READ A,B,C ELLIPSOID VALUES
  fscanf(fp, "%s", buff); //skip "a: "
  fscanf(fp, "%lf", &a); //read a value
  fscanf(fp, "%s", buff); //skip "b: "
  fscanf(fp, "%lf", &b); //read b value
  fscanf(fp, "%s", buff); //skip "c: "
  fscanf(fp, "%lf", &c);
  // READ CHARGE X,Y,Z AND CHARGE DATA
  for(int i=0; i<npts; ++i) {
    fscanf(fp, "%lf", xyz+3*i+0); //read x coordinate
    fscanf(fp, "%lf", xyz+3*i+1); //read y coordinate
    fscanf(fp, "%lf", xyz+3*i+2); //read z coordinate
    fscanf(fp, "%lf", chargeValues+i); //read charge
    printf("%lf %lf %lf %lf\n", xyz[3*i+0], xyz[3*i+1], xyz[3*i+2], chargeValues[i]);
  }
  fclose(fp);

  //initialize ellipsoidal system
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, a, b, c, 64);

  //set charge positions and calculation points to be the same
  Point *solPoints = (Point*) malloc(sizeof(Point)*npts);
  Point *chargePoints = (Point*) malloc(sizeof(Point)*npts);
  for(int i=0; i<npts; ++i) {
    solPoints[i].x1 = xyz[3*i+0];
    solPoints[i].x2 = xyz[3*i+1];
    solPoints[i].x3 = xyz[3*i+2];
    solPoints[i].type = 'c';
    chargePoints[i] = solPoints[i];
  }

  //convert to ellipsoidal coordinates
  for(int i=0; i<npts; ++i) {
    cartesianToEllipsoidal(&e, solPoints+i);
    cartesianToEllipsoidal(&e, chargePoints+i);
  }

  //initialize ellipsoidal problem context
  Problem prob;
  prob.e = &e;
  prob.positions = chargePoints;
  prob.charges = chargeValues;
  prob.nCharges = npts;
  prob.e1 = 4.0;
  prob.e2 = 80.0;

  //solution vector
  solution = (double*) malloc(sizeof(double)*npts);

  //calculate potential on solPoints
  calcCoulombEllipsoidalGrid(&prob, 10, solPoints, npts, solution);
  printf("wow!\n");
  for(int i=0; i<npts; ++i)
    printf("solution[%d] = %15.15f\n", i, solution[i]);
  for(int i=0; i<npts; ++i)
    printf("chargeValues[%d] = %15.15f\n", i, chargeValues[i]);

  //calculate free energy
  freeEnergy = 0;
  for(int i=0; i<npts; ++i) {
    freeEnergy += solution[i]*chargeValues[i];
    printf("freeEnergy[%d] = %15.15f\n", i, freeEnergy*.5*cf);
  }

  free(xyz);
  free(chargeValues);
  free(solPoints);
  free(solution);
  free(chargePoints);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "RunArgTester"
PetscErrorCode RunArgTester()
{
  PetscErrorCode ierr;
  char buff[64];
  int npts;
  double a, b, c;
  double *xyz;
  double *chargeValues;
  double *solution;
  double freeEnergy;

  PetscFunctionBegin;

  //constants
  const double q     = ELECTRON_CHARGE;
  const double Na    = AVOGADRO_NUMBER;
  const double JperC = 4.184; /* Jouled/Calorie */
  const double cf    = Na * (q*q/EPSILON_0)/JperC * (1e10/1000) * 1/4/PETSC_PI; /* kcal ang/mol */

  FILE *fp = fopen("../../pointbem/src/ellData.txt", "r");
  // READ NUMBER OF POINTS
  fscanf(fp, "%s", buff); //skip "size: "
  fscanf(fp, "%d", &npts);

  //initialize local charges/coordinates arrays
  xyz     = (double*) malloc(sizeof(double)*3*npts);
  chargeValues = (double*) malloc(sizeof(double)*npts);
  // READ A,B,C ELLIPSOID VALUES
  fscanf(fp, "%s", buff); //skip "a: "
  fscanf(fp, "%lf", &a); //read a value
  fscanf(fp, "%s", buff); //skip "b: "
  fscanf(fp, "%lf", &b); //read b value
  fscanf(fp, "%s", buff); //skip "c: "
  fscanf(fp, "%lf", &c);
  // READ CHARGE X,Y,Z AND CHARGE DATA
  for(int i=0; i<npts; ++i) {
    fscanf(fp, "%lf", xyz+3*i+0); //read x coordinate
    fscanf(fp, "%lf", xyz+3*i+1); //read y coordinate
    fscanf(fp, "%lf", xyz+3*i+2); //read z coordinate
    fscanf(fp, "%lf", chargeValues+i); //read charge
    printf("%lf %lf %lf %lf\n", xyz[3*i+0], xyz[3*i+1], xyz[3*i+2], chargeValues[i]);
  }
  fclose(fp);

  //initialize ellipsoidal system
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, a, b, c, 64);

  //set charge positions and calculation points to be the same
  Point *solPoints = (Point*) malloc(sizeof(Point)*npts);
  Point *chargePoints = (Point*) malloc(sizeof(Point)*npts);
  for(int i=0; i<npts; ++i) {
    solPoints[i].x1 = xyz[3*i+0];
    solPoints[i].x2 = xyz[3*i+1];
    solPoints[i].x3 = xyz[3*i+2];
    solPoints[i].type = 'c';
    chargePoints[i] = solPoints[i];
  }

  Vec valsbef; //testing
  PetscScalar* valsbefArray; //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, npts, &valsbef); CHKERRQ(ierr); //testing
  ierr = VecGetArray(valsbef, &valsbefArray); CHKERRQ(ierr); //testing
  //convert to ellipsoidal coordinates
  for(int i=0; i<npts; ++i) {
    cartesianToEllipsoidal(&e, solPoints+i);
    ierr = CalcSolidInteriorHarmonic(&e, solPoints[i].x1, solPoints[i].x2, solPoints[i].x3, 2, 3, valsbefArray+i); //testing
    cartesianToEllipsoidal(&e, chargePoints+i);
  }
  ierr = VecRestoreArray(valsbef, &valsbefArray); CHKERRQ(ierr); //testing
  Vec ellPoints, vals; //testing
  PetscScalar* ellPointsArray; //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*npts, &ellPoints); CHKERRQ(ierr); //testing
  ierr = VecGetArray(ellPoints, &ellPointsArray); CHKERRQ(ierr); //testing
  for(PetscInt i=0; i<npts; ++i) { //testing
    ellPointsArray[3*i+0] = solPoints[i].x1; //testing
    ellPointsArray[3*i+1] = solPoints[i].x2; //testing
    ellPointsArray[3*i+2] = solPoints[i].x3; //testing
  }
  ierr = VecRestoreArray(ellPoints, &ellPointsArray); CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, npts, &vals); CHKERRQ(ierr); //testing
  ierr = CalcSolidInteriorHarmonicVec(&e, ellPoints, 2, 3, vals); //testing
  PetscReal diff, val1, val2; //testing
  diff = 0;
  for(PetscInt i=0; i<npts; ++i) { //testing
    ierr = VecGetValues(vals, 1, &i, &val1); CHKERRQ(ierr); //testing
    ierr = VecGetValues(valsbef, 1, &i, &val2); CHKERRQ(ierr); //testing
    diff += PetscSqrtReal(PetscAbsReal(val1 - val2)); //testing
  }
  printf("Solid Harmonic Error: %3.3e\n", diff); //testing

  //initialize ellipsoidal problem context
  Problem prob;
  prob.e = &e;
  prob.positions = chargePoints;
  prob.charges = chargeValues;
  prob.nCharges = npts;
  prob.e1 = 4.0;
  prob.e2 = 80.0;

  /* ---------------------------- */
  /* NOW WE TEST THE Gnp FUNCTION */
  /* ---------------------------- */
  PetscInt NMAX = 5;
  PetscInt count = 0;

  for(PetscInt n=0; n<=NMAX; ++n) {
    for(PetscInt p=0; p<2*n+1; ++p)
      count++;
  }
  Vec srcCharges, GnpVals, GnpValsOld, BnpVals, BnpValsOld, CnpVals, CnpValsOld;
  Vec errorV;
  ierr = VecCreateSeq(PETSC_COMM_SELF, npts, &srcCharges);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &GnpValsOld);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &BnpVals);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &BnpValsOld);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &CnpVals);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &CnpValsOld);CHKERRQ(ierr); //testing
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, &errorV);CHKERRQ(ierr); //testing
  for(PetscInt i=0; i<npts; ++i)
    ierr = VecSetValues(srcCharges, 1, &i, chargeValues+i, INSERT_VALUES);CHKERRQ(ierr); //testing
  ierr = VecAssemblyBegin(srcCharges);CHKERRQ(ierr); ierr = VecAssemblyEnd(srcCharges);CHKERRQ(ierr); //testing

  //WOULD ACTUALLY USE CHARGE POINTS, NOT TARGET POINTS
  //Calculate values with new function
  ierr = CalcCoulombEllCoefs(&e, npts, ellPoints, srcCharges, NMAX, &GnpVals);CHKERRQ(ierr); //testing
  ierr = CalcReactAndExtCoefsFromCoulomb(&e, prob.e1, prob.e2, NMAX, GnpVals, BnpVals, CnpVals);CHKERRQ(ierr);

  //Calculate values with old function
  PetscReal GnpOld, BnpOld, CnpOld;
  count = 0;
  for(PetscInt n=0; n<=NMAX; ++n) {
    for(PetscInt p=0; p<2*n+1; ++p) {
      GnpOld = calcGnp(&e, chargePoints, chargeValues, npts, n, p);
      calcBnpAndCnpFromGnp(&prob, n, p, GnpOld, &BnpOld, &CnpOld);
      ierr = VecSetValues(GnpValsOld, 1, &count, &GnpOld, INSERT_VALUES);CHKERRQ(ierr);
      ierr = VecSetValues(BnpValsOld, 1, &count, &BnpOld, INSERT_VALUES);CHKERRQ(ierr);
      ierr = VecSetValues(CnpValsOld, 1, &count, &CnpOld, INSERT_VALUES);CHKERRQ(ierr);

      count++;
    }
  }
  ierr = VecAssemblyBegin(GnpValsOld);CHKERRQ(ierr); ierr = VecAssemblyEnd(GnpValsOld);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(BnpValsOld);CHKERRQ(ierr); ierr = VecAssemblyEnd(BnpValsOld);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(CnpValsOld);CHKERRQ(ierr); ierr = VecAssemblyEnd(CnpValsOld);CHKERRQ(ierr);

  PetscReal GnpErr, BnpErr, CnpErr;
  ierr = VecWAXPY(errorV, -1.0, GnpVals, GnpValsOld);CHKERRQ(ierr);
  ierr = VecNorm(errorV, NORM_2, &GnpErr);
  ierr = VecWAXPY(errorV, -1.0, BnpVals, BnpValsOld);CHKERRQ(ierr);
  ierr = VecNorm(errorV, NORM_2, &BnpErr);
  ierr = VecWAXPY(errorV, -1.0, CnpVals, CnpValsOld);CHKERRQ(ierr);
  ierr = VecNorm(errorV, NORM_2, &CnpErr);
  printf("GnpErr: %3.3e\n", GnpErr);
  printf("BnpErr: %3.3e\n", BnpErr);
  printf("CnpErr: %3.3e\n", CnpErr);

  //solution vector
  solution = (double*) malloc(sizeof(double)*npts);

  //calculate potential on solPoints
  calcCoulombEllipsoidalGrid(&prob, 10, solPoints, npts, solution);
  printf("wow!\n");
  for(int i=0; i<npts; ++i)
    printf("solution[%d] = %15.15f\n", i, solution[i]);
  for(int i=0; i<npts; ++i)
    printf("chargeValues[%d] = %15.15f\n", i, chargeValues[i]);

  //calculate free energy
  freeEnergy = 0;
  for(int i=0; i<npts; ++i) {
    freeEnergy += solution[i]*chargeValues[i];
    printf("freeEnergy[%d] = %15.15f\n", i, freeEnergy*.5*cf);
  }

  free(xyz);
  free(chargeValues);
  free(solPoints);
  free(solution);
  free(chargePoints);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CoulombExact"
PetscErrorCode CoulombExact(PetscReal eps1, Vec srcXYZ, Vec srcMag, Vec tarXYZ, Vec sol)
{
  PetscErrorCode ierr;
  PetscReal srcX, srcY, srcZ;
  PetscReal tarX, tarY, tarZ;
  PetscReal dx2, dy2, dz2;
  PetscReal val, mag;
  const PetscScalar *srcXYZArray;
  const PetscScalar *srcMagArray;
  const PetscScalar *tarXYZArray;
  PetscScalar *solArray;
  PetscInt nSrc, nTar;
  PetscInt i, k;
  PetscFunctionBeginUser;

  ierr = VecGetSize(srcMag, &nSrc);CHKERRQ(ierr);
  ierr = VecGetSize(sol, &nTar);CHKERRQ(ierr);

  ierr = VecGetArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecGetArray    (sol   , &solArray   );CHKERRQ(ierr);
  for(i=0; i<nTar; ++i) {
    /* get target X,Y,Z */
    tarX = tarXYZArray[3*i+0];
    tarY = tarXYZArray[3*i+1];
    tarZ = tarXYZArray[3*i+2];
    for(k=0; k<nSrc; ++k) {
      /* get source X,Y,Z */
      srcX = srcXYZArray[3*k+0];
      srcY = srcXYZArray[3*k+1];
      srcZ = srcXYZArray[3*k+2];
      mag  = srcMagArray[k];

      dx2 = (tarX-srcX)*(tarX-srcX);
      dy2 = (tarY-srcY)*(tarY-srcY);
      dz2 = (tarZ-srcZ)*(tarZ-srcZ);

      val = mag/(eps1*PetscSqrtReal(dx2 + dy2 + dz2));
      solArray[i] += val;
    }
  }
  ierr = VecRestoreArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecRestoreArray    (sol   , &solArray   );CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CoulombExactZeroInside"
PetscErrorCode CoulombExactZeroInside(PetscReal a, PetscReal b, PetscReal c, PetscReal eps1, Vec srcXYZ, Vec srcMag, Vec tarXYZ, Vec sol)
{
  PetscErrorCode ierr;
  PetscReal srcX, srcY, srcZ;
  PetscReal tarX, tarY, tarZ;
  PetscReal dx2, dy2, dz2;
  PetscReal val, mag;
  const PetscScalar *srcXYZArray;
  const PetscScalar *srcMagArray;
  const PetscScalar *tarXYZArray;
  PetscScalar *solArray;
  PetscInt nSrc, nTar;
  PetscInt i, k;
  PetscFunctionBeginUser;

  ierr = VecGetSize(srcMag, &nSrc);CHKERRQ(ierr);
  ierr = VecGetSize(sol, &nTar);CHKERRQ(ierr);

  ierr = VecZeroEntries(sol);CHKERRQ(ierr);

  ierr = VecGetArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecGetArray    (sol   , &solArray   );CHKERRQ(ierr);
  for(i=0; i<nTar; ++i) {
    /* get target X,Y,Z */
    tarX = tarXYZArray[3*i+0];
    tarY = tarXYZArray[3*i+1];
    tarZ = tarXYZArray[3*i+2];
    for(k=0; k<nSrc; ++k) {
      /* get source X,Y,Z */
      srcX = srcXYZArray[3*k+0];
      srcY = srcXYZArray[3*k+1];
      srcZ = srcXYZArray[3*k+2];
      mag  = srcMagArray[k];

      dx2 = (tarX-srcX)*(tarX-srcX);
      dy2 = (tarY-srcY)*(tarY-srcY);
      dz2 = (tarZ-srcZ)*(tarZ-srcZ);

      val = mag/(eps1*PetscSqrtReal(dx2 + dy2 + dz2));
      solArray[i] += val;
    }
    if((tarX*tarX)/(a*a) + (tarY*tarY)/(b*b) + (tarZ*tarZ)/(c*c) <= 1)
      solArray[i] = 0;
  }
  ierr = VecRestoreArrayRead(srcXYZ, &srcXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecRestoreArray    (sol   , &solArray   );CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "RandomEllipsoidPoints"
PetscErrorCode RandomEllipsoidPoints(PetscReal a, PetscReal b, PetscReal c, Vec xyz)
{
  PetscErrorCode ierr;
  PetscReal r, theta, phi;
  PetscReal x, y, z;
  const PetscReal theta_min = 0;
  const PetscReal theta_max = 2*PETSC_PI;
  const PetscReal phi_min   = 0;
  const PetscReal phi_max   = PETSC_PI;
  const PetscReal r_min     = 0;
  const PetscReal r_max     = 1;
  PetscScalar *xyzArray;
  PetscInt nvals;
  PetscFunctionBegin;

  /* petsc random object */
  PetscRandom rnd1, rnd2;
  ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rnd1);CHKERRQ(ierr);
  ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rnd2);CHKERRQ(ierr);
  ierr = PetscRandomSetInterval(rnd1, .1, 1.0);CHKERRQ(ierr);
  ierr = PetscRandomSetInterval(rnd2, .8, .9);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rnd1);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rnd2);CHKERRQ(ierr);
  srand(time(NULL));
  ierr = PetscRandomSetSeed(rnd1, rand());
  ierr = PetscRandomSetSeed(rnd2, rand());
  ierr = PetscRandomSeed(rnd1);CHKERRQ(ierr);
  ierr = PetscRandomSeed(rnd2);CHKERRQ(ierr);

  ierr = VecGetSize(xyz, &nvals);CHKERRQ(ierr);
  nvals = nvals/3;

  ierr = VecGetArray(xyz, &xyzArray);CHKERRQ(ierr);
  for(PetscInt i=0; i<nvals; ++i) {
    ierr = PetscRandomGetValue(rnd2, &r);CHKERRQ(ierr);
    ierr = PetscRandomGetValue(rnd1, &theta);CHKERRQ(ierr);
    ierr = PetscRandomGetValue(rnd1, &phi);CHKERRQ(ierr);

    //r = r_min + (r_max-r_min)*i/(nvals-1);
    //theta = theta_min + (theta_max-theta_min)*i/(nvals-1);
    //phi = phi_min + (phi_max-phi_min)*i/(nvals-1);
    r = r_min + (r_max-r_min)*r;
    theta = theta_min + (theta_max-theta_min)*theta;
    phi = phi_min + (phi_max-phi_min)*phi;

    x = a*r*PetscCosReal(theta)*PetscSinReal(phi);
    y = b*r*PetscSinReal(theta)*PetscSinReal(phi);
    z = c*r*PetscCosReal(phi);

    xyzArray[3*i+0] = x;
    xyzArray[3*i+1] = y;
    xyzArray[3*i+2] = z;
    if(nvals==1) {
      xyzArray[3*i+0] = .1;
      xyzArray[3*i+1] = .2;
      xyzArray[3*i+2] = .4;
    }
    /*if(nvals==2) {
      xyzArray[0] = .1;
      xyzArray[1] = .2;
      xyzArray[2] = .4;
      xyzArray[3] = .1;
      xyzArray[4] = .2;
      xyzArray[5] = .3;
      }*/
  }
  ierr = VecRestoreArray(xyz, &xyzArray);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "RandomMag"
PetscErrorCode RandomMag(Vec mags)
{
  PetscErrorCode ierr;
  PetscInt nMag, i;
  PetscReal val;
  PetscRandom rnd;
  PetscScalar *magsArray;
  PetscFunctionBegin;

  ierr = VecGetSize(mags, &nMag);CHKERRQ(ierr);

  /* petsc random object */
  ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rnd);CHKERRQ(ierr);
  ierr = PetscRandomSetInterval(rnd, -1.0, 1.0);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rnd);CHKERRQ(ierr);

  ierr = VecGetArray(mags, &magsArray);CHKERRQ(ierr);
  for(i=0; i<nMag; ++i) {
    ierr = PetscRandomGetValue(rnd, &val);CHKERRQ(ierr);
    magsArray[i] = val;
  }
  ierr = VecRestoreArray(mags, &magsArray);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SphereLimitExample"
PetscErrorCode SphereLimitExample()
{
  PetscErrorCode ierr;
  const PetscInt NUM_DELTAS = 14;
  const PetscReal DELTAS[14] = {1.3, 1.0, .5, .25, .15, .05, .01, .005, .001, .0005, .0001, .00005, .00001, .000005};
  const PetscInt  Nmax = 10;
  const PetscReal eps1 = 4.0;
  const PetscReal eps2 = 80.0;
  Vec srcXYZ[NUM_DELTAS], srcMag[NUM_DELTAS], solution[NUM_DELTAS];
  PetscScalar *srcXYZArray, *srcMagArray;
  PetscReal a, b, c;
  PetscReal delta;
  PetscReal *freeE;
  PetscInt i;
  PetscLogEvent events[NUM_DELTAS];
  char tText[40] = "Free energy with delta=%8.8f";
  char eText[40];
  EllipsoidalSystem ells[NUM_DELTAS];
  PetscFunctionBegin;

  ierr = PetscMalloc1(sizeof(PetscReal)*NUM_DELTAS, &freeE);

  for(i=0; i<NUM_DELTAS; ++i) {
    /* get size of delta and set a,b,c */
    delta = DELTAS[i];
    a = 1.0 + delta;
    b = 1.0 + (delta/5.0);
    c = 1.0 + (delta/10.0);
    /* init ellipsoidal system */
    ierr = initEllipsoidalSystem(ells+i, a, b, c, 64);CHKERRQ(ierr);
    /* initialize xyz and solution vectors */
    ierr = VecCreateSeq(PETSC_COMM_SELF, 3, srcXYZ+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, 1, srcMag+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, 1, solution+i);CHKERRQ(ierr);
    ierr = VecZeroEntries(srcXYZ[i]);CHKERRQ(ierr);
    ierr = VecGetArray(srcMag[i], &srcMagArray);CHKERRQ(ierr);
    srcMagArray[0] = 1.0; // set only charge to have magnitude 1.0
    ierr = VecRestoreArray(srcMag[i], &srcMagArray);CHKERRQ(ierr);
    /* register events for flop counting */
    sprintf(eText, tText, delta);
    ierr = PetscLogEventRegister((const char*) eText, 0, events+i);CHKERRQ(ierr);
  }
  /* solutions */
  for(i=0; i<NUM_DELTAS; ++i) {
    printf("calculating solution with delta=%8.8f\n", DELTAS[i]);
    ierr = PetscLogEventBegin(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = CalcEllipsoidFreeEnergy(ells+i, eps1, eps2, 1, srcXYZ[i], srcMag[i], 1e-5, Nmax, solution[i], freeE+i);CHKERRQ(ierr);
    printf("freeE: %15.15f\n", freeE[i]);

    ierr = PetscLogEventEnd(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
  }
  FILE *fp = fopen("out/sphConv.txt", "w");
  fprintf(fp, "delta error");
  PetscReal err;
  PetscReal exact = 1./eps2 - 1./eps1;
  for(i=0; i<NUM_DELTAS; ++i) {
    err = PetscAbsReal((exact - freeE[i])/exact);
    fprintf(fp, "%4.4e %4.4e\n", DELTAS[i], err);
    ierr = VecDestroy(srcXYZ+i);CHKERRQ(ierr);
    ierr = VecDestroy(srcMag+i);CHKERRQ(ierr);
    ierr = VecDestroy(solution+i);CHKERRQ(ierr);
  }
  fclose(fp);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "ChargeFlopsExample"
PetscErrorCode ChargeFlopsExample(PetscInt Nmax)
{
  PetscErrorCode ierr;

  const PetscInt NUM_SOLUTIONS = 10;
  const PetscInt CHARGE_START = 100;
  const PetscInt CHARGE_INC   = 50;
  Vec srcXYZ[NUM_SOLUTIONS], srcMag[NUM_SOLUTIONS], solution[NUM_SOLUTIONS];
  PetscInt i;
  PetscInt chargeNums[NUM_SOLUTIONS];
  PetscLogEvent events[NUM_SOLUTIONS];
  PetscInt charges;

  char tText[30] = "Free energy with %d charges";
  char eText[30];
  PetscInt tempNum;

  EllipsoidalSystem e;
  PetscFunctionBegin;

  const PetscReal eps1 = 4.0;
  const PetscReal eps2 = 4.0;

  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;

  ierr = initEllipsoidalSystem(&e, a, b, c, 64);CHKERRQ(ierr);

  /* create the vector with charge numbers to use */
  for(i=0; i < NUM_SOLUTIONS; ++i) {
    charges = CHARGE_START + CHARGE_INC*i;
    chargeNums[i] = charges;
    /* initialize xyz and solution vectors */
    ierr = VecCreateSeq(PETSC_COMM_SELF, 3*charges, srcXYZ+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, charges, srcMag+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, charges, solution+i);CHKERRQ(ierr);
    /* generate random points and magnitudes */
    ierr = RandomEllipsoidPoints(a, b, c, srcXYZ[i]);CHKERRQ(ierr);
    ierr = RandomMag(srcMag[i]);CHKERRQ(ierr);

    /* register events for flop counting */
    sprintf(eText, tText, charges);
    ierr = PetscLogEventRegister((const char*) eText, 0, events+i);CHKERRQ(ierr);
  }

  PetscReal freeE;
  for(i=0; i < NUM_SOLUTIONS; ++i) {
    printf("calculating solution %d/%d\n", i+1, NUM_SOLUTIONS);
    ierr = PetscLogEventBegin(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = CalcEllipsoidFreeEnergy(&e, eps1, eps2, chargeNums[i], srcXYZ[i], srcMag[i], 1e-5, Nmax, solution[i], &freeE);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
  }
  //ierr = EasyExample(Nmax, nSrc, nx, xl, xr, ny, yl, yr, zConst);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "WorkPrecExample"
PetscErrorCode WorkPrecExample()
{
  PetscErrorCode ierr;

  const PetscInt NUM_SOLUTIONS = 45;
  const PetscInt EXACT_NUM   = NUM_SOLUTIONS+10;
  const PetscInt NUM_CHARGES = 4;
  Vec srcXYZ, srcMag;
  Vec potential[NUM_SOLUTIONS];
  Vec potExact;
  Vec solution;
  PetscScalar *solutionArray;
  PetscReal exact;
  PetscInt i;
  PetscLogEvent events[NUM_SOLUTIONS];
  const PetscInt Nmin = 0;
  const PetscInt Nstep = 1;
  PetscInt N;
  char tText[30] = "Free energy with Nmax = %d";
  char eText[30];
  PetscInt tempNum;

  EllipsoidalSystem e;
  PetscFunctionBegin;

  const PetscReal eps1 = 4.0;
  const PetscReal eps2 = 80.0;

  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;

  ierr = initEllipsoidalSystem(&e, a, b, c, 64);CHKERRQ(ierr);

  /* create charge vectors and initialize with random values */
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*NUM_CHARGES, &srcXYZ);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,   NUM_CHARGES, &srcMag);CHKERRQ(ierr);
  ierr = RandomEllipsoidPoints(a, b, c, srcXYZ);CHKERRQ(ierr);
  printf("######\nXYZ\n#######\n");
  ierr = VecView(srcXYZ, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  ierr = RandomMag(srcMag);CHKERRQ(ierr);
  ierr = VecSet(srcMag, 1.0);CHKERRQ(ierr);
  //ierr = VecScale(srcMag, 1./NUM_CHARGES);CHKERRQ(ierr); //scale charge magnitudes so total energy is more normalized

  for(i=0; i<NUM_SOLUTIONS; ++i) {
    /* register events for flop counting */
    N = Nmin + Nstep*i;
    sprintf(eText, tText, N);
    ierr = PetscLogEventRegister((const char*) eText, 0, events+i);CHKERRQ(ierr);
    /* create potential vectors */
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_CHARGES, potential+i);CHKERRQ(ierr);
  }
  /* create free energy solution vector */
  ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOLUTIONS, &solution);CHKERRQ(ierr);

  ierr = VecGetArray(solution, &solutionArray);CHKERRQ(ierr);
  /* approximate solutions */
  for(i=0; i < NUM_SOLUTIONS; ++i) {
    N = Nmin + Nstep*i;
    printf("calculating solution %d/%d\n", i+1, NUM_SOLUTIONS);
    printf("NNNNN = %d\n", N);
    ierr = PetscLogEventBegin(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = CalcEllipsoidFreeEnergy(&e, eps1, eps2, NUM_CHARGES, srcXYZ, srcMag, 1e-16, N, potential[i], solutionArray+i);CHKERRQ(ierr);
    //ierr = VecView(solution, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
    ierr = VecView(potential[i], PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(events[i], 0, 0, 0, 0);CHKERRQ(ierr);
  }

  /* exact(ish) solution */
  ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_CHARGES, &potExact);CHKERRQ(ierr);
  printf("calculating exact solution\n");
  ierr = CalcEllipsoidFreeEnergy(&e, eps1, eps2, NUM_CHARGES, srcXYZ, srcMag, 1e-16, EXACT_NUM, potExact, &exact);CHKERRQ(ierr);

  FILE *fp = fopen("out/workprec.txt", "w");
  fprintf(fp, "work   precision\n");
  PetscReal errVal;
  PetscEventPerfInfo info;
  //* calculate error for approximate solutions */
  for(PetscInt i=0; i < NUM_SOLUTIONS; ++i) {
    errVal = PetscAbsReal((solutionArray[i] - exact)/(exact));
    ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, events[i], &info);CHKERRQ(ierr);
    fprintf(fp, "%4.4e %4.4e\n", info.flops, errVal);
    printf("approx: %15.15f\n", solutionArray[i]);
  }
  printf("exact: %15.15f\n", exact);
  fclose(fp);
  ierr = VecRestoreArray(solution, &solutionArray);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "InitGrid"
PetscErrorCode InitGrid(PetscInt nx, PetscReal xl, PetscReal xr, PetscInt ny, PetscReal yl, PetscReal yr, PetscReal zConst, Vec *xyz)
{
  PetscErrorCode ierr;
  PetscInt i, j;
  PetscReal xh, yh;
  PetscInt index;
  PetscScalar *xyzArray;
  PetscFunctionBegin;

  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nx*ny, xyz);CHKERRQ(ierr);

  xh = (xr - xl)/(nx - 1);
  yh = (yr - yl)/(ny - 1);

  ierr = VecGetArray(*xyz, &xyzArray);CHKERRQ(ierr);
  for(i=0; i < nx; ++i) {
    for(j=0; j<ny; ++j) {
      index = 3*(i*ny + j);
      xyzArray[index+0] = xl + xh*i;
      xyzArray[index+1] = yl + yh*j;
      xyzArray[index+2] = zConst;
    }
  }
  ierr = VecRestoreArray(*xyz, &xyzArray);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "GridSolution"
/*
  Example is of a some charges
*/
PetscErrorCode GridSolution(PetscInt Nmax, PetscInt nSrc, PetscInt nx, PetscReal xl, PetscReal xr, PetscInt ny, PetscReal yl, PetscReal yr, PetscReal zConst, Vec *xvals, Vec *yvals, Vec *sol)
{
  PetscErrorCode ierr;
  Vec sourceXYZ, sourceMag, tarXYZ, solution;
  PetscReal xh, yh;
  PetscReal xc, yc;
  PetscReal randVal;
  PetscInt ind;
  PetscFunctionBegin;

  const PetscReal eps1 = 4.0;
  const PetscReal eps2 = 4.0;

  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;

  // initialize XYZ and solution vectors
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nSrc , &sourceXYZ);CHKERRQ(ierr);
  //ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nx*ny, &tarXYZ);   CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nx*ny  , &solution); CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nSrc   , &sourceMag);CHKERRQ(ierr);

  xh = (xr - xl) / (nx - 1);
  yh = (yr - yl) / (ny - 1);

  ierr = RandomEllipsoidPoints(a, b, c, sourceXYZ);CHKERRQ(ierr);

  // populate source vec with points
  for(PetscInt k=0; k<nSrc; ++k) {
    ind = 3*k;
    ierr = VecSetValue(sourceMag, k  , 1.0      , INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(sourceMag);CHKERRQ(ierr); ierr = VecAssemblyEnd(sourceMag);CHKERRQ(ierr);

  ierr = InitGrid(nx, xl, xr, ny, yl, yr, zConst, &tarXYZ);

  // calculate the solvation potential
  ierr = CalcEllipsoidTester( a , b , c,
			      eps1, eps2,
			      nSrc, sourceXYZ, sourceMag,
			      nx*ny, tarXYZ,
			      Nmax, solution); CHKERRQ(ierr);

  /*
  FILE *fp = fopen("solvPlot.txt", "w");
  fprintf(fp, "xh: %4.4e\n", xh);
  fprintf(fp, "yh: %4.4e\n", yh);
  const PetscScalar *solutionArray;
  ierr = VecGetArrayRead(solution, &solutionArray);CHKERRQ(ierr);
  PetscInt index = 0;
  for(PetscInt i=0; i<nx; ++i) {
    for(PetscInt j=0; j<ny; ++j) {
      index = i*ny + j;
      fprintf(fp, "%4.4f ", solutionArray[index]);
    }
    fprintf(fp, "\n");
  }
  ierr = VecRestoreArrayRead(solution, &solutionArray);CHKERRQ(ierr);
  */
  if(xvals != NULL) {
    ierr = VecCreateSeq(PETSC_COMM_WORLD, nx, xvals);CHKERRQ(ierr);
    for(PetscInt i=0; i<nx; ++i)
      ierr = VecSetValue(*xvals, i, xl + i*xh, INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(*xvals);CHKERRQ(ierr); ierr = VecAssemblyEnd(*xvals);CHKERRQ(ierr);
  }
  if(yvals != NULL) {
    ierr = VecCreateSeq(PETSC_COMM_WORLD, nx, yvals);CHKERRQ(ierr);
    for(PetscInt i=0; i<ny; ++i)
      ierr = VecSetValue(*yvals, i, yl + i*yh, INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(*yvals);CHKERRQ(ierr); ierr = VecAssemblyEnd(*yvals);CHKERRQ(ierr);
  }
  if(sol != NULL) {
    ierr = VecDuplicate(solution, sol);CHKERRQ(ierr);
    ierr = VecCopy(solution, *sol);CHKERRQ(ierr);
  }
  ierr = VecDestroy(&sourceXYZ);CHKERRQ(ierr);
  ierr = VecDestroy(&tarXYZ);CHKERRQ(ierr);
  ierr = VecDestroy(&solution);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SubtractCoulombInside"
PetscErrorCode SubtractCoulombInside(PetscReal a, PetscReal b, PetscReal c, Vec tarXYZ, Vec solution, Vec solutionEx, Vec ncSolution)
{
  PetscErrorCode ierr;
  const PetscScalar *tarXYZArray, *solutionArray, *solutionExArray;
  PetscScalar *ncSolutionArray;
  PetscInt nTar;
  PetscInt i;
  PetscReal x, y, z;
  PetscFunctionBegin;

  ierr = VecGetSize(solution, &nTar);CHKERRQ(ierr);

  ierr = VecGetArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(solution, &solutionArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(solutionEx, &solutionExArray);CHKERRQ(ierr);
  ierr = VecGetArray(ncSolution, &ncSolutionArray);CHKERRQ(ierr);

  for(i=0; i<nTar; ++i) {
    x = tarXYZArray[3*i+0];
    y = tarXYZArray[3*i+1];
    z = tarXYZArray[3*i+2];

    if( ((x*x)/(a*a)) + ((y*y)/(b*b)) + ((z*z)/(c*c)) <= 1) {
      ncSolutionArray[i] = solutionArray[i];
    }
    else {
      ncSolutionArray[i] = solutionArray[i] - solutionExArray[i];
    }
  }

  ierr = VecRestoreArrayRead(tarXYZ, &tarXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(solution, &solutionArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(solutionEx, &solutionExArray);CHKERRQ(ierr);
  ierr = VecRestoreArray(ncSolution, &ncSolutionArray);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "GridAnimation"
PetscErrorCode GridAnimation()
{
  const PetscInt NUM_SOLUTIONS = 8;
  PetscErrorCode ierr;
  Vec sourceXYZ, sourceMag, tarXYZ;
  Vec solution[NUM_SOLUTIONS], solutionEx[NUM_SOLUTIONS], ncSolution[NUM_SOLUTIONS];
  Vec errorVec[NUM_SOLUTIONS];
  Vec xOut, yOut, zOut;
  PetscReal xh, yh;
  PetscReal xc, yc;
  PetscReal randVal;
  PetscInt ind;
  PetscFunctionBegin;

  PetscReal eps1 = 2.0;
  PetscReal eps2 = 20.0;
  PetscInt nSrc  = 5;
  PetscInt nx    = 15;
  PetscReal xl   = -5.13;
  PetscReal xr   = 5.13;
  PetscInt ny    = 15;
  PetscReal yl   = -4.02;
  PetscReal yr   =  4.02;
  PetscReal zConst = .1;

  xh = (xr-xl)/ny;
  yh = (yr-yl)/ny;

  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;

  // initialize source XYZ, mag
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nSrc , &sourceXYZ);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nSrc   , &sourceMag);CHKERRQ(ierr);
  // initialize solution vectors
  for(PetscInt i=0; i < NUM_SOLUTIONS; ++i) {
    ierr = VecCreateSeq(PETSC_COMM_SELF, nx*ny  , solution+i); CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, nx*ny  , solutionEx+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, nx*ny  , ncSolution+i);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, nx*ny  , errorVec+i);CHKERRQ(ierr);

  }

  ierr = RandomEllipsoidPoints(a, b, c, sourceXYZ);CHKERRQ(ierr); // on z=0

  //ierr = RandomMag(sourceMag); CHKERRQ(ierr);
  ierr = VecSet(sourceMag, 1.0);CHKERRQ(ierr);

  ierr = InitGrid(nx, xl, xr, ny, yl, yr, zConst, &tarXYZ);CHKERRQ(ierr); //zConst = 0

  //Vec tester;
  //ierr = VecDuplicate(tarXYZ, &tester);
  PetscReal err;
  // calculate the solvation potential
  for(PetscInt i=0; i < NUM_SOLUTIONS; ++i) {
    printf("wow\n");
    ierr = CoulombExact(eps1, sourceXYZ, sourceMag, tarXYZ, solutionEx[i]);
    //ierr = VecCopy(tarXYZ, tester);
    printf("eps: %15.15f\n", eps1);
    ierr = CalcEllipsoidTester( a , b , c,
				eps1, eps2,
				nSrc, sourceXYZ, sourceMag,
				nx*ny, tarXYZ,
				i, solution[i]); CHKERRQ(ierr);
    ierr = SubtractCoulombInside(a, b, c, tarXYZ, solution[i], solutionEx[i], ncSolution[i]);CHKERRQ(ierr);
    PetscInt zero = 0;
    PetscReal sol0, solEx0;
    ierr = VecGetValues(solution[i], 1, &zero, &sol0);CHKERRQ(ierr);
    ierr = VecGetValues(solutionEx[i], 1, &zero, &solEx0);CHKERRQ(ierr);
    //ierr = VecAXPY(tester, -1.0, tarXYZ);CHKERRQ(ierr);
    //ierr = VecView(solutionEx[i], PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
    //ierr = VecView(solution[i], PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
    ierr = VecCopy(solutionEx[i], errorVec[i]);CHKERRQ(ierr);
    ierr = VecAXPY(errorVec[i], -1.0, solution[i]);
    ierr = VecNorm(errorVec[i], NORM_2, &err);
    printf("\n");
    printf("the error on %d is %15.15f\n", i, err*xh*yh);
  }

  /* subtract coulomb from solution inside */

  char fname[20] = "out/sol%d.txt";
  char fnameEx[20] = "out/solEx%d.txt";
  char fnameErr[20] = "out/solErr%d.txt";
  char fnameS[20];
  char fnameExS[20];
  char fnameErrS[20];
  for(PetscInt i=0; i < NUM_SOLUTIONS; ++i) {
    sprintf(fnameS, fname, i);
    sprintf(fnameExS, fnameEx, i);
    sprintf(fnameErrS, fnameErr, i);
    ierr = WriteToFile((const char*) fnameS, ny, ncSolution[i]);CHKERRQ(ierr);
    ierr = WriteToFile((const char*) fnameExS, ny, solutionEx[i]);CHKERRQ(ierr);
    ierr = WriteToFile((const char*) fnameErrS, ny, errorVec[i]);CHKERRQ(ierr);
  }

  ierr = VecCreateSeq(PETSC_COMM_SELF, nx, &xOut);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, ny, &yOut);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, 1, &zOut);CHKERRQ(ierr);
  xh = (xr - xl) / (nx - 1);
  yh = (yr - yl) / (ny - 1);
  for(PetscInt i=0; i<nx; ++i) {
    ierr = VecSetValue(xOut, i, xl + xh*i, INSERT_VALUES);CHKERRQ(ierr);
  }
  for(PetscInt i=0; i<ny; ++i) {
    ierr = VecSetValue(yOut, i, yl + yh*i, INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecSetValue(zOut, 0, zConst, INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(xOut);CHKERRQ(ierr); ierr = VecAssemblyEnd(xOut);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(yOut);CHKERRQ(ierr); ierr = VecAssemblyEnd(yOut);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(zOut);CHKERRQ(ierr); ierr = VecAssemblyEnd(zOut);CHKERRQ(ierr);

  ierr = WriteToFile((const char*) "out/xVals.txt", 1, xOut);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/yVals.txt", 1, yOut);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/zVals.txt", 1, zOut);CHKERRQ(ierr);

  /* output source xyz values */
  ierr = WriteToFile((const char*) "out/chargeXYZ.txt", 3, sourceXYZ);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/chargeMag.txt", 1, sourceMag);CHKERRQ(ierr);
  FILE *fp = fopen("out/otherinfo.txt", "w");
  fprintf(fp, "%6.6f %6.6f %6.6f\n", a, b, c);
  fclose(fp);

  ierr = VecDestroy(&sourceXYZ);CHKERRQ(ierr);
  ierr = VecDestroy(&tarXYZ);CHKERRQ(ierr);
  for(PetscInt i=0; i < NUM_SOLUTIONS; ++i) {
    ierr = VecDestroy(solution+i);CHKERRQ(ierr);
    ierr = VecDestroy(solutionEx+i);CHKERRQ(ierr);
    ierr = VecDestroy(ncSolution+i);CHKERRQ(ierr);
  }
  ierr = VecDestroy(&sourceMag);CHKERRQ(ierr);
  ierr = VecDestroy(&xOut);CHKERRQ(ierr);
  ierr = VecDestroy(&yOut);CHKERRQ(ierr);
  ierr = VecDestroy(&zOut);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "numChargesPlot"
PetscErrorCode numChargesPlot(PetscInt nMin, PetscInt nMax, PetscInt nStep)
{
  PetscErrorCode ierr;
  PetscInt Nmax = 4;
  //PetscInt nSrc = 20;
  PetscInt nx   = 7;
  PetscReal xl  = -4.6;
  PetscReal xr  = 4.6;
  PetscInt ny   = 7;
  PetscReal yl  = -4.6;
  PetscReal yr  = 4.6;
  PetscReal zConst = .1;
  PetscInt nSource;
  PetscLogStage *chargeStages;
  PetscStageLog stageLog;
  PetscFunctionBegin;

  PetscInt count = 0;
  for(nSource = nMin; nSource <= nMax; nSource += nStep) {
    count++;
  }
  ierr = PetscMalloc1(sizeof(PetscLogStage)*count, &chargeStages);
  char bStr[40] = "Num Charges: %d";
  char indvStr[40];
  count = 0;
  for(nSource = nMin; nSource <= nMax; nSource += nStep) {
    sprintf(indvStr, bStr, nSource);
    ierr = PetscLogStageRegister(indvStr, chargeStages+count);CHKERRQ(ierr);
    count++;
  }
  count = 0;
  for(nSource = nMin; nSource <= nMax; nSource += nStep) {
    printf("wow\n");
    ierr = PetscLogStagePush(chargeStages[count]);CHKERRQ(ierr);
    ierr = GridSolution(Nmax, nSource, nx, xl, xr, ny, yl, yr, zConst, NULL, NULL, NULL);CHKERRQ(ierr);
    ierr = PetscLogStagePop();CHKERRQ(ierr);
    count++;
  }

  ierr = PetscLogGetStageLog(&stageLog);CHKERRQ(ierr);
  count = 0;
  for(nSource = nMin; nSource <= nMax; nSource += nStep) {
    ierr = PetscPrintf(PETSC_COMM_WORLD, "n = %d\nflops = %4.4e\n", nSource, stageLog->stageInfo[chargeStages[count]].perfInfo.flops);CHKERRQ(ierr);
    count++;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "WriteToFile"
PetscErrorCode WriteToFile(const char *fname, PetscInt rowsize, Vec values)
{
  PetscErrorCode ierr;
  PetscInt k;
  PetscInt nPts;
  const PetscScalar *valuesArray;
  PetscFunctionBegin;

  ierr = VecGetSize(values, &nPts);CHKERRQ(ierr);

  ierr = VecGetArrayRead(values, &valuesArray);CHKERRQ(ierr);
  FILE *fp = fopen(fname, "w");
  for(k=0; k < nPts; ++k) {
    fprintf(fp, "%4.4f ", valuesArray[k]);
    if((k+1)%rowsize == 0)
      fprintf(fp, "\n");
  }
  fclose(fp);
  ierr = VecRestoreArrayRead(values, &valuesArray);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SolutionAnimation"
PetscErrorCode SolutionAnimation()
{
  PetscErrorCode ierr;
  Vec x1, y1, s0, s1, s2, s3;
  PetscFunctionBegin;
  PetscInt nSrc = 10;
  PetscInt nx   = 15;
  PetscReal xl  = -4.6;
  PetscReal xr  = 4.6;
  PetscInt ny   = 15;
  PetscReal yl  = -4.6;
  PetscReal yr  = 4.6;
  PetscReal zConst = .1;

  ierr = GridSolution(0, nSrc, nx, xl, xr, ny, yl, yr, zConst, &x1, &y1, &s0);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/sol0.txt", ny, s0);CHKERRQ(ierr);
  ierr = GridSolution(1, nSrc, nx, xl, xr, ny, yl, yr, zConst, &x1, &y1, &s1);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/sol1.txt", ny, s1);CHKERRQ(ierr);
  ierr = GridSolution(2, nSrc, nx, xl, xr, ny, yl, yr, zConst, &x1, &y1, &s2);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/sol2.txt", ny, s2);CHKERRQ(ierr);
  ierr = GridSolution(3, nSrc, nx, xl, xr, ny, yl, yr, zConst, &x1, &y1, &s3);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/sol3.txt", ny, s3);CHKERRQ(ierr);

  ierr = WriteToFile((const char*) "out/xVals.txt", 1, x1);CHKERRQ(ierr);
  ierr = WriteToFile((const char*) "out/yVals.txt", 1, y1);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "testLame"
PetscErrorCode testLame()
{
  PetscErrorCode ierr;
  PetscInt n, p;
  PetscInt Nmax;
  PetscReal Fnp;
  PetscReal Enp;
  PetscReal a, b, c;
  PetscReal cval;
  PetscInt signm, signn;
  EllipsoidalSystem e;
  PetscFunctionBegin;

  cval = .3;
  signm = -1;
  signn = 1;

  a = 3.0; b = 2.0; c = 1.0;
  initEllipsoidalSystem(&e, a, b, c, 64);
  FILE *fp = fopen("newellout.txt", "w");
  Nmax = 6;
  for(n=0; n < Nmax; ++n) {
    for(p=0; p < 2*n+1; ++p) {
      ierr = calcLame(&e, n, p, cval, signm, signn, &Enp);CHKERRQ(ierr);
      fprintf(fp, "%d %d %15.15f\n", n, p, Enp);
    }
  }
  fclose(fp);

  PetscFunctionReturn(0);
}

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

  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  //initialize Petsc
  ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr);
  ierr = PetscLogDefaultBegin(); CHKERRQ(ierr);

  //ierr = EllVsSphConvergence();CHKERRQ(ierr);

  //testSphericalCoordinates();
  //runTest1();
  //testLegendre();
  //runTest4();
  //test4sphere();
  //RunArg();
  //ierr = RunArgTester(); CHKERRQ(ierr);
  //PetscErrorCode GridSolution(PetscInt Nmax, PetscInt nSrc, PetscInt nx, PetscReal xl, PetscReal xr, PetscInt ny, PetscReal yl, PetscReal yr, PetscReal zConst)
  /*
  PetscInt Nmax = 4;
  PetscInt nSrc = 25;
  PetscInt nx   = 15;
  PetscReal xl  = -5.23;
  PetscReal xr  = 5.23;
  PetscInt ny   = 15;
  PetscReal yl  = -4.45;
  PetscReal yr  = 4.45;
  PetscReal zConst = .1;
  PetscReal eps1 = 4.0;
  PetscReal eps2 = 80.0;
  */
  //ierr = GridSolution(Nmax, nSrc, nx, xl, xr, ny, yl, yr, zConst, NULL, NULL);CHKERRQ(ierr);
  //ierr = SolutionAnimation(); //<---- total crap
  //runTest1();
  //testLame();

//  ierr = runTest1();CHKERRQ(ierr);
  //ierr = RunArg();
  //ierr = GridAnimation();CHKERRQ(ierr);
  //ierr = ChargeFlopsExample(10);CHKERRQ(ierr);
  ierr = WorkPrecExample();CHKERRQ(ierr);
  //ierr = SphereLimitExample();CHKERRQ(ierr);

  //ierr = numChargesPlot(100, 500, 100);
  ierr = PetscFinalize();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Andreas Bertsch $
// $Authors: $
// --------------------------------------------------------------------------
//
#ifndef OPENMS_MATH_STATISTICS_GAMMADISTRIBUTIONFITTER_H
#define OPENMS_MATH_STATISTICS_GAMMADISTRIBUTIONFITTER_H

#include <OpenMS/DATASTRUCTURES/String.h>
#include <OpenMS/DATASTRUCTURES/DPosition.h>

#include <vector>

// gsl includes
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlin.h>

namespace OpenMS
{
  namespace Math
  {
    /**
      @brief Implements a fitter for the Gamma distribution.

      This class fits a Gamma distribution to a number of data points.
      The results as well as the initial guess are specified using the struct
          GammaDistributionFitResult.

      The formula with the fitted parameters can be transformed into a
      gnuplot formula using getGnuplotFormula() after fitting.

          The implementation is done using GSL fitting algorithms.

          @ingroup Math
      */
    class OPENMS_DLLAPI GammaDistributionFitter
    {
public:

      /// struct to represent the parameters of a gamma distribution
      struct GammaDistributionFitResult
      {
public:

        GammaDistributionFitResult() :
          b(1.0),
          p(5.0)
        {
        }

        GammaDistributionFitResult(const GammaDistributionFitResult & rhs) :
          b(rhs.b),
          p(rhs.p)
        {
        }

        GammaDistributionFitResult & operator=(const GammaDistributionFitResult & rhs)
        {
          if (this != &rhs)
          {
            b = rhs.b;
            p = rhs.p;
          }
          return *this;
        }

        /// parameter b of the gamma distribution
        double b;

        /// parameter p of the gamma distribution
        double p;
      };

      /// Default constructor
      GammaDistributionFitter();
      /// Destructor
      virtual ~GammaDistributionFitter();

      /// sets the gamma distribution start parameters b and p for the fitting
      void setInitialParameters(const GammaDistributionFitResult & result);

      /**
          @brief Fits a gamma distribution to the given data points

          @param points Input parameter which represents the point used for the fitting

          @exception Exception::UnableToFit is thrown if fitting cannot be performed
      */
      GammaDistributionFitResult fit(std::vector<DPosition<2> > & points);

      /// returns the gnuplot formula of the fitted gamma distribution
      const String & getGnuplotFormula() const;

protected:

      static int gammaDistributionFitterf_(const gsl_vector * x, void * params, gsl_vector * f);

      static int gammaDistributionFitterdf_(const gsl_vector * x, void * params, gsl_matrix * J);

      static int gammaDistributionFitterfdf_(const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J);

      void printState_(size_t iter, gsl_multifit_fdfsolver * s);

      GammaDistributionFitResult init_param_;

      String gnuplot_formula_;

private:
      /// Copy constructor (not implemented)
      GammaDistributionFitter(const GammaDistributionFitter & rhs);
      /// assignment operator (not implemented)
      GammaDistributionFitter & operator=(const GammaDistributionFitter & rhs);
    };
  }
}

#endif

ext: h

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
#pragma once
#include "utility/types.h"
#include <gsl/gsl>
#include <array>

namespace cws80 {

extern std::array<gsl::span<const u8>, 11> extra_bank_data;

}  // namespace cws80


text:
/* Simple Travelling Salesman Problem heuristics.
 *
 * Krzysztof Chalupka, 2017.
 */
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include "../utils/utils.h"
#include "../datastructures/datastructures.h"

int *tsp_bs(gsl_vector **, int); /* brute-force search */
int *tsp_nn(gsl_vector **, int); /* nearest-neighbor */
int *tsp_cp(gsl_vector **, int); /* closest-pair */
double pathlen(gsl_vector **, int *, int); /* path length of a permutation. */

double pathlen(gsl_vector **data, int *path, int n){
  /* Compute the length of the path defined by the vectors
   * in `data`, traversed in order given in `path`,
   * and closing the loop by traversing path[n-1] -- path[0]
   * at the end.
   */
  int i;
  double total = 0.;
  for (i = 0; i < n-1; i++)
    total += vector_dist(data[path[i]], data[path[i+1]]);
  if (n > 0)
    total += vector_dist(data[path[n-1]], data[path[0]]);
  return total;
}

int *tsp_nn(gsl_vector **data, int n)
{
  /* Solve the TSP using the nearest-neighbor heuristic. */
  if (n == 0)
    return NULL;
  int i, j, best_id;
  int usedup = 0; /* Flag already-used elements. */
  gsl_vector *current = data[0];
  double cur_dist;
  double best_dist;
  int *nn_path = (int *) calloc(n, sizeof(int));

  nn_path[0] = 0;
  usedup |= 1;
  for (i = 1; i < n; i++){
    best_dist = 0;
    for (j = 1; j < n; j++){
      if (!(usedup & (2 << (j-1)))){
    	cur_dist = vector_dist(current, data[j]);
    	if (cur_dist < best_dist || best_dist == 0){
    	  best_dist = cur_dist;
	  best_id = j;
    	}
      }
    }
    current = data[best_id];
    nn_path[i] = best_id;
    usedup |= (2 << (best_id-1));
  }

  return nn_path;
}

int *tsp_cp(gsl_vector **data, int n)
{
  /* Solve the TSP using the closest pair heuristic. */
  int i, j, ch1_id, ch2_id, best_ch1, best_ch2;
  int match_flag;
  double dist, best_dist;
  dlist_node *node;
  int *best_path = (int *) calloc(n, sizeof(int));

  /* Create an array where each element corresponds to a chain. */
  dlist **chains = (dlist **) calloc(n, sizeof(dlist *));
  for (i = 0; i < n; i++){
    chains[i] = dlist_new(data[i]);
  }

  /* Merge the chains according to the heuristic. */
  for (i = 0; i < n-1; i++){
    /* Find the smallest-distance chain ends pair. */
    best_dist = 0.;
    for (ch1_id = 0; ch1_id < n; ch1_id++){
      if (chains[ch1_id] == NULL){
	continue;
      }
      for (ch2_id = 0; ch2_id < n; ch2_id++){
	if (chains[ch2_id] == NULL || ch1_id == ch2_id){
	  continue;
	}
	dist = vector_dist(chains[ch1_id]->head->item,
			   chains[ch2_id]->head->item);
	if (dist < best_dist || best_dist == 0){
	  /* Heads match. */
	  best_dist = dist;
	  best_ch1 = ch1_id;
	  best_ch2 = ch2_id;
	  match_flag = 0;
	}

	dist = vector_dist(chains[ch1_id]->head->item,
			   chains[ch2_id]->tail->item);
	if (dist < best_dist || best_dist == 0){
	  /* Head1 matches tail2. */
	  best_dist = dist;
	  best_ch1 = ch1_id;
	  best_ch2 = ch2_id;
	  match_flag = 1;
	}

	dist = vector_dist(chains[ch1_id]->tail->item,
			   chains[ch2_id]->head->item);
	if (dist < best_dist || best_dist == 0){
	  /* Tail1 matches head2. */
	  best_dist = dist;
	  best_ch1 = ch1_id;
	  best_ch2 = ch2_id;
	  match_flag = 2;
	}

	dist = vector_dist(chains[ch1_id]->tail->item,
			   chains[ch2_id]->head->item);
	if (dist < best_dist || best_dist == 0){
	  /* Tail1 matches tail2. */
	  best_dist = dist;
	  best_ch1 = ch1_id;
	  best_ch2 = ch2_id;
	  match_flag = 3;
	}
      }
    }
    printf("Iter %d, best match (%d, %d) %f\n",
	   i, best_ch1, best_ch2, best_dist);

    /* Merge the smallest-distance chain ends. */
    if (match_flag == 0){
      /* Heads match: push ch2 to the front of ch1 in reverse-order. */
      while(chains[best_ch2]->head){
	dlist_insert(chains[best_ch1],
		     chains[best_ch2]->head->item);
	dlist_delete(chains[best_ch2], chains[best_ch2]->head);
      }
    }
    else if (match_flag == 1){
      /* Head1 to tail2 match: push ch2 to the front of ch1 in-order. */
      while(chains[best_ch2]->tail){
	dlist_insert(chains[best_ch1],
		     chains[best_ch2]->tail->item);
	dlist_delete(chains[best_ch2], chains[best_ch2]->tail);
      }
    }
    else if (match_flag == 2){
      /* Head2 to tail1 match: push ch2 to the end of ch1 in-order. */
      while(chains[best_ch2]->head){
	dlist_insert_tail(chains[best_ch1],
			  chains[best_ch2]->head->item);
	dlist_delete(chains[best_ch2], chains[best_ch2]->head);
      }
    }
    else if (match_flag == 3){
      /* Tails match: push ch2 to the end of ch1 in reverse-order. */
      while(chains[best_ch2]->tail){
	dlist_insert_tail(chains[best_ch1],
			  chains[best_ch2]->tail->item);
	dlist_delete(chains[best_ch2], chains[best_ch2]->tail);
      }
    }
    chains[best_ch2] = NULL;
  }

  /* All the chains except for one should be null. The non-null
     chain contains the best (heuristic) path. */
  for (i = 0; i < n; i++)
    if (chains[i])
      node = chains[i]->head;

  for(i = 0; node; i++){
    for (j = 0; j < n; j++)
      if (data[j] == node->item)
	best_path[i] = j;
    node = node->next;
  }
  return best_path;
}

int *tsp_bs(gsl_vector **data, int n)
{
  /* Find the shortest loop between `n` points. stored in `data`. */
  if (n == 0)
    return NULL;
  double best_dist = 0;
  double cur_dist = 0;
  int i_perm, i_iter;
  int *cur_path;
  int *best_path = calloc(n, sizeof(int));

  int **perm = permutations(n);
  for (i_iter = 0; i_iter < factorial(n); i_iter++){
    /* Retrieve the permutation into an array. */
    cur_path = perm[i_iter];

    /* Compute the path length corresponding to this permutation. */
    printf("Iter %d, permutation ", i_iter);
    print_intarray(cur_path, n);
    cur_dist = pathlen(data, cur_path, n);
    printf(" dist %.10f.\n", cur_dist);

    /* Save best pathlen to-date. */
    if (best_dist > cur_dist || best_dist == 0){
      best_dist = cur_dist;
      for (i_perm = 0; i_perm < n; i_perm++)
	best_path[i_perm] = cur_path[i_perm];
    }
  }
  printf("\n");
  free(cur_path);
  free(perm);
  return best_path;
}

int main(){
  char *fname;
  int npoints;
  int *bestpath;
  gsl_vector **vecs;

  /* Load data into a vector array. */
  fname = "tsp.dat";
  npoints = get_nlines(fname);
  vecs = load_vectors_from_file("tsp.dat", npoints);

  /* Compute the best path using the brute-force method. */
  bestpath = tsp_bs(vecs, npoints);
  printf("Best path (brute-force) len = %f, for path ",
	 pathlen(vecs, bestpath, npoints));
  print_intarray(bestpath, npoints);
  printf("\n");

  /* Compute the best path using the nearest-neighbor heuristic. */
  bestpath = tsp_nn(vecs, npoints);
  printf("Best path (nearest-neighbor) len = %f, for path ",
	 pathlen(vecs, bestpath, npoints));
  print_intarray(bestpath, npoints);
  printf("\n");

  /* Compute the best path using the closest-pairs heuristic. */
  bestpath = tsp_cp(vecs, npoints);
  printf("Best path (closest-pair) len = %f, for path ",
	 pathlen(vecs, bestpath, npoints));
  print_intarray(bestpath, npoints);
  printf("\n");

  return 0;
}

ext: c


text:
#include "lah.h"

#ifdef HAVE_LAPACK /* Use a LAPACK package to do the heavy work */
#include <cblas.h>

lah_Return lah_matMul(lah_MatOp transA, lah_MatOp transB,
                      lah_value alpha, lah_value beta,
                      lah_mat *C, lah_mat const *A, lah_mat const *B)
{

    CBLAS_TRANSPOSE transA_cblas = CblasNoTrans;
    CBLAS_TRANSPOSE transB_cblas = CblasNoTrans;

    if (C == NULL || A == NULL || B == NULL )
        return lahReturnParameterError;

    /*
    order = LAH_CBLAS_LAYOUT(C);

    if (order == CblasColMajor)
    {
        lda = A->nR;
        ldb = B->nR;
        ldc = C->nR;
    }
    else
    {
        lda = A->nC;
        ldb = B->nC;
        ldc = C->nC;
    }
    */

    if (transA != lahNorm)
        transA_cblas = CblasTrans;
    if (transB != lahNorm)
        transB_cblas = CblasTrans;

    GEMM(LAH_CBLAS_LAYOUT, transA_cblas, transB_cblas, C->nR, C->nC, A->nC, beta,
         A->data, LAH_LEADING_DIM(A), B->data, LAH_LEADING_DIM(B), alpha, C->data, LAH_LEADING_DIM(C));

    return lahReturnOk;
}
#else /* fallback to primitive linear algebra subroutines */
#include <stdlib.h>
/* Computes the general Matrix mulitplikation of form:
 * C = alpha*C + beta*op(A)*op(B)
 * similar to the BLAS routine ?gemm.
 * op(A) = A for transA=lahNorm or op(A) = At for transA = lahTrans */
lah_Return lah_matMul(lah_MatOp transA, lah_MatOp transB,
                      lah_value alpha, lah_value beta,
                      lah_mat *C, const lah_mat *A, const lah_mat *B)
{
    lah_index c, r, i;
    lah_value *res, *value1, *value2;
    lah_value temp;

    if (C == NULL || A == NULL || B == NULL)
    {
        return lahReturnParameterError;
    }

    /* NOTE: lah_matMul calls itself with an transposed view of
     * the corresponding kop_mat and sets transA resp. transB
     * to lahNorm (else we would have an endless recursion)
    */
    if (transA == lahTrans)
    {
        /* Call lah_matMul again with transposed */
        lah_mat *At = lah_matTrans(A);
        lah_Return result = lah_matMul(lahNorm, transB, alpha,
                                         beta, C, At, B);
        free(At);
        return result;
    }
    if (transB == lahTrans)
    {
        /* Call lah_matMul again with transposed of B */
        lah_mat *Bt = lah_matTrans(B);
        lah_Return result = lah_matMul(transA, lahNorm, alpha,
                                       beta, C, A, Bt);
        free(Bt);
        return result;
    }

    if (C->nR != A->nR || C->nC != B->nC || A->nC != B->nR)
    {
        return lahReturnParameterError;
    }

    for (c = 0; c < C->nC; c++)
    {
        for (r = 0; r < C->nR; r++)
        {
            res = C->data + C->incRow * r + c * C->incCol;
            value1 = A->data + r * A->incRow;
            value2 = B->data + c * B->incCol;
            temp = 0;
            for (i = 0; i < A->nC; i++, value1 += A->incCol, value2 += B->incRow)
            {
                temp += *value1 * *value2;
            }
            *res = alpha * *res + beta * temp;
        }
    }

    return lahReturnOk;
}
#endif /* endif primitive implementation */


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "../ellipsoid/ellipsoid.h"
#include "../ellipsoid/ellSolv.h"
#include "../sphere/sphere.h"
#include "../sphere/sphSolv.h"
#include "../constants.h"
#include "testFunctions.h"

//compares ellipsoidal and spherical
//has ellipsoidal distribution of charges just inside ellipsoidal boundary
//computes the solution at two point, one inside smallest enclosing sphere where spherical diverges
//other solution point is outside of sphere and ellipse to compare convergence
#undef __FUNCT__
#define __FUNCT__ "main"
PetscErrorCode main(int argc, char **argv)
{
  PetscErrorCode ierr;
  //semi-axes of ellipsoid
  PetscReal a = 3.0;
  PetscReal b = 2.5;
  PetscReal c = 1.2;
  //sphere radius
  PetscReal sphRad = 3;
  //dielecric permitivities both set to 1
  const PetscReal eps1 = 1.0;
  const PetscReal eps2 = 1.0;
  //charges placed on ellipsoid with sem-axes smaller by epsilon
  PetscReal epsilon = .5;
  //interior charge XYZ and magnitude
  const PetscInt SQRT_NUM_CHARGES = 3;
  Vec chargeXYZ;
  Vec chargeMag;
  //maximum expansion order for test
  const PetscInt MAX_N = 25;
  //solutions vector
  const PetscInt NUM_SOL_PTS = 2;
  Vec ellSolutions[MAX_N], sphSolutions[MAX_N];
  Vec ellError[MAX_N], sphError[MAX_N];
  Vec exactSolution;
  Vec solXYZ;
  //solution point inside of sphere
  PetscReal sol1X = .1;
  PetscReal sol1Y = .1;
  PetscReal sol1Z = 1.5;
  //solution point outside of sphere
  PetscReal sol2X = 2.1;
  PetscReal sol2Y = 2.6;
  PetscReal sol2Z = 1;
  //pointer for accessing petsc vectors
  PetscScalar *vecPt, *vecPt2;
  PetscFunctionBeginUser;

  //initialize petsc
  ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);
  ierr = PetscLogDefaultBegin();CHKERRQ(ierr);

  //create, set charge magnitudes to 1
  ierr = VecCreateSeq(PETSC_COMM_SELF, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, &chargeMag);CHKERRQ(ierr);
  ierr = VecSet(chargeMag, 1.0);CHKERRQ(ierr);

  //create, set charge XYZ values to be on a-eps,b-eps,c-eps ellipsoid
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, &chargeXYZ);CHKERRQ(ierr);
  ierr = VecGetArray(chargeXYZ, &vecPt);CHKERRQ(ierr);
  PetscReal pu, pv;
  for(PetscInt i=0; i<SQRT_NUM_CHARGES; ++i) {
    for(PetscInt j=0; j<SQRT_NUM_CHARGES; ++j) {
      pu = i*2.0*PETSC_PI/(SQRT_NUM_CHARGES-1);
      pv = j*PETSC_PI/(SQRT_NUM_CHARGES-1);
      PetscInt ind = i*SQRT_NUM_CHARGES + j;
      //there is a slight issue with on-axes coordinate transforms which is what the +.01 is for
      vecPt[3*ind + 0] = (a-epsilon)*PetscCosReal(pu)*PetscSinReal(pv) +.01;
      vecPt[3*ind + 1] = (b-epsilon)*PetscSinReal(pu)*PetscSinReal(pv) + .01;
      vecPt[3*ind + 2] = (c-epsilon)*PetscCosReal(pv) + .01;
      printf("pu: %3.3f pv: %3.3f\n", pu, pv);
    }
  }
  ierr = VecRestoreArray(chargeXYZ, &vecPt);CHKERRQ(ierr);

  //set solution points vector
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*NUM_SOL_PTS, &solXYZ);CHKERRQ(ierr);
  ierr = VecGetArray(solXYZ, &vecPt);CHKERRQ(ierr);
  vecPt[0] = sol1X; vecPt[1] = sol1Y; vecPt[2] = sol1Z;
  vecPt[3] = sol2X; vecPt[4] = sol2Y; vecPt[5] = sol2Z;
  ierr = VecRestoreArray(solXYZ, &vecPt);CHKERRQ(ierr);

  //caulculate exact solution
  ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, &exactSolution);CHKERRQ(ierr);
  ierr = CoulombExact(eps1, chargeXYZ, chargeMag, solXYZ, exactSolution);CHKERRQ(ierr);

  //initialize solutions and error vectors
  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, ellSolutions+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, sphSolutions+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, ellError+k);CHKERRQ(ierr);
    ierr = VecCreateSeq(PETSC_COMM_SELF, NUM_SOL_PTS, sphError+k);CHKERRQ(ierr);
  }

  //compute solutions for expansions up to order MAX_N
  for(PetscInt k=1; k<=MAX_N; ++k) {
    printf("N = %d\n", k);
    printf("calculating ellipsoidal solution...\n");
    ierr = CalcEllipsoidCoulombPotential(a, b, c, eps1, eps2, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, chargeXYZ, chargeMag, NUM_SOL_PTS, solXYZ, k, ellSolutions[k-1]);CHKERRQ(ierr);
    printf("calculating spherical solution...\n");
    ierr = CalcSphericalCoulombPotential(1.0, eps1, eps2, SQRT_NUM_CHARGES*SQRT_NUM_CHARGES, chargeXYZ, chargeMag, NUM_SOL_PTS, solXYZ, k, sphSolutions[k-1]);CHKERRQ(ierr);
  }

  //compute errors
  printf("Computing errors...\n");
  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecCopy(exactSolution, ellError[k]);CHKERRQ(ierr);
    ierr = VecCopy(exactSolution, sphError[k]);CHKERRQ(ierr);
    ierr = VecAXPY(ellError[k], -1.0, ellSolutions[k]);CHKERRQ(ierr);
    ierr = VecAXPY(sphError[k], -1.0, sphSolutions[k]);CHKERRQ(ierr);
    ierr = VecAbs(ellError[k]);CHKERRQ(ierr);
    ierr = VecAbs(sphError[k]);CHKERRQ(ierr);
  }

  FILE *fp = fopen("out/ex1.txt", "w");

  for(PetscInt k=0; k<MAX_N; ++k) {
    ierr = VecGetArrayRead(ellError[k], &vecPt);
    ierr = VecGetArrayRead(sphError[k], &vecPt2);
    printf("Error[%d] = %4.4e\n", k, vecPt[0]);
    fprintf(fp, "%4.4e %4.4e %4.4e %4.4e\n", vecPt[0], vecPt2[0], vecPt[1], vecPt2[1]);
    ierr = VecRestoreArrayRead(ellError[k], &vecPt);
    ierr = VecRestoreArrayRead(sphError[k], &vecPt2);
  }
  printf("Exact: \n");
  ierr = VecView(exactSolution, PETSC_VIEWER_STDOUT_SELF);
  printf("Ell: \n");
  ierr = VecView(ellSolutions[MAX_N-1], PETSC_VIEWER_STDOUT_SELF);

  ierr = PetscFinalize();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}


text:
 /**
 *
 * @file core_dormqr.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:46 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dormqr overwrites the general complex M-by-N tile C with
 *
 *                    SIDE = 'L'     SIDE = 'R'
 *    TRANS = 'N':      Q * C          C * Q
 *    TRANS = 'C':      Q**T * C       C * Q**T
 *
 *  where Q is a complex unitary matrix defined as the product of k
 *  elementary reflectors
 *
 *    Q = H(1) H(2) . . . H(k)
 *
 *  as returned by CORE_dgeqrt. Q is of order M if SIDE = 'L' and of order N
 *  if SIDE = 'R'.
 *
 *******************************************************************************
 *
 * @param[in] side
 *         @arg PlasmaLeft  : apply Q or Q**T from the Left;
 *         @arg PlasmaRight : apply Q or Q**T from the Right.
 *
 * @param[in] trans
 *         @arg PlasmaNoTrans   :  No transpose, apply Q;
 *         @arg PlasmaTrans :  Transpose, apply Q**T.
 *
 * @param[in] M
 *         The number of rows of the tile C.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile C.  N >= 0.
 *
 * @param[in] K
 *         The number of elementary reflectors whose product defines
 *         the matrix Q.
 *         If SIDE = PlasmaLeft,  M >= K >= 0;
 *         if SIDE = PlasmaRight, N >= K >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in] A
 *         Dimension:  (LDA,K)
 *         The i-th column must contain the vector which defines the
 *         elementary reflector H(i), for i = 1,2,...,k, as returned by
 *         CORE_dgeqrt in the first k columns of its array argument A.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.
 *         If SIDE = PlasmaLeft,  LDA >= max(1,M);
 *         if SIDE = PlasmaRight, LDA >= max(1,N).
 *
 * @param[in] T
 *         The IB-by-K triangular factor T of the block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] LDT
 *         The leading dimension of the array T. LDT >= IB.
 *
 * @param[in,out] C
 *         On entry, the M-by-N tile C.
 *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 *
 * @param[in] LDC
 *         The leading dimension of the array C. LDC >= max(1,M).
 *
 * @param[in,out] WORK
 *         On exit, if INFO = 0, WORK(1) returns the optimal LDWORK.
 *
 * @param[in] LDWORK
 *         The dimension of the array WORK.
 *         If SIDE = PlasmaLeft,  LDWORK >= max(1,N);
 *         if SIDE = PlasmaRight, LDWORK >= max(1,M).
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dormqr = PCORE_dormqr
#define CORE_dormqr PCORE_dormqr
#endif
int CORE_dormqr(PLASMA_enum side, PLASMA_enum trans,
                int M, int N, int K, int IB,
                const double *A, int LDA,
                const double *T, int LDT,
                double *C, int LDC,
                double *WORK, int LDWORK)
{
    int i, kb;
    int i1, i3;
    int nq, nw;
    int ic = 0;
    int jc = 0;
    int ni = N;
    int mi = M;

    /* Check input arguments */
    if ((side != PlasmaLeft) && (side != PlasmaRight)) {
        coreblas_error(1, "Illegal value of side");
        return -1;
    }
    /*
     * NQ is the order of Q and NW is the minimum dimension of WORK
     */
    if (side == PlasmaLeft) {
        nq = M;
        nw = N;
    }
    else {
        nq = N;
        nw = M;
    }

    if ((trans != PlasmaNoTrans) && (trans != PlasmaTrans)) {
        coreblas_error(2, "Illegal value of trans");
        return -2;
    }
    if (M < 0) {
        coreblas_error(3, "Illegal value of M");
        return -3;
    }
    if (N < 0) {
        coreblas_error(4, "Illegal value of N");
        return -4;
    }
    if ((K < 0) || (K > nq)) {
        coreblas_error(5, "Illegal value of K");
        return -5;
    }
    if ((IB < 0) || ( (IB == 0) && ((M > 0) && (N > 0)) )) {
        coreblas_error(6, "Illegal value of IB");
        return -6;
    }
    if ((LDA < max(1,nq)) && (nq > 0)) {
        coreblas_error(8, "Illegal value of LDA");
        return -8;
    }
    if ((LDC < max(1,M)) && (M > 0)) {
        coreblas_error(12, "Illegal value of LDC");
        return -12;
    }
    if ((LDWORK < max(1,nw)) && (nw > 0)) {
        coreblas_error(14, "Illegal value of LDWORK");
        return -14;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (K == 0))
        return PLASMA_SUCCESS;

    if (((side == PlasmaLeft) && (trans != PlasmaNoTrans))
        || ((side == PlasmaRight) && (trans == PlasmaNoTrans))) {
        i1 = 0;
        i3 = IB;
    }
    else {
        i1 = ( ( K-1 ) / IB )*IB;
        i3 = -IB;
    }

    for(i = i1; (i >- 1) && (i < K); i+=i3 ) {
        kb = min(IB, K-i);

        if (side == PlasmaLeft) {
            /*
             * H or H' is applied to C(i:m,1:n)
             */
            mi = M - i;
            ic = i;
        }
        else {
            /*
             * H or H' is applied to C(1:m,i:n)
             */
            ni = N - i;
            jc = i;
        }
        /*
         * Apply H or H'
         */
        LAPACKE_dlarfb_work(LAPACK_COL_MAJOR,
            lapack_const(side),
            lapack_const(trans),
            lapack_const(PlasmaForward),
            lapack_const(PlasmaColumnwise),
            mi, ni, kb,
            &A[LDA*i+i], LDA,
            &T[LDT*i], LDT,
            &C[LDC*jc+ic], LDC,
            WORK, LDWORK);
    }
    return PLASMA_SUCCESS;
}


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

#include <config.h>
#include <gsl/gsl_multiroots.h>

int
gsl_multiroot_fdjacobian (gsl_multiroot_function * F,
                           const gsl_vector * x, const gsl_vector * f,
                           double epsrel, gsl_matrix * jacobian)
{
  const size_t n = x->size;
  const size_t m = f->size;
  const size_t n1 = jacobian->size1;
  const size_t n2 = jacobian->size2;

  if (m != n1 || n != n2)
    {
      GSL_ERROR ("function and jacobian are not conformant", GSL_EBADLEN);
    }

  {
    size_t i,j;
    gsl_vector *x1, *f1;

    x1 = gsl_vector_alloc (n);

    if (x1 == 0)
      {
	GSL_ERROR ("failed to allocate space for x1 workspace", GSL_ENOMEM);
      }

    f1 = gsl_vector_alloc (m);

    if (f1 == 0)
      {
	gsl_vector_free (x1);

	GSL_ERROR ("failed to allocate space for f1 workspace", GSL_ENOMEM);
      }

    gsl_vector_memcpy (x1, x);	/* copy x into x1 */

    for (j = 0; j < n; j++)
      {
	double xj = gsl_vector_get (x, j);
	double dx = epsrel * fabs (xj);

	if (dx == 0)
	  {
	    dx = epsrel;
	  }

	gsl_vector_set (x1, j, xj + dx);

        {
          int status = GSL_MULTIROOT_FN_EVAL (F, x1, f1);

          if (status != GSL_SUCCESS)
            {
              return GSL_EBADFUNC;
            }
        }

	gsl_vector_set (x1, j, xj);

	for (i = 0; i < m; i++)
	  {
	    double g1 = gsl_vector_get (f1, i);
	    double g0 = gsl_vector_get (f, i);
	    gsl_matrix_set (jacobian, i, j, (g1 - g0) / dx);
	  }
      }

    gsl_vector_free (x1);
    gsl_vector_free (f1);
  }

  return GSL_SUCCESS;
}


text:
//        Copyright The Authors 2018.
//    Distributed under the 3-Clause BSD License.
//    (See accompanying file LICENSE or copy at
//   https://opensource.org/licenses/BSD-3-Clause)

#pragma once

#include <array>     // for std::array
#include <bitset>    // for std::bitset
#include <cstdint>   // for standard int types
#include <cstring>   // for std::memcpy
#include <iosfwd>    // for implementation of operator<<
#include <iterator>  // for std::reverse_iterator
#include <string>    // for std::string

#include <gsl/gsl>  // for gsl::span

#include <codec/base16.h>   // for hex enc/decoding
#include <common/assert.h>  // for assertions
#include <crypto/random.h>  // for random block generation

namespace blocxxi {
namespace crypto {

namespace detail {
/// Forward declarations for utility conversion function between host and
/// network (big endian) byte orders.
/// Avoid including header files of boost and pulling all its bloat here.
std::uint32_t HostToNetwork(std::uint32_t);
std::uint32_t NetworkToHost(std::uint32_t);

/// Count the number of leading zero bits in a buffer of contiguous 32-bit
/// integers.
int CountLeadingZeroBits(gsl::span<std::uint32_t const> buf);

}  // namespace detail

/*!
   Represents a N bits hash digest or any other kind of N bits sequence.
   The structure is 32-bit aligned and must be at least 32-bits wide. It
   has almost the same semantics of std::array<> except that it will never be
   zero-length (size() is always > 0). It offers a number of additional
   convenience methods, such as bitwise arithemtics, comparisons etc.
 */
template <unsigned int BITS>
class Hash {
  static_assert(BITS % 32 == 0, "Hash size in bits must be a multiple of 32");
  static_assert(BITS > 0, "Hash size in bits must be greater than 0");

 public:
  using value_type = std::uint8_t;
  using pointer = uint8_t *;
  using const_pointer = std::uint8_t const *;
  using size_type = std::size_t;
  using difference_type = std::ptrdiff_t;
  using reference = std::uint8_t &;
  using const_reference = std::uint8_t const &;
  using iterator = pointer;
  using const_iterator = const_pointer;
  using reverse_iterator = std::reverse_iterator<iterator>;
  using const_reverse_iterator = std::reverse_iterator<const_iterator>;

  ///@name Constructors and Factory methods
  //@{
  /// Initializes the hash with all '0' bits
  Hash() : storage_({}) {}

  Hash(const Hash &other) = default;

  Hash &operator=(const Hash &rhs) = default;

  Hash(Hash &&) noexcept = default;

  Hash &operator=(Hash &&rhs) noexcept = default;

  virtual ~Hash() = default;

  /*!
  \brief Construct a Hash by assigning content from the provided sequence of
  bytes.

  > The source sequence can be smaller than the hash size but it __cannot be
  > bigger__.

  In the case it is smaller, the hash will be padded (leading bytes)
  with 0. The last byte in the source sequence will be the LSB of the hash
  (last byte in the hash).

  In the case it is bigger, result is undefined and, if assertions are enabled,
  this would be considered as a fatal error.

  \param buf source sequence of bytes.
  */
  explicit Hash(gsl::span<std::uint8_t const> buf) noexcept : storage_({}) {
    auto src_size = buf.size();
    ASAP_ASSERT_PRECOND(src_size <= Size());
    auto start = begin();
    if (Size() > src_size) {
      // Pad with zeros and adjust the start
      std::memset(Data(), 0, Size() - src_size);
      start += Size() - src_size;
    }
    Assign(buf, start);
  }

  /// \brief Returns an all-1 hash, i.e. the biggest value representable by an N
  /// bit number (N/8 bytes).
  static Hash Max() noexcept {
    Hash h;
    h.storage_.fill(0xffffffffU);
    return h;
  }

  /// \brief Returns an all-0 hash, i.e. the smallest value representable by an
  /// N bit number (N/8 bytes).
  static Hash Min() noexcept {
    return Hash();
    // all bits are already 0
  }

  static Hash<BITS> FromHex(const std::string &src, bool reverse = false) {
    Hash<BITS> hash;
    codec::hex::Decode(gsl::make_span<>(src),
                       gsl::make_span<>(hash.Data(), hash.Size()), reverse);
    return hash;
  }

  static Hash<BITS> RandomHash() {
    Hash<BITS> hash;
    hash.Randomize();
    return hash;
  }
  //@}

  /// @name Checked element access
  /// access specified element with bounds checking.
  //@{
  /*!
    Returns a reference to the element at specified location `pos`, with
    bounds checking.

    \param pos position of the element to return
    \exception std::out_of_range if !(pos < size()).
   */
  reference At(size_type pos) {
    if (Size() <= pos) {
      throw std::out_of_range("Hash<BITS> index out of range");
    }
    return reinterpret_cast<pointer>(storage_.data())[pos];
  }
  const_reference At(size_type pos) const {
    if (Size() <= pos) {
      throw std::out_of_range("Hash<BITS> index out of range");
    }
    return reinterpret_cast<const_pointer>(storage_.data())[pos];
  };
  //@}

  /// @name Unchecked element access
  /// access specified element __without__ bounds checking.
  //@{
  /// Returns a reference to the element at specified location pos. No bounds
  /// checking is performed.
  reference operator[](size_type pos) noexcept {
    ASAP_ASSERT_PRECOND(pos < Size());
    return reinterpret_cast<pointer>(storage_.data())[pos];
  }
  const_reference operator[](size_type pos) const noexcept {
    ASAP_ASSERT_PRECOND(pos < Size());
    return reinterpret_cast<pointer>(storage_.data())[pos];
  }

  /// Returns a reference to the first element in the container.
  reference Front() { return reinterpret_cast<pointer>(storage_.data())[0]; }
  const_reference Front() const {
    return reinterpret_cast<const_pointer>(storage_.data())[0];
  }

  /// Returns reference to the last element in the container.
  reference Back() {
    return reinterpret_cast<pointer>(storage_.data())[Size() - 1];
  }
  const_reference Back() const {
    return reinterpret_cast<const_pointer>(storage_.data())[Size() - 1];
  }

  /// Returns pointer to the underlying array serving as element storage. The
  /// pointer is such that range [Data(); Data() + Size()) is always a valid
  /// range.
  pointer Data() noexcept { return reinterpret_cast<pointer>(storage_.data()); }
  const_pointer Data() const noexcept {
    return reinterpret_cast<const_pointer>(storage_.data());
  }
  //@}

  /// @name Iterators
  //@{
  /// Returns an iterator to the first element of the container.
  iterator begin() noexcept {
    return reinterpret_cast<pointer>(storage_.data());
  }
  const_iterator begin() const noexcept {
    return reinterpret_cast<const_pointer>(storage_.data());
  }
  const_iterator cbegin() const noexcept { return begin(); }

  /// Returns an iterator to the element following the last element of the
  /// container. This element acts as a placeholder; attempting to access it
  /// results in undefined behavior.
  iterator end() noexcept {
    return reinterpret_cast<pointer>(storage_.data()) + Size();
  }
  const_iterator end() const noexcept {
    return reinterpret_cast<const_pointer>(storage_.data()) + Size();
  }
  const_iterator cend() const noexcept { return end(); }

  /// Returns a reverse iterator to the first element of the reversed container.
  /// It corresponds to the last element of the non-reversed container.
  reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
  const_reverse_iterator rbegin() const noexcept {
    return const_reverse_iterator(end());
  }
  const_reverse_iterator crbegin() const noexcept { return rbegin(); }

  /// Returns a reverse iterator to the element following the last element of
  /// the reversed container. It corresponds to the element preceding the first
  /// element of the non-reversed container. This element acts as a placeholder,
  /// attempting to access it results in undefined behavior.
  reverse_iterator rend() noexcept { return reverse_iterator(begin()); }
  const_reverse_iterator rend() const noexcept {
    return const_reverse_iterator(begin());
  }
  const_reverse_iterator crend() const noexcept { return rend(); }
  //@}

  // Capacity

  /// The size of the hash in bytes.
  constexpr static std::size_t Size() { return BITS / 8; };

  /// The size of the hash in bytes.
  constexpr static std::size_t BitSize() { return BITS; };

  // Operations

  /// Set all bits to 0.
  void Clear() { storage_.fill(0); }

  /// Return true if all bits are set to 0.
  bool IsAllZero() const {
    for (auto v : storage_)
      if (v != 0) return false;
    return true;
  }

  /// Exchanges the contents of the container with those of other. Does not
  /// cause iterators and references to associate with the other container.
  void swap(Hash &other) noexcept { storage_.swap(other.storage_); }

  /// Count leading zero bits.
  int LeadingZeroBits() const {
    return detail::CountLeadingZeroBits(gsl::make_span(storage_));
  }

  /// @name Comparison operators
  /// The important thing to note here is that only two of these operators
  /// actually do anything, the others are just forwarding their arguments to
  /// either of these two to do the actual work (== and <). Therefore only these
  /// two are implemented as friends.
  //@{
  /// \brief Checks if the contents of `lhs` and `rhs` are equal, that is,
  /// whether each element in `lhs` compares equal with the element in `rhs` at
  /// the same position.
  friend bool operator==(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
    return std::equal(lhs.begin(), lhs.end(), rhs.begin());
  }
  /// \brief Compares the contents of `lhs` and `rhs` lexicographically. The
  /// comparison is performed by a function equivalent to
  /// `std::lexicographical_compare`.
  friend bool operator<(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
    for (std::size_t i = 0; i < SIZE_DWORD; ++i) {
      std::uint32_t const native_lhs = detail::NetworkToHost(lhs.storage_[i]);
      std::uint32_t const native_rhs = detail::NetworkToHost(rhs.storage_[i]);
      if (native_lhs < native_rhs) return true;
      if (native_lhs > native_rhs) return false;
    }
    return false;
  }
  //@}

  /// In-place bitwise XOR with the given other hash.
  Hash &operator^=(Hash const &other) noexcept {
    for (std::size_t i = 0; i < SIZE_DWORD; ++i)
      storage_[i] ^= other.storage_[i];
    return *this;
  }

  /*!
  Starting at the given position in the hash, assign contents from the given
  sequence of bytes.

  > Data represented by the span is expected to be organized in network byte
  > order (i.e. the byte returned by span.first() is the MSB (Most Significant
  > Byte) and therefore will go into the MSB of the Hash).

  The _span_ size can be smaller than the remaining space after _start_ in
  the hash but it __cannot be bigger__. If this is not true, result is undefined
  and, if assertions are enabled, this would be considered as a fatal error.

  \param buf source sequence of bytes.
  \param start position at which to assign the sequence in the hash.
  */
  void Assign(gsl::span<std::uint8_t const> buf, iterator start) noexcept {
    size_type src_size = buf.size();
    size_type dst_size = end() - start;
    ASAP_ASSERT_PRECOND(src_size <= dst_size);
    std::memcpy(start, buf.data(), std::min(src_size, dst_size));
  }

  // TODO: Replace the parameters with a span
  void Randomize() { random::GenerateBlock(Data(), Size()); }

  const std::string ToHex() const {
    return codec::hex::Encode(gsl::make_span<>(Data(), Size()), false, true);
  }

  std::bitset<BITS> ToBitSet() const {
    std::bitset<BITS> bs;  // [0,0,...,0]
    int shift_left = BITS;
    for (uint32_t num_part : storage_) {
      shift_left -= 32;
      num_part = detail::HostToNetwork(num_part);
      std::bitset<BITS> bs_part(num_part);
      bs |= bs_part << shift_left;
    }
    return bs;
  }

  std::string ToBitStringShort(std::size_t length = 32U) const {
    auto truncate = false;
    if (length < BITS) { truncate = true; length -= 3; }
    auto str = ToBitSet().to_string().substr(0, length);
    if (truncate) str.append("...");
    return str;
  }

 private:
  /// The size of the underlying data buffer in terms of 32-bit DWORD
  static constexpr std::size_t SIZE_DWORD = BITS / 32;

  /// The underlying data buffer is built as an array of unsigned 32-bit
  /// integers to facilitate the math operations. Note however that iteration
  /// over the Hash and random access to the data should always be presented as
  /// if it were a byte array layed out in network order (big endian).
  ///
  /// > First DWORD is the most significant 4 bytes.
  std::array<std::uint32_t, SIZE_DWORD> storage_;
};

/// @name Comparison operators
/// Additional Comparison operators (not needing to be friend)
//@{

template <unsigned int BITS>
inline bool operator!=(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
  return !(lhs == rhs);
}
template <unsigned int BITS>
inline bool operator>(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
  return rhs < lhs;
}
template <unsigned int BITS>
inline bool operator<=(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
  return !(lhs > rhs);
}
template <unsigned int BITS>
inline bool operator>=(const Hash<BITS> &lhs, const Hash<BITS> &rhs) {
  return !(lhs < rhs);
}
//@}

/// Bitwise XOR
template <unsigned int BITS>
inline Hash<BITS> operator^(Hash<BITS> lhs, Hash<BITS> const &rhs) noexcept {
  return lhs.operator^=(rhs);
}

// Specialized algorithms

/// Exchanges the contents of one with other. Does not cause
/// iterators and references to associate with the other container.
template <unsigned int BITS>
inline void swap(Hash<BITS> &lhs, Hash<BITS> &rhs) {
  lhs.swap(rhs);
}

/// Dump the hex-encoded contents of the hash to the stream.
template <unsigned int BITS>
inline std::ostream &operator<<(std::ostream &out, Hash<BITS> const &hash) {
  out << hash.ToHex();
  return out;
}

using Hash512 = Hash<512>;  // 64 bytes
using Hash256 = Hash<256>;  // 32 bytes
using Hash160 = Hash<160>;  // 20 bytes

}  // namespace crypto
}  // namespace blocxxi