created a git repository for my thesis code.

17 files changed, 4726 insertions, 0 deletions

diff --git a/src/bifurchaser.cpp b/src/bifurchaser.cpp
new file mode 100644
index 0000000..1ef595b
--- /dev/null
+++ b/src/bifurchaser.cpp
@@ -0,0 +1,359 @@
+/* bifurcation_chaser.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* A program which facilitates automated mapping of bifurcation points in the
+ * energy of a system where the Hessian is available.  Currently, only a one
+ * dimensional parameter space is supported.
+ */
+
+#include "domain_energy.h"
+#include "domain_minimize.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+void geteigenvalues(gsl_vector *eigenvalues, unsigned n, const gsl_vector *z, double c) {
+	gsl_matrix *hess;
+	hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3);
+
+	domain_energy_fixedHessian(hess, n, z, c);
+
+	domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess);
+	gsl_matrix_free(hess);
+}
+
+
+int domain_eigen_perturb(gsl_vector *z, unsigned k, unsigned n,
+		unsigned eigen_num, double a0,double a_fact, double c, double eps,
+		double g, double N, double energy_thres) {
+
+	printf("Beginning perturbation.\n");
+
+	double a, temp_eigenval, eigenval;
+	unsigned kk;
+
+	gsl_vector *temp_z, *eigenvalues, *eigenvector;
+	gsl_permutation *eigenorder;
+	gsl_matrix *hess;
+
+	eigenvalues = gsl_vector_alloc(3 * n + 3);
+	eigenvector = gsl_vector_alloc(3 * n + 3);
+	temp_z = gsl_vector_alloc(3 * n + 3);
+	hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3);
+
+	eigenorder = gsl_permutation_alloc(3 * n + 3);
+
+	domain_energy_fixedHessian(hess, n, z, c);
+
+	domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess);
+	domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues);
+
+	kk = gsl_permutation_get(eigenorder, k);
+
+	eigenval = gsl_vector_get(eigenvalues, kk);
+
+	printf("Getting eigenvector.\n");
+	domain_energy_fixedHessian(hess, n, z, c);
+	domain_eigen_vector(eigenvector, 3 * n + 3, 2 * n, kk, hess);
+
+	a = a0;
+	int failed = 0;
+
+	printf("Starting loop.\n");
+	while (true) {
+		gsl_vector_memcpy(temp_z, z);
+		gsl_blas_daxpy(a, eigenvector, temp_z);
+		failed = domain_minimize_fixed(z, n, c, eps, N, 0.9, g, 0.9);
+
+		if (failed) {
+			printf("Relaxation failed, reducing perturb size.\n");
+			a *= 0.1;
+		} else {
+
+	domain_energy_fixedHessian(hess, n, z, c);
+
+	domain_eigen_values(eigenvalues, 3 * n + 3, 2 * n, hess);
+		domain_eigen_sort(eigenorder,  3 * n + 3, eigen_num, eigenvalues);
+
+		kk = gsl_permutation_get(eigenorder, k);
+
+		temp_eigenval = gsl_vector_get(eigenvalues, kk);
+
+		printf("BIFUR: Perturbing %i, %e, %e\n", k, eigenval, temp_eigenval);
+
+		if (GSL_SIGN(temp_eigenval) != GSL_SIGN(eigenval)) {
+			gsl_vector_memcpy(z, temp_z);
+			break;
+		}
+
+		a *= a_fact;
+		}
+	}
+
+	gsl_vector_free(eigenvector);
+	gsl_vector_free(temp_z);
+	gsl_vector_free(eigenvalues);
+
+	gsl_permutation_free(eigenorder);
+
+	return 0;
+}
+
+
+bool bifur_consent() {
+	printf(" (y/n): ");
+	char in;
+	in = getchar();
+	getchar();
+	if (in == 'y') return true;
+	else return false;
+}
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt, min_fails, eigen_follow, eigen_num, examining;
+	unsigned n, N, j, a, last_pert, ii, old_ii;
+	double c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, state, old_state;
+	char *filename, str[19], in;
+	bool subcrit, reset;
+
+	// Setting default values.
+	eps = 0;
+	eigen_thres = 1e-13;
+	approach_thres = 1e-6;
+	eigen_follow = -1;
+	examining = -1;
+	eigen_num = 25;
+	last_pert = 0;
+	subcrit = false;
+	reset = false;
+	dc = 0;
+
+	j = 0;
+
+	gsl_vector *z, *old_z, *eigenvalues, *eigenstate, *old_eigenstate, *eigenchanges;
+	gsl_permutation *eigenorder, *old_eigenorder;
+
+	while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'N':
+				N = atoi(optarg);
+				break;
+			case 'j':
+				j = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'd':
+				dc0 = atof(optarg);
+				break;
+			case 'h':
+				dc = atof(optarg);
+				break;
+			case 'g':
+				g0 = atof(optarg);
+				break;
+			case 'i':
+				filename = optarg;
+				break;
+			case 'm':
+				eigen_follow = atof(optarg);
+				break;
+			case 'e':
+				eps = atof(optarg);
+				break;
+			case 's':
+				subcrit = true;
+				break;
+			case 't':
+				approach_thres = atof(optarg);
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	z = gsl_vector_alloc(3 * n + 3);
+	old_z = gsl_vector_alloc(3 * n + 3);
+	eigenvalues = gsl_vector_alloc(3 * n + 3);
+	eigenorder = gsl_permutation_alloc(3 * n + 3);
+	old_eigenorder = gsl_permutation_alloc(3 * n + 3);
+	old_eigenstate = gsl_vector_alloc(3 * n + 3);
+	eigenstate = gsl_vector_alloc(3 * n + 3);
+
+	FILE *f = fopen(filename, "r");
+	gsl_vector_fscanf(f, z);
+	fclose(f);
+
+	g = g0;
+	if (dc == 0) dc = dc0;
+
+	min_fails = domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9);
+
+	if (min_fails) {
+		printf("BIFUR: Initial relaxation failed, exiting.\n");
+		return 1;
+	}
+
+	geteigenvalues(eigenvalues, n, z, c);
+	domain_eigen_state(old_eigenstate, eigenvalues, n, eigen_thres);
+	domain_eigen_sort(old_eigenorder, 3 * n + 3, eigen_num, eigenvalues);
+
+	while (true) {
+		j += 1;
+		c += dc;
+		reset = false;
+
+		gsl_vector_memcpy(old_z, z);
+
+		printf("BIFUR: Step %05d, starting with c = %f.\n", j, c);
+
+		min_fails = domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9);
+
+		if (min_fails) {
+			printf("BIFUR: Newton's method failed to converge, reducing step size.\n");
+			c -= dc;
+			j -= 1;
+			last_pert = 0;
+			gsl_vector_memcpy(z, old_z);
+			dc *= 0.1;
+			reset = true;
+		} else {
+
+		geteigenvalues(eigenvalues, n, z, c);
+		domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num, eigenvalues);
+		domain_eigen_state(eigenstate, eigenvalues, n, eigen_thres);
+
+		if (eigen_follow > -1) examining = eigen_follow;
+
+		for (unsigned i = 0; i < eigen_num; i++) {
+			ii = gsl_permutation_get(eigenorder, i);
+			old_ii = gsl_permutation_get(old_eigenorder, i);
+
+			state = gsl_vector_get(eigenstate, ii);
+			old_state = gsl_vector_get(old_eigenstate, old_ii);
+
+			if (state != old_state) {
+				if (i == examining) {
+					c -= dc;
+					gsl_vector_memcpy(z, old_z);
+					gsl_vector_memcpy(eigenstate, old_eigenstate);
+					gsl_permutation_memcpy(eigenorder, old_eigenorder);
+					j -= 1;
+					dc *= 0.1;
+					reset = true;
+					last_pert = 0;
+				} if (examining == -1 && state != 0 && old_state != 0) {
+					printf("BIFUR: Eigenvalue %i changed sign past threshold to %e.  Examine?", i,
+						gsl_vector_get(eigenvalues, ii));
+					if (bifur_consent()) {
+
+						examining = i;
+						c -= dc;
+						gsl_vector_memcpy(z, old_z);
+						gsl_vector_memcpy(eigenstate, old_eigenstate);
+						gsl_permutation_memcpy(eigenorder, old_eigenorder);
+						j -= 1;
+						dc *= 0.1;
+						reset = true;
+						last_pert = 0;
+
+						break;
+					}
+				}
+			}
+		}
+
+		if (!reset && examining > -1 && fabs(dc) < approach_thres) {
+
+			if (!subcrit) {
+				c += GSL_SIGN(dc) * approach_thres;
+				domain_minimize_fixed(z, n, c, eps, N, 0.9, 1, 0.9);
+			}
+
+			printf("BIFUR: Perturbing at c = %.8f.\n", c);
+
+			domain_eigen_perturb(z, examining, n, eigen_num, 1, 1.1, c, eps, g, N, 0);
+			geteigenvalues(eigenvalues, n, z, c);
+			domain_eigen_sort(eigenorder, 3 * n + 3, eigen_num,  eigenvalues);
+			domain_eigen_state(eigenstate, eigenvalues, n, eigen_thres);
+
+			if (subcrit) dc = - GSL_SIGN(dc) * approach_thres;
+			else dc = GSL_SIGN(dc) * approach_thres;
+
+			examining = -1;
+			last_pert = 0;
+		}
+
+		if (!reset) {
+
+			gsl_vector_memcpy(old_eigenstate, eigenstate);
+			gsl_permutation_memcpy(old_eigenorder, eigenorder);
+
+			if (last_pert > 10 && fabs(dc) < fabs(dc0)) {
+				last_pert = 0;
+				dc = GSL_SIGN(dc) * fmin(fabs(dc) * 10, fabs(dc0));
+			}
+
+			last_pert += 1;
+
+			double energy = domain_energy_fixedEnergy(n, z, c);
+
+			sprintf(str, "output/out-%05d.dat", j);
+			FILE *fout = fopen(str, "w");
+			fprintf(fout, "%.10e\t%.10e\n", c, energy);
+			for (unsigned i = 0; i < eigen_num; i++) {
+				ii = gsl_permutation_get(eigenorder, i);
+
+				fprintf(fout, "%.10e\t", gsl_vector_get(eigenvalues, ii));
+			}
+			fprintf(fout, "\n");
+			gsl_vector_fprintf(fout, z, "%.10e");
+			fclose(fout);
+		}
+	}
+	}
+
+	gsl_vector_free(z);
+
+}
+
diff --git a/src/domain_check.cpp b/src/domain_check.cpp
new file mode 100644
index 0000000..5270663
--- /dev/null
+++ b/src/domain_check.cpp
@@ -0,0 +1,274 @@
+/* domain_improve.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* A program which facilitates automated mapping of bifurcation points in the
+ * energy of a system where the Hessian is available.  Currently, only a one
+ * dimensional parameter space is supported.
+ */
+
+#include "domain_energy.h"
+#include "domain_minimize.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+void bifur_eigenvalues(gsl_vector *eigenvalues, unsigned n,
+		const gsl_vector *z, double c) {
+
+	double eigenvalue;
+
+	gsl_vector *beta;
+	gsl_vector_complex *alpha;
+	gsl_matrix *hess, *modI;
+	gsl_eigen_gen_workspace *w;
+
+	alpha = gsl_vector_complex_alloc(3 * n + 3);
+	beta = gsl_vector_alloc(3 * n + 3);
+	hess = gsl_matrix_alloc(3 * n + 3, 3 * n + 3);
+	modI = gsl_matrix_alloc(3 * n + 3, 3 * n + 3);
+	w = gsl_eigen_gen_alloc(3 * n + 3);
+
+	gsl_matrix_set_zero(modI);
+	for (unsigned i = 0; i < 2 * n; i++) gsl_matrix_set(modI, i, i, 1);
+
+	domain_energy_hessian(hess, n, z, c);
+
+	gsl_eigen_gen(hess, modI, alpha, beta, w);
+
+	for (unsigned i = 0; i < 3 * n + 3; i++) {
+		eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i);
+		gsl_vector_set(eigenvalues, i, eigenvalue);
+	}
+
+	gsl_vector_free(beta);
+	gsl_vector_complex_free(alpha);
+	gsl_matrix_free(modI);
+	gsl_matrix_free(hess);
+	gsl_eigen_gen_free(w);
+}
+
+void bifur_trueEigenvalues(gsl_vector *eigenvalues, unsigned n,
+		const gsl_vector *z, double c) {
+
+	double eigenvalue;
+
+	gsl_vector *beta;
+	gsl_vector_complex *alpha;
+	gsl_matrix *hess, *modI;
+	gsl_eigen_gen_workspace *w;
+
+	alpha = gsl_vector_complex_alloc(3 * n + 4);
+	beta = gsl_vector_alloc(3 * n + 4);
+	hess = gsl_matrix_alloc(3 * n + 4, 3 * n + 4);
+	modI = gsl_matrix_alloc(3 * n + 4, 3 * n + 4);
+	w = gsl_eigen_gen_alloc(3 * n + 4);
+
+	gsl_matrix_set_zero(modI);
+	for (unsigned i = 0; i < 2 * n + 1; i++) gsl_matrix_set(modI, i, i, 1);
+
+	domain_energy_truehessian(hess, n, z, c);
+
+	gsl_eigen_gen(hess, modI, alpha, beta, w);
+
+	for (unsigned i = 0; i < 3 * n + 4; i++) {
+		eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i);
+		gsl_vector_set(eigenvalues, i, eigenvalue);
+	}
+
+	gsl_vector_free(beta);
+	gsl_vector_complex_free(alpha);
+	gsl_matrix_free(modI);
+	gsl_matrix_free(hess);
+	gsl_eigen_gen_free(w);
+}
+
+
+void bifur_eigensort(gsl_permutation *eigenorder, unsigned n, unsigned eigen_num,
+		const gsl_vector *eigenvalues) {
+
+	unsigned ii;
+
+	gsl_vector *abs_eigenvalues;
+
+	abs_eigenvalues = gsl_vector_alloc(3 * n + 3);
+	
+	for (unsigned i = 0; i < 3 * n + 3; i++) {
+		gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i)));
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_memcpy(abs_eigenvalues, eigenvalues);
+
+	for (unsigned i = eigen_num; i < 3 * n + 3; i++) {
+		ii = gsl_permutation_get(eigenorder, i);
+
+		gsl_vector_set(abs_eigenvalues, ii, INFINITY);
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_free(abs_eigenvalues);
+}
+
+void bifur_trueEigensort(gsl_permutation *eigenorder, unsigned n, unsigned eigen_num,
+		const gsl_vector *eigenvalues) {
+
+	unsigned ii;
+
+	gsl_vector *abs_eigenvalues;
+
+	abs_eigenvalues = gsl_vector_alloc(3 * n + 4);
+	
+	for (unsigned i = 0; i < 3 * n + 4; i++) {
+		gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i)));
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_memcpy(abs_eigenvalues, eigenvalues);
+
+	for (unsigned i = eigen_num; i < 3 * n + 4; i++) {
+		ii = gsl_permutation_get(eigenorder, i);
+
+		gsl_vector_set(abs_eigenvalues, ii, INFINITY);
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_free(abs_eigenvalues);
+}
+
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt, min_fails;
+	unsigned n, N, num;
+	double c, g0, g, eps, energy;
+	char *filename;
+	bool eigenpres = true;
+
+	// Setting default values.
+	eps = 0;
+	num = 25;
+
+	gsl_vector *z, *old_z, *eigenvalues, *trueEigenvalues;
+	gsl_permutation *eigenorder, *trueEigenorder;
+
+	while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'N':
+				N = atoi(optarg);
+				break;
+			case 'g':
+				g0 = atof(optarg);
+				break;
+			case 'i':
+				filename = optarg;
+				break;
+			case 'e':
+				eps = atof(optarg);
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	z = gsl_vector_alloc(3 * n + 3);
+	old_z = gsl_vector_alloc(3 * n + 3);
+	eigenvalues = gsl_vector_alloc(3 * n + 3);
+	trueEigenvalues = gsl_vector_alloc(3 * n + 4);
+	eigenorder = gsl_permutation_alloc(3 * n + 3);
+	trueEigenorder = gsl_permutation_alloc(3 * n + 4);
+
+	g = g0;
+
+	char ch;
+	double throwaway;
+
+	FILE *f = fopen(filename, "r+");
+	while (ch != '\n') ch = fgetc(f);
+	ch = 'a';
+	while (ch != '\n' && ch != '\t') ch = fgetc(f);
+	if (ch == '\n') eigenpres = false;
+
+	rewind(f);
+
+	fscanf(f, "%le\t", &c);
+	fscanf(f, "%le\n", &energy);
+
+	if (eigenpres) {
+		ch = 'a';
+		while (ch != '\n') ch = fgetc(f);
+	}
+	gsl_vector_fscanf(f, z);
+	fclose(f);
+
+	min_fails = domain_minimize(z, n, c, eps, g, N, 4, 2, 0.9);
+
+	if (min_fails) {
+		printf("BIFUR: Initial relaxation failed, exiting.\n");
+		return 1;
+	}
+
+	bifur_eigenvalues(eigenvalues, n, z, c);
+	bifur_eigensort(eigenorder, n, num, eigenvalues);
+	bifur_trueEigenvalues(trueEigenvalues, n, z, c);
+	bifur_trueEigensort(trueEigenorder, n, num, trueEigenvalues);
+
+	energy = domain_energy_energy(n, z, c);
+	unsigned ii;
+
+	FILE *newf = fopen(filename, "w");
+	fprintf(newf, "%.12le\t%.12le\n", c, energy);
+	for (unsigned i = 0; i < num; i++) {
+		ii = gsl_permutation_get(eigenorder, i);
+		fprintf(newf, "%.12le\t", gsl_vector_get(eigenvalues, ii));
+	}
+	fprintf(newf, "\n");
+	for (unsigned i = 0; i < num; i++) {
+		ii = gsl_permutation_get(trueEigenorder, i);
+		fprintf(newf, "%.12le\t", gsl_vector_get(trueEigenvalues, ii));
+	}
+	fprintf(newf, "\n");
+	for (unsigned i = 0; i < 3 * n + 3; i++) {
+		fprintf(newf, "%.12le\t", gsl_vector_get(z, i));
+	}
+	fclose(newf);
+}
diff --git a/src/domain_eigen.cpp b/src/domain_eigen.cpp
new file mode 100644
index 0000000..c02240e
--- /dev/null
+++ b/src/domain_eigen.cpp
@@ -0,0 +1,187 @@
+/* domain_eigen.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* A set of utilities for find the generalized eigenvalues and eigenvectors of
+ * modulated domains.
+ */
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+// Finds the generalized eigenvalues of the Hessian for the state vector z when
+// Lambda = c.
+void domain_eigen_values(gsl_vector *eigenvalues, unsigned size, unsigned params, gsl_matrix *hess) {
+
+	double eigenvalue;
+
+	gsl_vector *beta;
+	gsl_vector_complex *alpha;
+	gsl_matrix *modI;
+	gsl_eigen_gen_workspace *w;
+
+	alpha = gsl_vector_complex_alloc(size);
+	beta = gsl_vector_alloc(size);
+	modI = gsl_matrix_alloc(size, size);
+	w = gsl_eigen_gen_alloc(size);
+
+	gsl_matrix_set_zero(modI);
+	for (unsigned i = 0; i < params; i++) gsl_matrix_set(modI, i, i, 1);
+
+	gsl_eigen_gen(hess, modI, alpha, beta, w);
+
+	for (unsigned i = 0; i < size; i++) {
+		eigenvalue = gsl_vector_complex_get(alpha, i).dat[0] / gsl_vector_get(beta, i);
+		gsl_vector_set(eigenvalues, i, eigenvalue);
+	}
+
+	gsl_vector_free(beta);
+	gsl_vector_complex_free(alpha);
+	gsl_matrix_free(modI);
+	gsl_eigen_gen_free(w);
+}
+
+
+void domain_eigen_sort(gsl_permutation *eigenorder, unsigned size, unsigned eigen_num,
+		const gsl_vector *eigenvalues) {
+
+	unsigned ii;
+
+	gsl_vector *abs_eigenvalues;
+
+	abs_eigenvalues = gsl_vector_alloc(size);
+	
+	for (unsigned i = 0; i < size; i++) {
+		gsl_vector_set(abs_eigenvalues, i, fabs(gsl_vector_get(eigenvalues, i)));
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_memcpy(abs_eigenvalues, eigenvalues);
+
+	for (unsigned i = eigen_num; i < size; i++) {
+		ii = gsl_permutation_get(eigenorder, i);
+
+		gsl_vector_set(abs_eigenvalues, ii, INFINITY);
+	}
+
+	gsl_sort_vector_index(eigenorder, abs_eigenvalues);
+
+	gsl_vector_free(abs_eigenvalues);
+}
+
+
+void domain_eigen_state(gsl_vector *eigenstate, const gsl_vector *eigenvalues,
+		unsigned n, double thres) {
+
+	double eigenvalue;
+
+	gsl_vector_set_zero(eigenstate);
+
+	for (unsigned i = 0; i < 3 * n + 3; i++) {
+		eigenvalue = gsl_vector_get(eigenvalues, i);
+		if (eigenvalue > fabs(thres)) gsl_vector_set(eigenstate, i, 1);
+		if (eigenvalue < -fabs(thres)) gsl_vector_set(eigenstate, i, -1);
+	}
+}
+
+
+void domain_eigen_vector(gsl_vector *eigenvector, unsigned size, unsigned params, unsigned k, gsl_matrix *hess) {
+
+	gsl_vector *beta;
+	gsl_vector_complex *alpha;
+	gsl_matrix *modI;
+	gsl_matrix_complex *evec;
+	gsl_eigen_genv_workspace *w;
+
+	alpha = gsl_vector_complex_alloc(size);
+	beta = gsl_vector_alloc(size);
+	modI = gsl_matrix_alloc(size, size);
+	evec = gsl_matrix_complex_alloc(size, size);
+	w = gsl_eigen_genv_alloc(size);
+
+	gsl_matrix_set_zero(modI);
+
+	for (unsigned i = 0; i < params; i++) gsl_matrix_set(modI, i, i, 1);
+
+	gsl_eigen_genv(hess, modI, alpha, beta, evec, w);
+
+	for (unsigned i = 0; i < size; i++) {
+		gsl_vector_set(eigenvector, i,
+			gsl_matrix_complex_get(evec, i, k).dat[0]);
+	}
+
+	gsl_vector_free(beta);
+	gsl_vector_complex_free(alpha);
+	gsl_matrix_free(modI);
+	gsl_matrix_complex_free(evec);
+	gsl_eigen_genv_free(w);
+}
+
+
+void domain_eigen_ortho(gsl_vector *eigenvector, unsigned n, const gsl_vector *z) {
+	gsl_vector *rotation, *translation_x, *translation_y;
+	double x, y, prod;
+
+	rotation = gsl_vector_alloc(3 * n + 3);
+	translation_x = gsl_vector_alloc(3 * n + 3);
+	translation_y = gsl_vector_alloc(3 * n + 3);
+
+	for (unsigned i = 0; i < n; i++) {
+		x = gsl_vector_get(z, i);
+		y = 0;
+
+		gsl_vector_set(translation_x, i, 1.0 / n);
+
+		if (n != 0) {
+			y = gsl_vector_get(z, n + i - 1);
+			gsl_vector_set(translation_y, n + i - 1, 1.0 / n);
+		}
+
+		gsl_vector_set(rotation, i, - y / (gsl_pow_2(x) + gsl_pow_2(y)));
+		if (n != 0) gsl_vector_set(rotation, n + i - 1, x / (n * (gsl_pow_2(x) + gsl_pow_2(y))));
+	}
+
+	gsl_blas_ddot(rotation, eigenvector, &prod);
+	prod = prod / gsl_blas_dnrm2(rotation);
+	gsl_vector_memcpy(rotation, eigenvector);
+	gsl_blas_daxpy(-prod, rotation, eigenvector);
+
+	gsl_blas_ddot(translation_x, eigenvector, &prod);
+	prod = prod / gsl_blas_dnrm2(translation_x);
+	gsl_vector_memcpy(translation_x, eigenvector);
+	gsl_blas_daxpy(-prod, translation_x, eigenvector);
+
+	gsl_blas_ddot(translation_y, eigenvector, &prod);
+	prod = prod / gsl_blas_dnrm2(translation_y);
+	gsl_vector_memcpy(translation_y, eigenvector);
+	gsl_blas_daxpy(-prod, translation_y, eigenvector);
+
+	return;
+}
+
+
diff --git a/src/domain_eigen.h b/src/domain_eigen.h
new file mode 100644
index 0000000..38366b0
--- /dev/null
+++ b/src/domain_eigen.h
@@ -0,0 +1,28 @@
+#ifndef DOMAIN_EIGEN_H
+#define DOMAIN_EIGEN_H
+
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+void domain_eigen_values(gsl_vector *eigenvalues, unsigned size, unsigned params,
+		gsl_matrix *hess);
+
+void domain_eigen_sort(gsl_permutation *eigenorder, unsigned size, unsigned eigen_num,
+		const gsl_vector *eigenvalues);
+
+void domain_eigen_state(gsl_vector *eigenstate, const gsl_vector *eigenvalues,
+		unsigned n, double thres);
+
+void domain_eigen_vector(gsl_vector *eigenvector, unsigned size, unsigned params, unsigned k, gsl_matrix *hess);
+
+void domain_eigen_ortho(gsl_vector *eigenvector, unsigned n, const gsl_vector *z);
+
+#endif
diff --git a/src/domain_energy.cpp b/src/domain_energy.cpp
new file mode 100644
index 0000000..11d6bfe
--- /dev/null
+++ b/src/domain_energy.cpp
@@ -0,0 +1,2103 @@
+/* domain_energy.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* A lightweight and efficient model of two dimensional phase-modulated domains
+ * using Newton's method to minimize a Lagrangian.
+ *
+ * Best viewed in an 80 character terminal with two character hard tabs.
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_blas.h>
+
+
+void domain_energy_x(gsl_vector *x, unsigned n, const gsl_vector *z) {
+// Gets the full set of x coordinates from the state array.
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+
+}
+
+
+void domain_energy_y(gsl_vector *y, unsigned n, const gsl_vector *z) {
+// Gets the full set of y coordinates from the state array.
+
+	gsl_vector_set(y, 0, 0);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < n; i++) {
+		gsl_vector_set(y, i, gsl_vector_get(z, n + i - 1));
+	}
+
+}
+
+
+void gsl_permutation_over(unsigned n, gsl_permutation *perm, bool right) {
+// Shifts a GSL permutation object circularly.  If right is true, then the
+// permutation is shifted to the right; if false, then it is shifted to the
+// left.
+
+	gsl_permutation_swap(perm, 0, n - 1);
+
+	if (right) {
+		for (unsigned i = 0; i < n - 2; i++) {
+			gsl_permutation_swap(perm, n - 1 - i, n - 2 - i);
+		}
+	}
+
+	else {
+		for (unsigned i = 0; i < n - 2; i++) {
+			gsl_permutation_swap(perm, i, i + 1);
+		}
+	}
+}
+
+
+double domain_energy_area(unsigned n, const gsl_vector *x, const gsl_vector *y) {
+// Computes the area of a domain.
+
+	double area, x_i, y_i, x_ii, y_ii;
+	unsigned ii;
+
+	gsl_permutation *indices_left;
+
+	indices_left = gsl_permutation_alloc(n);
+	gsl_permutation_init(indices_left);
+	gsl_permutation_over(n, indices_left, false);
+
+	area = 0;
+
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_left, i);
+
+		x_i = gsl_vector_get(x, i);
+		y_i = gsl_vector_get(y, i);
+		x_ii = gsl_vector_get(x, ii);
+		y_ii = gsl_vector_get(y, ii);
+
+		area += (x_i * y_ii - x_ii * y_i) / 2;
+	}
+
+	gsl_permutation_free(indices_left);
+
+	return area;
+}
+
+
+void domain_energy_rt(gsl_vector *rt, unsigned n, const gsl_vector *x,
+		double p) {
+// Converts x and y coordinates to r_x, r_y, t_x, or t_y, depending on input.
+
+	double x_i, x_ii;
+	unsigned ii;
+
+	gsl_permutation *indices_left;
+
+	indices_left = gsl_permutation_alloc(n);
+	gsl_permutation_init(indices_left);
+	gsl_permutation_over(n, indices_left, false);
+
+	#pragma omp parallel for private(ii, x_i, x_ii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_left, i);
+
+		x_i = gsl_vector_get(x, i);
+		x_ii = gsl_vector_get(x, ii);
+
+		gsl_vector_set(rt, i, x_ii + p * x_i);
+	}
+
+	gsl_permutation_free(indices_left);
+}
+
+
+void domain_energy_tdots(gsl_matrix *tdots, unsigned n, const gsl_vector *tx,
+		const gsl_vector *ty) {
+// Creates a matrix of dotted tangent vectors.
+
+	gsl_matrix_set_zero(tdots);
+
+	gsl_blas_dger(1, tx, tx, tdots);
+	gsl_blas_dger(1, ty, ty, tdots);
+}
+
+
+void domain_energy_dists(gsl_matrix *dists, unsigned n, const gsl_vector *rx,
+		const gsl_vector *ry) {
+// Creates a matrix of distances between points on the domain boundary..
+
+	double rx_i, rx_j, ry_i, ry_j;
+
+	#pragma omp parallel for private(rx_i, rx_j, ry_i, ry_j)
+	for (unsigned i = 0; i < n; i++) {
+		for (unsigned j = 0; j < n; j++) {
+			rx_i = gsl_vector_get(rx, i);
+			rx_j = gsl_vector_get(rx, j);
+
+			ry_i = gsl_vector_get(ry, i);
+			ry_j = gsl_vector_get(ry, j);
+
+			gsl_matrix_set(dists, i, j,
+					2 / sqrt(gsl_pow_2(rx_i - rx_j) + gsl_pow_2(ry_i - ry_j)));
+		}
+	}
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		gsl_matrix_set(dists, i, i, 0);
+	}
+}
+
+
+double domain_energy_dconst(unsigned n, double L, double *ld,
+		const gsl_vector *tx, const gsl_vector *ty) {
+// Computes the value of the Lagrangian constraint on the distances.
+
+	double dconst, tx_i, ty_i;
+
+	dconst = 0;
+
+	for (unsigned i = 0; i < n; i++) {
+		tx_i = gsl_vector_get(tx, i);
+		ty_i = gsl_vector_get(ty, i);
+
+		dconst += ld[i] * (gsl_pow_2(L / n)
+			- (gsl_pow_2(tx_i) + gsl_pow_2(ty_i)));
+	}
+
+	return dconst;
+}
+
+
+double domain_energy_dval(unsigned n, double L, double *ld,
+		const gsl_vector *tx, const gsl_vector *ty) {
+// Computes the value of the Lagrangian constraint on the distances.
+
+	double dconst, tx_i, ty_i;
+
+	dconst = 0;
+
+	for (unsigned i = 0; i < n; i++) {
+		tx_i = gsl_vector_get(tx, i);
+		ty_i = gsl_vector_get(ty, i);
+
+		dconst += gsl_pow_2(gsl_pow_2(L / n)
+			- (gsl_pow_2(tx_i) + gsl_pow_2(ty_i)));
+	}
+
+	return dconst;
+}
+
+
+double domain_energy_init(unsigned n, const gsl_vector *z, gsl_vector *rx,
+		gsl_vector *ry, gsl_vector *tx, gsl_vector *ty, gsl_matrix *tdots,
+		gsl_matrix *dists) {
+// Get useful objects from the state vector.  Fills all GSL matrix and vector
+// objects, and returns the value of the area.
+
+	double area;
+
+	gsl_vector *x, *y;
+	
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	domain_energy_x(x, n, z);
+	domain_energy_y(y, n, z);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	area = domain_energy_area(n, x, y);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+
+	return area;
+}
+
+
+double domain_energy_energy(unsigned n,
+		double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, gsl_matrix *tdots, gsl_matrix *dists, double L) {
+// Computes the Lagrangian.
+
+	double energy_mu, energy, lagrangian, H;
+ 
+	gsl_vector *v_temp_a, *v_temp_b;
+
+	v_temp_a = gsl_vector_alloc(n);
+	v_temp_b = gsl_vector_alloc(n);
+
+	gsl_vector_set_all(v_temp_a, 1);
+	gsl_vector_set_all(v_temp_b, 1);
+
+	gsl_matrix_mul_elements(tdots, dists);
+	gsl_blas_dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, tdots, v_temp_a);
+	gsl_blas_ddot(v_temp_a, v_temp_b, &energy_mu);
+
+	H = log(0.5 * (n - 1)) +  M_EULER + 1.0 / (n - 1)
+		- 1.0 / (12 * gsl_pow_2(0.5 * (n - 1)));
+
+	L = fabs(L);
+
+	energy = c * L - L * log(L) - energy_mu + L * H;
+
+	gsl_vector_free(v_temp_a);
+	gsl_vector_free(v_temp_b);
+
+	return energy;
+}
+
+double domain_energy_lagrangian(unsigned n,
+		double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, gsl_matrix *tdots, gsl_matrix *dists, double L, double area, double la, double *ld) {
+// Computes the Lagrangian.
+
+	double energy, lagrangian;
+
+	energy = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L);
+
+	L = fabs(L);
+
+	lagrangian = energy - la * (area - M_PI) - domain_energy_dconst(n, L, ld, tx, ty);
+
+	return lagrangian;
+}
+
+
+void domain_energy_gradient(gsl_vector *grad, unsigned n,
+		double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, const gsl_matrix *tdots, const gsl_matrix *dists, double L, double area, double la, double *ld) {
+
+// Computes the gradient of the Lagrangian.
+	double rx_i, rx_ii, rx_j, tx_i, tx_ii, tx_j, ry_i, ry_ii,
+	ry_j, ty_i, ty_ii, ty_j, d_ij, d_iij, tdt_ij, tdt_iij, d_ij3, d_iij3;
+	unsigned ii, jj;
+
+	gsl_vector *v_ones, *v_storage;
+	gsl_matrix *m_dx, *m_dy;
+	gsl_permutation *indices_right;
+	gsl_permutation *indices_left;
+
+	v_ones = gsl_vector_alloc(n);
+	v_storage = gsl_vector_alloc(n);
+
+	m_dx = gsl_matrix_alloc(n, n);
+	m_dy = gsl_matrix_alloc(n, n);
+
+	indices_right = gsl_permutation_alloc(n);
+	indices_left = gsl_permutation_alloc(n);
+
+	gsl_vector_set_zero(grad);
+	gsl_vector_set_all(v_ones, 1);
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+	gsl_permutation_init(indices_left);
+	gsl_permutation_over(n, indices_left, false);
+
+	#pragma omp parallel for private(rx_i, rx_j, tx_j, tdt_ij, d_ij, rx_ii, ry_i, ry_j, ty_j, ry_ii, tdt_iij, d_iij, d_ij3, d_iij3, ii)
+	for (unsigned i = 0; i < n; i++) {
+		for (unsigned j = 0; j < n; j++) {
+			ii = gsl_permutation_get(indices_right, i);
+
+			rx_i = gsl_vector_get(rx, i);
+			rx_ii = gsl_vector_get(rx, ii);
+			rx_j = gsl_vector_get(rx, j);
+			tx_j = gsl_vector_get(tx, j);
+
+			ry_i = gsl_vector_get(ry, i);
+			ry_ii = gsl_vector_get(ry, ii);
+			ry_j = gsl_vector_get(ry, j);
+			ty_j = gsl_vector_get(ty, j);
+
+			d_ij = gsl_matrix_get(dists, i, j);
+			d_iij = gsl_matrix_get(dists, ii, j);
+			tdt_ij = gsl_matrix_get(tdots, i, j);
+			tdt_iij = gsl_matrix_get(tdots, ii, j);
+
+			d_ij3 = gsl_pow_3(d_ij);
+			d_iij3 = gsl_pow_3(d_iij);
+
+			gsl_matrix_set(m_dx, i, j,
+					- tx_j * d_ij - (rx_i - rx_j) * tdt_ij * d_ij3 / 4
+					+ tx_j * d_iij - (rx_ii - rx_j) * tdt_iij * d_iij3 / 4);
+
+			gsl_matrix_set(m_dy, i, j,
+					- ty_j * d_ij - (ry_i - ry_j) * tdt_ij * d_ij3 / 4
+					+ ty_j * d_iij - (ry_ii - ry_j) * tdt_iij * d_iij3 / 4);
+		}
+	}
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+		gsl_matrix_set(m_dx, i, ii, 0);
+		gsl_matrix_set(m_dx, i, i, 0);
+		gsl_matrix_set(m_dy, i, ii, 0);
+		gsl_matrix_set(m_dy, i, i, 0);
+	}
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dx, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, tx_i, tx_ii, d_ij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+
+		gsl_vector_set(grad, i,
+				- (gsl_vector_get(v_storage, i) + (tx_i - tx_ii) * d_ij));
+	}
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dy, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, ty_i, ty_ii, d_ij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+
+		gsl_vector_set(grad, i + n,
+				- (gsl_vector_get(v_storage, i) + (ty_i - ty_ii) * d_ij));
+	}
+
+	// darea/dx_i or y_i
+
+	#pragma omp parallel for private(ii, jj)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+		jj = gsl_permutation_get(indices_left, i);
+
+		gsl_vector_set(grad, i, gsl_vector_get(grad, i) -
+				0.5 * (gsl_vector_get(ty, i) + gsl_vector_get(ty, ii)) * la);
+
+		gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) +
+				0.5 * (gsl_vector_get(tx, i) + gsl_vector_get(tx, ii)) * la);
+	}
+
+	#pragma omp parallel for private(ii, jj, tx_i, tx_ii, ty_i, ty_ii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+		jj = gsl_permutation_get(indices_left, i);
+
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+
+		gsl_vector_set(grad, i, gsl_vector_get(grad, i) -
+				2 * (ld[i] * tx_i - ld[ii] * tx_ii));
+
+		gsl_vector_set(grad, i + n, gsl_vector_get(grad, i + n) -
+				2 * (ld[i] * ty_i - ld[ii] * ty_ii));
+	}
+
+	// The gradient with respect to L.
+
+	L = fabs(L);
+	double gradLDist = 0;
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+		jj = gsl_permutation_get(indices_left, i);
+
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+
+		gradLDist += 2 * L / gsl_pow_2(n) * ld[i];
+	}
+
+	double H = log(0.5 * (n - 1)) +  M_EULER + 1.0 / (n - 1)
+		- 1.0 / (12 * gsl_pow_2(0.5 * (n - 1)));
+
+	double gradL = GSL_SIGN(L) * (c - (1 + log(L) - H) - gradLDist);
+
+	gsl_vector_set(grad, 2 * n, gradL);
+
+	// The gradients with respect to the undetermined coefficients are simply
+	// the constraints.
+
+	double gradla = M_PI - area;
+
+	gsl_vector_set(grad, 2 * n + 1, gradla);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		gsl_vector_set(grad, 2 * n + 2 + i, gsl_pow_2(gsl_vector_get(tx, i)) +
+				gsl_pow_2(gsl_vector_get(ty, i)) - gsl_pow_2(L / n));
+	}
+
+	gsl_vector_free(v_ones);
+	gsl_vector_free(v_storage);
+
+	gsl_matrix_free(m_dx);
+	gsl_matrix_free(m_dy);
+
+	gsl_permutation_free(indices_right);
+	gsl_permutation_free(indices_left);
+}
+
+
+void domain_energy_halfHessian(gsl_matrix *hess, unsigned n,
+		double c, const gsl_vector *rx, const gsl_vector *ry, const gsl_vector *tx, const gsl_vector *ty, const gsl_matrix *tdots, const gsl_matrix *dists, double L, double la, double *ld) {
+/* Computes the Hessian of the Lagrangian without the center of mass
+ * constraints and fixed point.
+ */
+
+	gsl_matrix *m_dxidxj, *m_dyidyj, *m_dxidxi, *m_dyidyi, *m_dxidxii,
+						 *m_dyidyii, *m_dxidyj, *m_dxidyi, *m_dxiidyi, *m_dxidyii;
+	gsl_vector *v_ones, *v_storage;
+	gsl_permutation *indicesRight, *indicesLeft, *indices2Right;
+
+	unsigned ii, jj, iii;
+	double rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj,
+				 ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij,
+				 tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii,
+				 d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5;
+
+	m_dxidxj = gsl_matrix_alloc(n, n);
+	m_dyidyj = gsl_matrix_alloc(n, n);
+	m_dxidxi = gsl_matrix_alloc(n, n);
+	m_dyidyi = gsl_matrix_alloc(n, n);
+	m_dxidxii = gsl_matrix_alloc(n, n);
+	m_dyidyii = gsl_matrix_alloc(n, n);
+	m_dxidyj = gsl_matrix_alloc(n, n);
+	m_dxidyi = gsl_matrix_alloc(n, n);
+	m_dxiidyi = gsl_matrix_alloc(n, n);
+	m_dxidyii = gsl_matrix_alloc(n, n);
+
+	v_ones = gsl_vector_alloc(n);
+	v_storage = gsl_vector_alloc(n);
+
+	indicesRight = gsl_permutation_alloc(n);
+	indicesLeft = gsl_permutation_alloc(n);
+	indices2Right = gsl_permutation_alloc(n);
+
+	gsl_matrix_set_zero(hess);
+	gsl_vector_set_all(v_ones, 1);
+
+	gsl_permutation_init(indicesRight);
+	gsl_permutation_init(indicesLeft);
+	gsl_permutation_over(n, indicesRight, true);
+	gsl_permutation_over(n, indicesLeft, false);
+	gsl_permutation_memcpy(indices2Right, indicesRight);
+	gsl_permutation_over(n, indices2Right, true);
+
+	#pragma omp parallel for private(rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj, ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij, tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii, d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5, ii, jj)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+		for (unsigned j = 0; j < n; j++) {
+			jj = gsl_permutation_get(indicesRight, j);
+
+			rx_i = gsl_vector_get(rx, i);
+			rx_j = gsl_vector_get(rx, j);
+			tx_i = gsl_vector_get(tx, i);
+			tx_j = gsl_vector_get(tx, j);
+
+			ry_i = gsl_vector_get(ry, i);
+			ry_j = gsl_vector_get(ry, j);
+			ty_i = gsl_vector_get(ty, i);
+			ty_j = gsl_vector_get(ty, j);
+
+			rx_ii = gsl_vector_get(rx, ii);
+			rx_jj = gsl_vector_get(rx, jj);
+			tx_ii = gsl_vector_get(tx, ii);
+			tx_jj = gsl_vector_get(tx, jj);
+
+			ry_ii = gsl_vector_get(ry, ii);
+			ry_jj = gsl_vector_get(ry, jj);
+			ty_ii = gsl_vector_get(ty, ii);
+			ty_jj = gsl_vector_get(ty, jj);
+
+			d_ij = gsl_matrix_get(dists, i, j);
+			tdt_ij = gsl_matrix_get(tdots, i, j);
+
+			d_iij = gsl_matrix_get(dists, ii, j);
+			tdt_iij = gsl_matrix_get(tdots, ii, j);
+
+			d_ijj = gsl_matrix_get(dists, i, jj);
+			tdt_ijj = gsl_matrix_get(tdots, i, jj);
+
+			d_iijj = gsl_matrix_get(dists, ii, jj);
+			tdt_iijj = gsl_matrix_get(tdots, ii, jj);
+
+			d_ij3 = gsl_pow_3(d_ij);
+			d_iij3 = gsl_pow_3(d_iij);
+			d_ijj3 = gsl_pow_3(d_ijj);
+			d_iijj3 = gsl_pow_3(d_iijj);
+			d_ij5 = gsl_pow_5(d_ij);
+			d_iij5 = gsl_pow_5(d_iij);
+			d_ijj5 = gsl_pow_5(d_ijj);
+			d_iijj5 = gsl_pow_5(d_iijj);
+
+			// d^2E / dx_i dx_j for i != j, j - 1.
+			gsl_matrix_set(m_dxidxj, i, j,
+					(rx_i - rx_j) * (tx_i - tx_j) * d_ij3 / 4 + d_ij
+					- 3 * gsl_pow_2(rx_i - rx_j) * tdt_ij * d_ij5 / 16
+					+ tdt_ij * d_ij3 / 4
+
+					+ (rx_ii - rx_j) * (tx_ii + tx_j) * d_iij3 / 4 - d_iij
+					- 3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16
+					+ tdt_iij * d_iij3 / 4
+
+					- (rx_i - rx_jj) * (tx_i + tx_jj) * d_ijj3 / 4 - d_ijj
+					- 3 * gsl_pow_2(rx_i - rx_jj) * tdt_ijj * d_ijj5 / 16
+					+ tdt_ijj * d_ijj3 / 4
+
+					- (rx_ii - rx_jj) * (tx_ii - tx_jj) * d_iijj3 / 4 + d_iijj
+					- 3 * gsl_pow_2(rx_ii - rx_jj) * tdt_iijj * d_iijj5 / 16
+					+ tdt_iijj * d_iijj3 / 4
+					);
+
+			// d^2E / dy_i dy_j for i != j, j - 1.
+			gsl_matrix_set(m_dyidyj, i, j,
+					(ry_i - ry_j) * (ty_i - ty_j) * d_ij3 / 4 + d_ij
+					- 3 * gsl_pow_2(ry_i - ry_j) * tdt_ij * d_ij5 / 16
+					+ tdt_ij * d_ij3 / 4
+
+					+ (ry_ii - ry_j) * (ty_ii + ty_j) * d_iij3 / 4 - d_iij
+					- 3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+					+ tdt_iij * d_iij3 / 4
+
+					- (ry_i - ry_jj) * (ty_i + ty_jj) * d_ijj3 / 4 - d_ijj
+					- 3 * gsl_pow_2(ry_i - ry_jj) * tdt_ijj * d_ijj5 / 16
+					+ tdt_ijj * d_ijj3 / 4
+
+					- (ry_ii - ry_jj) * (ty_ii - ty_jj) * d_iijj3 / 4 + d_iijj
+					- 3 * gsl_pow_2(ry_ii - ry_jj) * tdt_iijj * d_iijj5 / 16
+					+ tdt_iijj * d_iijj3 / 4
+					);
+
+			// d^2E / dx_i^2
+			gsl_matrix_set(m_dxidxi, i, j,
+					(rx_i - rx_j) * tx_j * d_ij3 / 2
+					+ 3 * gsl_pow_2(rx_i - rx_j) * tdt_ij * d_ij5 / 16
+					- tdt_ij * d_ij3 / 4
+
+					- (rx_ii - rx_j) * tx_j * d_iij3 / 2
+					+ 3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16
+					- tdt_iij * d_iij3 / 4
+					);
+
+			// d^2E / dy_i^2
+			gsl_matrix_set(m_dyidyi, i, j,
+					(ry_i - ry_j) * ty_j * d_ij3 / 2
+					+ 3 * gsl_pow_2(ry_i - ry_j) * tdt_ij * d_ij5 / 16
+					- tdt_ij * d_ij3 / 4
+
+					- (ry_ii - ry_j) * ty_j * d_iij3 / 2
+					+ 3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+					- tdt_iij * d_iij3 / 4
+					);
+
+			// d^2E / dx_i dx_(i-1)
+			gsl_matrix_set(m_dxidxii, i, j,
+					3 * gsl_pow_2(rx_ii - rx_j) * tdt_iij * d_iij5 / 16
+					- tdt_iij * d_iij3 / 4
+					);
+
+			// d^2E / dy_i dy_(i-1)
+			gsl_matrix_set(m_dyidyii, i, j,
+					3 * gsl_pow_2(ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+					- tdt_iij * d_iij3 / 4
+					);
+
+			gsl_matrix_set(m_dxidyj, i, j,
+					(rx_i - rx_j) * ty_i * d_ij3 / 4
+					- (ry_i - ry_j) * tx_j * d_ij3 / 4
+					- 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16
+
+					+ (rx_ii - rx_j) * ty_ii * d_iij3 / 4
+					+ (ry_ii - ry_j) * tx_j * d_iij3 / 4
+					- 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+
+					- (rx_i - rx_jj) * ty_i * d_ijj3 / 4
+					- (ry_i - ry_jj) * tx_jj * d_ijj3 / 4
+					- 3 * (rx_i - rx_jj) * (ry_i - ry_jj) * tdt_ijj * d_ijj5 / 16
+
+					- (rx_ii - rx_jj) * ty_ii * d_iijj3 / 4
+					+ (ry_ii - ry_jj) * tx_jj * d_iijj3 / 4
+					- 3 * (rx_ii - rx_jj) * (ry_ii - ry_jj) * tdt_iijj * d_iijj5 / 16
+					);
+
+			gsl_matrix_set(m_dxidyi, i, j,
+					(ry_i - ry_j) * tx_j * d_ij3 / 4
+					+ (rx_i - rx_j) * ty_j * d_ij3 / 4
+					+ 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16
+
+					- (ry_ii - ry_j) * tx_j * d_iij3 / 4
+					- (rx_ii - rx_j) * ty_j * d_iij3 / 4
+					+ 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+					);
+
+			gsl_matrix_set(m_dxiidyi, i, j,
+					(ry_ii - ry_j) * tx_j * d_iij3 / 4
+					- (rx_ii - rx_j) * ty_j * d_iij3 / 4
+					+ 3 * (rx_ii - rx_j) * (ry_ii - ry_j) * tdt_iij * d_iij5 / 16
+					);
+
+			gsl_matrix_set(m_dxidyii, i, j,
+					- (ry_i - ry_j) * tx_j * d_ij3 / 4
+					+ (rx_i - rx_j) * ty_j * d_ij3 / 4
+					+ 3 * (rx_i - rx_j) * (ry_i - ry_j) * tdt_ij * d_ij5 / 16
+					);
+
+		}
+	}
+
+		// Setting terms of d^2E / dy_i dy_j  and d^2E / dx_i dx_j in the Hessian.
+	#pragma omp parallel for
+	for (unsigned i = 2; i < n; i++) {
+		for (unsigned j = 0; j < i - 1; j++) {
+
+			gsl_matrix_set(hess, i, j, - gsl_matrix_get(m_dxidxj, i, j));
+
+			gsl_matrix_set(hess, n + i, n + j, - gsl_matrix_get(m_dyidyj, i, j));
+		}
+	}
+
+	// Zeroing out terms which aren't supposed to appear in the sums.
+	#pragma omp parallel for private(ii, iii, jj)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+		iii = gsl_permutation_get(indices2Right, i);
+		jj = gsl_permutation_get(indicesLeft, i);
+		gsl_matrix_set(m_dxidxi, i, i, 0);
+		gsl_matrix_set(m_dxidxi, i, ii, 0);
+		gsl_matrix_set(m_dyidyi, i, i, 0);
+		gsl_matrix_set(m_dyidyi, i, ii, 0);
+		gsl_matrix_set(m_dxidxii, i, i, 0);
+		gsl_matrix_set(m_dxidxii, i, ii, 0);
+		gsl_matrix_set(m_dxidxii, i, iii, 0);
+		gsl_matrix_set(m_dyidyii, i, i, 0);
+		gsl_matrix_set(m_dyidyii, i, ii, 0);
+		gsl_matrix_set(m_dyidyii, i, iii, 0);
+		gsl_matrix_set(m_dxidyi, i, i, 0);
+		gsl_matrix_set(m_dxidyi, i, ii, 0);
+		gsl_matrix_set(m_dxiidyi, i, i, 0);
+		gsl_matrix_set(m_dxiidyi, i, ii, 0);
+		gsl_matrix_set(m_dxiidyi, i, iii, 0);
+		gsl_matrix_set(m_dxidyii, i, ii, 0);
+		gsl_matrix_set(m_dxidyii, i, jj, 0);
+		gsl_matrix_set(m_dxidyii, i, i, 0);
+	}
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dxidxi, v_ones, 0, v_storage);
+
+	// Setting terms of d^2E / dx_i^2 in the Hessian.
+	#pragma omp parallel for private(ii, d_iij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		d_iij = gsl_matrix_get(dists, i, ii);
+
+		gsl_matrix_set(hess, i, i,
+				- (gsl_vector_get(v_storage, i) - 2 * d_iij)
+				+ 2 * (ld[i] + ld[ii])
+				);
+	}
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dyidyi, v_ones, 0, v_storage);
+
+	// Setting terms of d^2E / dy_i^2 in the Hessian.
+	#pragma omp parallel for private(ii, d_iij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		d_iij = gsl_matrix_get(dists, i, ii);
+
+		gsl_matrix_set(hess, n + i, n + i,
+				- (gsl_vector_get(v_storage, i) - 2 * d_iij)
+				+ 2 * (ld[i] + ld[ii])
+				);
+	}
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dxidxii, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, rx_i, rx_ii, tx_i, tx_ii, d_ij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		rx_i = gsl_vector_get(rx, i);
+		rx_ii = gsl_vector_get(rx, ii);
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- (rx_ii - rx_i) * (tx_i - tx_ii) * gsl_pow_3(d_ij) / 4
+				+ d_ij);
+	}
+
+	#pragma omp parallel for private(ii, iii, rx_i, rx_ii, rx_iii, tx_i, tx_ii, tx_iii, d_ij, d_iij, d_iijj, tdt_ij)
+	for (unsigned i = 0; i < n; i++) {
+		iii = gsl_permutation_get(indices2Right, i);
+		ii = gsl_permutation_get(indicesRight, i);
+
+		rx_i = gsl_vector_get(rx, i);
+		rx_ii = gsl_vector_get(rx, ii);
+		rx_iii = gsl_vector_get(rx, iii);
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+		tx_iii = gsl_vector_get(tx, iii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+		d_iij = gsl_matrix_get(dists, i, iii);
+		d_iijj = gsl_matrix_get(dists, ii, iii);
+		tdt_ij = gsl_matrix_get(tdots, i, iii);
+
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- 3 * gsl_pow_2(rx_iii - rx_i) * tdt_ij * gsl_pow_5(d_iij) / 16
+				+ (tx_iii + tx_i) * (rx_iii - rx_i) * gsl_pow_3(d_iij) / 4
+				+ tdt_ij * gsl_pow_3(d_iij) / 4 - d_iij
+				- (tx_ii - tx_iii) * (rx_ii - rx_iii) * gsl_pow_3(d_iijj) / 4 + d_iijj);
+	}
+
+	gsl_matrix_set(hess, n - 1, 0,
+			- gsl_vector_get(v_storage, 0) - 2 * ld[n - 1]);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < n; i++) {
+		gsl_matrix_set(hess, i,  i - 1,
+				- gsl_vector_get(v_storage, i) - 2 * ld[i - 1]
+				);
+	}
+
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dyidyii, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, ry_i, ry_ii, ty_i, ty_ii, d_ij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		ry_i = gsl_vector_get(ry, i);
+		ry_ii = gsl_vector_get(ry, ii);
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- (ry_ii - ry_i) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4
+				+ d_ij);
+	}
+
+	#pragma omp parallel for private(ii, iii, ry_i, ry_ii, ry_iii, ty_i, ty_ii, ty_iii, d_ij, d_iij, d_iijj, tdt_ij)
+	for (unsigned i = 0; i < n; i++) {
+		iii = gsl_permutation_get(indices2Right, i);
+		ii = gsl_permutation_get(indicesRight, i);
+
+		ry_i = gsl_vector_get(ry, i);
+		ry_ii = gsl_vector_get(ry, ii);
+		ry_iii = gsl_vector_get(ry, iii);
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+		ty_iii = gsl_vector_get(ty, iii);
+		d_ij = gsl_matrix_get(dists, i, ii);
+		d_iij = gsl_matrix_get(dists, i, iii);
+		d_iijj = gsl_matrix_get(dists, ii, iii);
+		tdt_ij = gsl_matrix_get(tdots, i, iii);
+
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- 3 * gsl_pow_2(ry_iii - ry_i) * tdt_ij * gsl_pow_5(d_iij) / 16
+				+ (ty_iii + ty_i) * (ry_iii - ry_i) * gsl_pow_3(d_iij) / 4
+				+ tdt_ij * gsl_pow_3(d_iij) / 4 - d_iij
+				- (ty_ii - ty_iii) * (ry_ii - ry_iii) * gsl_pow_3(d_iijj) / 4 + d_iijj);
+	}
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < n; i++) {
+		gsl_matrix_set(hess, n + i,  n + i - 1,
+				- gsl_vector_get(v_storage, i) - 2 * ld[i - 1]
+				);
+	}
+
+	gsl_matrix_set(hess, 2 * n -1, n,
+			- gsl_vector_get(v_storage, 0) - 2 * ld[n-1]);
+
+
+	// dxdy boring style
+	#pragma omp parallel for
+	for (unsigned j = 2; j < n; j++) {
+		for (unsigned i = 0; i < j - 1; i++) {
+			gsl_matrix_set(hess, n + j, i,
+					- gsl_matrix_get(m_dxidyj, i, j));
+		}
+	}
+
+	#pragma omp parallel for
+	for (unsigned j = 0; j < n - 2; j++) {
+		for (unsigned i = j + 2; i < n; i++) {
+			gsl_matrix_set(hess, n + j, i,
+					- gsl_matrix_get(m_dxidyj, i, j));
+		}
+	}
+
+	// d^2E / dx_i dy_i
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dxidyi, v_ones, 0, v_storage);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		gsl_matrix_set(hess, n + i, i,
+				- gsl_vector_get(v_storage, i)
+				);
+	}
+
+	// d^2 E / dx_ii dy_i
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dxiidyi, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, iii, rx_i, rx_ii, rx_iii, tx_i, tx_ii, tx_iii, ry_i, ry_ii, ry_iii, ty_i, ty_ii, ty_iii, d_ij, d_iij, d_iijj, tdt_iij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+		iii = gsl_permutation_get(indices2Right, i);
+
+		rx_i = gsl_vector_get(rx, i);
+		rx_ii = gsl_vector_get(rx, ii);
+		rx_iii = gsl_vector_get(rx, iii);
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+		tx_iii = gsl_vector_get(tx, iii);
+		ry_i = gsl_vector_get(ry, i);
+		ry_ii = gsl_vector_get(ry, ii);
+		ry_iii = gsl_vector_get(ry, iii);
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+		ty_iii = gsl_vector_get(ty, iii);
+
+		d_ij = gsl_matrix_get(dists, i, ii);
+		d_iij = gsl_matrix_get(dists, i, iii);
+		d_iijj = gsl_matrix_get(dists, ii, iii);
+		tdt_iij = gsl_matrix_get(tdots, i, iii);
+
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- (rx_ii - rx_i) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4
+				- 3 * (rx_iii - rx_i) * (ry_iii - ry_i) * tdt_iij * gsl_pow_5(d_iij) / 16
+				+ (rx_iii - rx_i) * ty_iii * gsl_pow_3(d_iij) / 4
+				+ (ry_iii - ry_i) * tx_i * gsl_pow_3(d_iij) / 4
+				- (tx_ii - tx_iii) * (ry_ii - ry_iii) * gsl_pow_3(d_iijj) / 4
+				);
+	}
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n - 1; i++) {
+		gsl_matrix_set(hess, n + i + 1, i,
+				- gsl_vector_get(v_storage, i + 1) - la / 2
+				);
+	}
+
+	gsl_matrix_set(hess, n, n - 1,
+			-gsl_vector_get(v_storage, 0) - la / 2);
+
+	// Upper off-diagonal of dxdy submatrix.
+
+	gsl_blas_dgemv(CblasNoTrans, 1, m_dxidyii, v_ones, 0, v_storage);
+
+	#pragma omp parallel for private(ii, jj, rx_i, rx_ii, rx_jj, tx_i, tx_ii, tx_jj, ry_i, ry_ii, ry_jj, ty_i, ty_ii, ty_jj, d_ij, d_ijj, d_iijj, tdt_iijj)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+		jj = gsl_permutation_get(indicesLeft, i);
+
+		rx_i = gsl_vector_get(rx, i);
+		rx_ii = gsl_vector_get(rx, ii);
+		rx_jj = gsl_vector_get(rx, jj);
+		tx_i = gsl_vector_get(tx, i);
+		tx_ii = gsl_vector_get(tx, ii);
+		tx_jj = gsl_vector_get(tx, jj);
+		ry_i = gsl_vector_get(ry, i);
+		ry_ii = gsl_vector_get(ry, ii);
+		ry_jj = gsl_vector_get(ry, jj);
+		ty_i = gsl_vector_get(ty, i);
+		ty_ii = gsl_vector_get(ty, ii);
+		ty_jj = gsl_vector_get(ty, jj);
+
+		d_ij = gsl_matrix_get(dists, i, ii);
+		d_ijj = gsl_matrix_get(dists, i, jj);
+		d_iijj = gsl_matrix_get(dists, ii, jj);
+		tdt_iijj = gsl_matrix_get(tdots, ii, jj);
+
+		gsl_vector_set(v_storage, i, gsl_vector_get(v_storage, i)
+				- (rx_i - rx_ii) * (ty_i - ty_ii) * gsl_pow_3(d_ij) / 4
+				- (tx_jj - tx_i) * (ry_i - ry_jj) * gsl_pow_3(d_ijj) / 4
+				+ tx_ii * (ry_ii - ry_jj) * gsl_pow_3(d_iijj) / 4
+				- (rx_jj - rx_ii) * ty_jj * gsl_pow_3(d_iijj) / 4
+				+ 3 * (rx_jj - rx_ii) * (ry_ii - ry_jj) * tdt_iijj * gsl_pow_5(d_iijj) / 16
+				);
+	}
+
+	gsl_matrix_set(hess, 2 * n - 1, 0,
+			- gsl_vector_get(v_storage, n - 1) + la / 2
+			);
+
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n-1; i++) {
+		gsl_matrix_set(hess, n + i, i + 1,
+				-  gsl_vector_get(v_storage, i) + la / 2
+				);
+	}
+
+
+	// dLdL
+
+	double gradLDist = 0;
+	for (unsigned i = 0; i < n; i++) {
+		gradLDist += 2 * ld[i] / gsl_pow_2(n);
+	}
+
+	L = fabs(L);
+
+	double gradL = - 1 / L - gradLDist;
+
+	gsl_matrix_set(hess, 2 * n, 2 * n, gradL);
+
+	// dLdlambdad
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		gsl_matrix_set(hess, i + 2 * n + 2, 2 * n,  - 2 * L / gsl_pow_2(n));
+	}
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		gsl_matrix_set(hess, 2 * n + 1, i, -
+				0.5 * (gsl_vector_get(ty, i) + gsl_vector_get(ty, ii)));
+
+		gsl_matrix_set(hess, 2 * n + 1,  n + i, 
+				0.5 * (gsl_vector_get(tx, i) + gsl_vector_get(tx, ii)));
+	}
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		gsl_matrix_set(hess, 2 * n + 2 + i, i, -
+				2 * (gsl_vector_get(tx, i)));
+
+		gsl_matrix_set(hess, 2 * n + 2 + i, n + i, -
+				2 * (gsl_vector_get(ty, i)));
+	}
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i = 1; i < n; i++) {
+		ii = gsl_permutation_get(indicesRight, i);
+
+		gsl_matrix_set(hess, 2 * n + i + 1, i,
+				2 * gsl_vector_get(tx, ii));
+
+		gsl_matrix_set(hess, 2 * n + i + 1, n + i,
+				2 * gsl_vector_get(ty, ii));
+	}
+
+	gsl_matrix_set(hess, 3 * n + 1, 0,
+			2 * gsl_vector_get(tx, n-1));
+	gsl_matrix_set(hess, 3 * n + 1, n,
+			2 * gsl_vector_get(ty, n-1));
+
+
+	gsl_vector_free(v_ones);
+	gsl_vector_free(v_storage);
+
+	gsl_matrix_free(m_dxidxj);
+	gsl_matrix_free(m_dyidyj);
+	gsl_matrix_free(m_dxidxi);
+	gsl_matrix_free(m_dyidyi);
+	gsl_matrix_free(m_dxidxii);
+	gsl_matrix_free(m_dyidyii);
+	gsl_matrix_free(m_dxidyj);
+	gsl_matrix_free(m_dxidyi);
+	gsl_matrix_free(m_dxiidyi);
+	gsl_matrix_free(m_dxidyii);
+
+	gsl_permutation_free(indicesLeft);
+	gsl_permutation_free(indicesRight);
+	gsl_permutation_free(indices2Right);
+}
+
+
+double domain_energy_nakedEnergy(unsigned n, const gsl_vector *z, double c) {
+	double lagrangian;
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n);
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	lagrangian = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+
+	return lagrangian;
+};
+
+
+double domain_energy_nakedLagrangian(unsigned n, const gsl_vector *z, double c) {
+	double lagrangian;
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n);
+	double la = gsl_vector_get(z, 2 * n + 1);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 2 + i);
+	}
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	double area = domain_energy_area(n, x, y);
+
+	lagrangian = domain_energy_lagrangian(n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+
+	return lagrangian;
+};
+
+
+void domain_energy_nakedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c) {
+
+	gsl_permutation *indices_right;
+	indices_right = gsl_permutation_alloc(n);
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	// Setting pointers to give the elements of z more convenient names.
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n);
+	double la = gsl_vector_get(z, 2 * n + 1);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 2 + i);
+	}
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	double area = domain_energy_area(n, x, y);
+
+	domain_energy_gradient(grad, n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+	gsl_permutation_free(indices_right);
+
+}
+
+
+void domain_energy_nakedHalfHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c) {
+
+	gsl_permutation *indices_right;
+	indices_right = gsl_permutation_alloc(n);
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	// Setting pointers to give the elements of z more convenient names.
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n);
+	double la = gsl_vector_get(z, 2 * n + 1);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 2 + i);
+	}
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	double area = domain_energy_area(n, x, y);
+
+	domain_energy_halfHessian(hess, n, c, rx, ry, tx, ty, tdots, dists, L, la, ld);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+	gsl_permutation_free(indices_right);
+}
+
+
+void domain_energy_nakedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c) {
+
+	domain_energy_nakedHalfHessian(hess, n, z, c);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < 3 * n + 2; i++) {
+		for (unsigned j = 0; j < i; j++) {
+			gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j));
+		}
+	}
+}
+
+
+double domain_energy_fixedEnergy(unsigned n, const gsl_vector *z, double c) {
+	double lagrangian;
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	gsl_vector_set(y, 0, 0);
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n - 1; i++) gsl_vector_set(y, i + 1, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n - 1);
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	lagrangian = domain_energy_energy(n, c, rx, ry, tx, ty, tdots, dists, L);
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+
+	return lagrangian;
+};
+
+
+// The fixed functions.
+
+double domain_energy_fixedLagrangian(unsigned n, const gsl_vector *z, double c) {
+	double lagrangian;
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	gsl_vector_set(y, 0, 0);
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n - 1; i++) gsl_vector_set(y, i + 1, gsl_vector_get(z, i + n));
+	double L = gsl_vector_get(z, 2 * n - 1);
+	double la = gsl_vector_get(z, 2 * n);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 1 + i);
+	}
+	double lx = gsl_vector_get(z, 3 * n + 1);
+	double ly = gsl_vector_get(z, 3 * n + 2);
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	domain_energy_rt(rx, n, x, 1);
+	domain_energy_rt(ry, n, y, 1);
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	domain_energy_tdots(tdots, n, tx, ty);
+	domain_energy_dists(dists, n, rx, ry);
+
+	double area = domain_energy_area(n, x, y);
+
+	lagrangian = domain_energy_lagrangian(n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld);
+
+	double xtot = 0;
+	for (unsigned i = 0; i < n; i++) xtot += gsl_vector_get(z, i);
+
+	double ytot = 0;
+	for (unsigned i = 1; i < n; i++) ytot += gsl_vector_get(z, i + n - 1);
+
+	lagrangian += - lx * xtot - ly * ytot;
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+
+	return lagrangian;
+};
+
+
+void domain_energy_fixedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z,
+		double c) {
+
+	// Setting pointers to give the elements of z more convenient names.
+	double L = gsl_vector_get(z, 2 * n - 1);
+	double la = gsl_vector_get(z, 2 * n);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 1 + i);
+	}
+	double lx = gsl_vector_get(z, 3 * n + 1);
+	double ly = gsl_vector_get(z, 3 * n + 2);
+
+	unsigned ii;
+
+	gsl_vector *rx, *ry, *tx, *ty, *freegrad;
+	gsl_matrix *tdots, *dists;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	freegrad = gsl_vector_alloc(3 * n + 2);
+
+	double area = domain_energy_init(n, z, rx, ry, tx, ty, tdots, dists);
+
+	domain_energy_gradient(freegrad, n, c, rx, ry, tx, ty, tdots, dists, L, area, la, ld);
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i =0; i < 3 * n + 2; i++) {
+		if (i != n) {
+			if (i < n) ii = i;
+			if (i > n) ii = i - 1;
+
+			gsl_vector_set(grad, ii, gsl_vector_get(freegrad, i));
+		}
+	}
+
+
+	#pragma omp parallel for private(ii)
+	for (unsigned i = 0; i < n; i++) {
+		gsl_vector_set(grad, i, gsl_vector_get(grad, i) - lx);
+		if (i != 0) gsl_vector_set(grad, n + i - 1, gsl_vector_get(grad, i + n - 1) -ly);
+	}
+
+
+	double xtot = 0;
+	for (unsigned i = 0; i < n; i++) xtot += gsl_vector_get(z, i);
+	double ytot = 0;
+	for (unsigned i = 1; i < n; i++) ytot += gsl_vector_get(z, i + n - 1);
+
+	gsl_vector_set(grad, 3 * n + 1, -xtot);
+
+	gsl_vector_set(grad, 3 * n + 2, -ytot);
+
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_vector_free(freegrad);
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+
+}
+
+
+void domain_energy_fixedHalfHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z,
+		double c) {
+
+	// Setting pointers to give the elements of z more convenient names.
+	double L = gsl_vector_get(z, 2 * n - 1);
+	double la = gsl_vector_get(z, 2 * n);
+	double ld[n];
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		ld[i] = gsl_vector_get(z, 2 * n + 1 + i);
+	}
+	double lx = gsl_vector_get(z, 3 * n + 1);
+	double ly = gsl_vector_get(z, 3 * n + 2);
+
+	gsl_vector *rx, *ry, *tx, *ty;
+	gsl_matrix *tdots, *dists, *freehess;
+
+	rx = gsl_vector_alloc(n);
+	ry = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	tdots = gsl_matrix_alloc(n, n);
+	dists = gsl_matrix_alloc(n, n);
+
+	freehess = gsl_matrix_alloc(3 * n + 2, 3 * n + 2);
+
+	double area = domain_energy_init(n, z, rx, ry, tx, ty, tdots, dists);
+
+	domain_energy_halfHessian(freehess, n, c, rx, ry, tx, ty, tdots, dists, L, la, ld);
+
+	gsl_matrix *m_dxidxj, *m_dyidyj, *m_dxidxi, *m_dyidyi, *m_dxidxii,
+						 *m_dyidyii, *m_dxidyj, *m_dxidyi, *m_dxiidyi, *m_dxidyii;
+	gsl_vector *v_ones, *v_storage;
+	gsl_permutation *indicesRight, *indicesLeft, *indices2Right;
+
+	unsigned ii, jj, iii;
+	double rx_i, rx_j, tx_i, tx_j, tdt_ij, d_ij, rx_ii, rx_jj, tx_ii, tx_jj,
+				 ry_i, ry_j, ty_i, ty_j, ry_ii, ry_jj, ty_ii, ty_jj, tdt_iij, d_iij,
+				 tdt_ijj, d_ijj, tdt_iijj, d_iijj, rx_iii, tx_iii, ry_iii, ty_iii,
+				 d_ij3, d_iij3, d_ijj3, d_iijj3, d_ij5, d_iij5, d_ijj5, d_iijj5;
+
+	gsl_matrix_set_zero(hess);
+
+	#pragma omp parallel for private(ii, jj)
+	for (unsigned i = 0; i < 3 * n + 2; i++) {
+		if (i != n) {
+			if (i < n) ii = i;
+			if (i > n) ii = i - 1;
+			for (unsigned j = 0; j <= i; j++) {
+				if (j != n) {
+					if (j < n) jj = j;
+					if (j > n) jj = j - 1;
+
+					gsl_matrix_set(hess, ii, jj, gsl_matrix_get(freehess, i, j));
+				}
+			}
+		}
+	}
+
+	indicesRight = gsl_permutation_alloc(n);
+	indicesLeft = gsl_permutation_alloc(n);
+	indices2Right = gsl_permutation_alloc(n);
+
+	gsl_permutation_init(indicesRight);
+	gsl_permutation_init(indicesLeft);
+	gsl_permutation_over(n, indicesRight, true);
+	gsl_permutation_over(n, indicesLeft, false);
+	gsl_permutation_memcpy(indices2Right, indicesRight);
+	gsl_permutation_over(n, indices2Right, true);
+
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) {
+		gsl_matrix_set(hess, 3 * n + 1, i, -1);
+		if (i != 0) gsl_matrix_set(hess, 3 * n + 2, n + i - 1, -1);
+	}
+
+	gsl_vector_free(rx);
+	gsl_vector_free(ry);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_matrix_free(tdots);
+	gsl_matrix_free(dists);
+	gsl_matrix_free(freehess);
+
+	gsl_permutation_free(indicesLeft);
+	gsl_permutation_free(indicesRight);
+	gsl_permutation_free(indices2Right);
+}
+
+
+void domain_energy_fixedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z,
+		double c) {
+
+	domain_energy_fixedHalfHessian(hess, n, z, c);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < 3 * n + 3; i++) {
+		for (unsigned j = 0; j < i; j++) {
+			gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j));
+		}
+	}
+}
+
+
+// The random functions.
+
+double domain_energy_randEnergy(unsigned n, const gsl_vector *z,
+	unsigned ord, const gsl_vector *k, const gsl_vector *a,
+	const gsl_vector *phi) {
+
+	double energy;
+
+	gsl_vector *x, *y;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+
+	unsigned ii;
+
+	gsl_vector *tx, *ty;
+
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	energy = 0;
+
+	for (unsigned i = 0; i < n; i++) {
+		for (unsigned j = 0; j < ord; j++) {
+			energy += gsl_vector_get(a, j) * gsl_sf_sin(gsl_vector_get(k, j) * gsl_vector_get(x, i) + gsl_vector_get(k, j + ord) * gsl_vector_get(y, i) + gsl_vector_get(phi, j)) * (gsl_vector_get(ty, i) / gsl_vector_get(k, j) - gsl_vector_get(tx, i) / gsl_vector_get(k, j + ord)) / 2;
+		}
+	}
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	return energy;
+};
+
+
+void domain_energy_randGradient(gsl_vector *grad, unsigned n,
+	const gsl_vector *z, unsigned ord, const gsl_vector *k,
+	const gsl_vector *a, const gsl_vector *phi) {
+
+	gsl_permutation *indices_right;
+	indices_right = gsl_permutation_alloc(n);
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	gsl_vector *x, *y, *tx, *ty;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+
+	unsigned ii;
+
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	double thesumx, thesumy, aj, kxj, kyj, phij, xi, yi, xii, yii, txi, tyi;
+
+	#pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, thesumx, thesumy, aj, kxj, kyj, phij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		xi = gsl_vector_get(x, i);
+		yi = gsl_vector_get(y, i);
+		xii = gsl_vector_get(x, ii);
+		yii = gsl_vector_get(y, ii);
+		txi = gsl_vector_get(tx, i);
+		tyi = gsl_vector_get(ty, i);
+
+		thesumx = 0;
+		thesumy = 0;
+
+		for (unsigned j = 0; j < ord; j++) {
+			aj = gsl_vector_get(a, j);
+			kxj = gsl_vector_get(k, j);
+			kyj = gsl_vector_get(k, ord + j);
+			phij = gsl_vector_get(phi, j);
+
+			thesumx += aj * (kxj * gsl_sf_cos(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) +
+					(gsl_sf_sin(kxj * xi + kyj * yi + phij) - gsl_sf_sin(kxj * xii + kyj * yii + phij)) / kyj);
+
+			thesumy += aj * (kyj * gsl_sf_cos(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) -
+					(gsl_sf_sin(kxj * xi + kyj * yi + phij) - gsl_sf_sin(kxj * xii + kyj * yii + phij)) / kxj);
+		}
+
+		gsl_vector_set(grad, i, gsl_vector_get(grad, i) + thesumx / 2);
+		gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) + thesumy / 2);
+	}
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_permutation_free(indices_right);
+}
+
+
+void domain_energy_randHalfHessian(gsl_matrix *hess, unsigned n,
+	const gsl_vector *z, unsigned ord, const gsl_vector *k, const gsl_vector *a,
+	const gsl_vector *phi) {
+
+	gsl_permutation *indices_right;
+	indices_right = gsl_permutation_alloc(n);
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	gsl_vector *x, *y, *tx, *ty;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+
+	unsigned ii;
+
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	double thesumx, thesumy, thesumxy, thesumxx, thesumyy, thesumxxy, thesumxyy, aj, kxj, kyj, phij, xi, yi, xii, yii, txi, tyi;
+
+	#pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, thesumx, thesumy, thesumxy, thesumxx, thesumyy, thesumxxy, thesumxyy, aj, kxj, kyj, phij)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		xi = gsl_vector_get(x, i);
+		yi = gsl_vector_get(y, i);
+		xii = gsl_vector_get(x, ii);
+		yii = gsl_vector_get(y, ii);
+		txi = gsl_vector_get(tx, i);
+		tyi = gsl_vector_get(ty, i);
+
+		thesumx = 0;
+		thesumy = 0;
+		thesumxy = 0;
+		thesumxx = 0;
+		thesumyy = 0;
+		thesumxxy = 0;
+		thesumxyy = 0;
+
+		for (unsigned j = 0; j < ord; j++) {
+			aj = gsl_vector_get(a, j);
+			kxj = gsl_vector_get(k, j);
+			kyj = gsl_vector_get(k, ord + j);
+			phij = gsl_vector_get(phi, j);
+
+			thesumx += aj * ( - gsl_pow_2(kxj) * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) +
+					2 * gsl_sf_cos(kxj * xi + kyj * yi + phij) * kxj / kyj);
+
+			thesumy += aj * ( - gsl_pow_2(kyj) * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj) -
+					2 * gsl_sf_cos(kxj * xi + kyj * yi + phij) * kyj / kxj);
+
+			thesumxy += - aj * kxj * kyj * gsl_sf_sin(kxj * xi + kyj * yi + phij) * (tyi / kxj - txi / kyj);
+
+			thesumxx += - aj * gsl_sf_cos(kxj * xii + kyj * yii + phij) * kxj / kyj;
+
+			thesumyy += aj * gsl_sf_cos(kxj * xii + kyj * yii + phij) * kyj / kxj;
+
+			thesumxyy += - aj * gsl_sf_cos(kxj * xii + kyj * yii + phij);
+
+			thesumxxy += aj * gsl_sf_cos(kxj * xii + kyj * yii + phij);
+		}
+
+		gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + thesumx / 2);
+		gsl_matrix_set(hess, n + i, n + i, gsl_matrix_get(hess, n + i, n + i) + thesumy / 2);
+		gsl_matrix_set(hess, i + n, i, gsl_matrix_get(hess, n + i, i) + thesumxy / 2);
+		if (i == 0) {
+			gsl_matrix_set(hess, n - 1, 0, gsl_matrix_get(hess, n - 1, 0) + thesumxx / 2);
+			gsl_matrix_set(hess, 2 * n - 1, n, gsl_matrix_get(hess, 2 * n - 1, n) + thesumyy / 2);
+		} else {
+			gsl_matrix_set(hess, i, ii, gsl_matrix_get(hess, i, ii) + thesumxx / 2);
+			gsl_matrix_set(hess, n + i, n + ii, gsl_matrix_get(hess, n + i, n + ii) + thesumyy / 2);
+		}
+		gsl_matrix_set(hess, n + i, ii, gsl_matrix_get(hess, n + i, ii) + thesumxxy / 2);
+		gsl_matrix_set(hess, n + ii, i, gsl_matrix_get(hess, n + ii, i) + thesumxyy / 2);
+	}
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_permutation_free(indices_right);
+}
+
+
+// The random naked functions.
+
+double domain_energy_nakedRandLagrangian(unsigned n, const gsl_vector *z,
+	double c, unsigned ord,  const gsl_vector *k, const gsl_vector *a,
+	const gsl_vector *phi) {
+
+	double lagrangian, randEnergy;
+
+	lagrangian = domain_energy_nakedLagrangian(n, z, c);
+	randEnergy = domain_energy_randEnergy(n, z, ord, k, a, phi);
+
+	return lagrangian + randEnergy;
+}
+
+
+void domain_energy_nakedRandGradient(gsl_vector *grad, unsigned n,
+	const gsl_vector *z, double c, unsigned ord,  const gsl_vector *k,
+	const gsl_vector *a, const gsl_vector *phi) {
+
+	domain_energy_nakedGradient(grad, n, z, c);
+	domain_energy_randGradient(grad, n, z, ord, k, a, phi);
+}
+
+
+void domain_energy_nakedRandHessian(gsl_matrix *hess, unsigned n,
+	const gsl_vector *z, double c, unsigned ord,  const gsl_vector *k,
+	const gsl_vector *a, const gsl_vector *phi) {
+
+	domain_energy_nakedHalfHessian(hess, n, z, c);
+	domain_energy_randHalfHessian(hess, n, z, ord, k, a, phi);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < 3 * n + 2; i++) {
+		for (unsigned j = 0; j < i; j++) {
+			gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j));
+		}
+	}
+}
+
+
+// The well functions.
+
+double domain_energy_wellEnergy(unsigned n, const gsl_vector *z, double w,
+	double s) {
+
+	double xi, yi, txi, tyi, energy;
+	gsl_vector *x, *y, *tx, *ty;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i,
+		gsl_vector_get(z, i + n));
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	energy = 0;
+
+	for (unsigned i = 0; i < n; i++) {
+		xi = gsl_vector_get(x, i);
+		yi = gsl_vector_get(y, i);
+		txi = gsl_vector_get(tx, i);
+		tyi = gsl_vector_get(ty, i);
+
+		energy += ((gsl_sf_exp(s * (xi - w)) - gsl_sf_exp(-s * (xi + w))) * tyi -
+			(gsl_sf_exp(s * (yi - w)) - gsl_sf_exp(-s * (yi + w))) * txi) / s;
+	}
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+
+	return energy;
+};
+
+
+void domain_energy_wellGradient(gsl_vector *grad, unsigned n,
+	const gsl_vector *z, double w, double s) {
+
+	unsigned ii;
+	double xi, yi, xii, yii, txi, tyi, grad_xi, grad_yi;
+	gsl_vector *x, *y, *tx, *ty;
+	gsl_permutation *indices_right;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+	indices_right = gsl_permutation_alloc(n);
+
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	#pragma omp parallel for
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+
+	#pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, grad_xi,\
+		grad_yi)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		xi = gsl_vector_get(x, i);
+		yi = gsl_vector_get(y, i);
+		xii = gsl_vector_get(x, ii);
+		yii = gsl_vector_get(y, ii);
+		txi = gsl_vector_get(tx, i);
+		tyi = gsl_vector_get(ty, i);
+
+		grad_xi = tyi * (gsl_sf_exp(s * (xi - w)) + gsl_sf_exp(-s * (xi + w))) +
+			(gsl_sf_exp(s * (yi - w)) - gsl_sf_exp(-s * (yi + w)) -
+			 gsl_sf_exp(s * (yii - w)) + gsl_sf_exp(-s * (yii + w))) / s;
+
+		grad_yi = - txi * (gsl_sf_exp(s * (yi - w)) + gsl_sf_exp(-s * (yi + w))) +
+			(- gsl_sf_exp(s * (xi - w)) + gsl_sf_exp(-s * (xi + w)) +
+			 gsl_sf_exp(s * (xii - w)) - gsl_sf_exp(-s * (xii + w))) / s;
+
+		gsl_vector_set(grad, i, gsl_vector_get(grad, i) + grad_xi);
+
+		gsl_vector_set(grad, n + i, gsl_vector_get(grad, n + i) + grad_yi);
+	}
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_permutation_free(indices_right);
+}
+
+
+void domain_energy_wellHalfHessian(gsl_matrix *hess, unsigned n,
+	const gsl_vector *z, double w, double s) {
+
+	unsigned ii;
+	double xi, yi, xii, yii, txi, tyi, hess_xi, hess_yi, hess_xiyi, hess_xiiyi,
+		hess_xiyii, exp_mxi, exp_pxi, exp_myi, exp_pyi, exp_mxii, exp_myii,
+		exp_pxii, exp_pyii;
+	gsl_vector *x, *y, *tx, *ty;
+	gsl_permutation *indices_right;
+
+	x = gsl_vector_alloc(n);
+	y = gsl_vector_alloc(n);
+	tx = gsl_vector_alloc(n);
+	ty = gsl_vector_alloc(n);
+	indices_right = gsl_permutation_alloc(n);
+
+	for (unsigned i = 0; i < n; i++) gsl_vector_set(x, i, gsl_vector_get(z, i));
+	for (unsigned i = 0; i < n; i++) {
+		gsl_vector_set(y, i, gsl_vector_get(z, i + n));
+	}
+
+	domain_energy_rt(tx, n, x, -1);
+	domain_energy_rt(ty, n, y, -1);
+
+	gsl_permutation_init(indices_right);
+	gsl_permutation_over(n, indices_right, true);
+
+	#pragma omp parallel for private(ii, xi, yi, xii, yii, txi, tyi, hess_xi,\
+		hess_yi, hess_xiyi, hess_xiiyi, hess_xiyii, exp_mxi, exp_pxi, exp_myi,\
+		exp_pyi, exp_mxii, exp_myii, exp_pxii, exp_pyii)
+	for (unsigned i = 0; i < n; i++) {
+		ii = gsl_permutation_get(indices_right, i);
+
+		xi = gsl_vector_get(x, i);
+		yi = gsl_vector_get(y, i);
+		xii = gsl_vector_get(x, ii);
+		yii = gsl_vector_get(y, ii);
+		txi = gsl_vector_get(tx, i);
+		tyi = gsl_vector_get(ty, i);
+		exp_mxi = gsl_sf_exp(-s * (xi + w));
+		exp_pxi = gsl_sf_exp(s * (xi - w));
+		exp_myi = gsl_sf_exp(-s * (yi + w));
+		exp_pyi = gsl_sf_exp(s * (yi - w));
+		exp_mxii = gsl_sf_exp(-s * (xii + w));
+		exp_pxii = gsl_sf_exp(s * (xii - w));
+		exp_myii = gsl_sf_exp(-s * (yii + w));
+		exp_pyii = gsl_sf_exp(s * (yii - w));
+
+		hess_xi = tyi * s * (exp_pxi - exp_mxi);
+		hess_yi = - txi * s * (exp_pyi - exp_myi);
+
+		hess_xiyi = - exp_pxi - exp_mxi + exp_pyi + exp_myi;
+
+		hess_xiyii = - exp_pyii - exp_myii;
+		hess_xiiyi = exp_pxii + exp_mxii;
+
+		gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + hess_xi);
+		gsl_matrix_set(hess, n + i, n + i, gsl_matrix_get(hess, n + i, n + i) +
+			hess_yi);
+		gsl_matrix_set(hess, i + n, i, gsl_matrix_get(hess, n + i, i) + hess_xiyi);
+		gsl_matrix_set(hess, n + i, ii, gsl_matrix_get(hess, n + i, ii) +
+			hess_xiiyi);
+		gsl_matrix_set(hess, n + ii, i, gsl_matrix_get(hess, n + ii, i) +
+			hess_xiyii);
+	}
+
+	gsl_vector_free(x);
+	gsl_vector_free(y);
+	gsl_vector_free(tx);
+	gsl_vector_free(ty);
+	gsl_permutation_free(indices_right);
+}
+
+
+// The naked well functions.
+
+double domain_energy_nakedWellLagrangian(unsigned n, const gsl_vector *z,
+	double c, double w, double s) {
+
+	double nakedLagrangian, wellEnergy;
+
+	nakedLagrangian = domain_energy_nakedLagrangian(n, z, c);
+	wellEnergy = domain_energy_wellEnergy(n, z, w, s);
+
+	return nakedLagrangian + wellEnergy;
+}
+
+
+void domain_energy_nakedWellGradient(gsl_vector *grad, unsigned n,
+	const gsl_vector *z, double c, double w, double s) {
+
+	domain_energy_nakedGradient(grad, n, z, c);
+	domain_energy_wellGradient(grad, n, z, w, s);
+}
+
+
+void domain_energy_nakedWellHessian(gsl_matrix *hess, unsigned n,
+	const gsl_vector *z, double c, double w, double s) {
+
+	domain_energy_nakedHalfHessian(hess, n, z, c);
+	domain_energy_wellHalfHessian(hess, n, z, w, s);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < 3 * n + 2; i++) {
+		for (unsigned j = 0; j < i; j++) {
+			gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j));
+		}
+	}
+}
+
+
+// The random well functions.
+
+double domain_energy_randWellLagrangian(unsigned n, const gsl_vector *z,
+	double c, unsigned ord,  const gsl_vector *k, const gsl_vector *a,
+	const gsl_vector *phi, double w, double s) {
+
+	double lagrangian, randEnergy, wellEnergy;
+
+	lagrangian = domain_energy_nakedLagrangian(n, z, c);
+	randEnergy = domain_energy_randEnergy(n, z, ord, k, a, phi);
+	wellEnergy = domain_energy_wellEnergy(n, z, w, s);
+
+	return lagrangian + randEnergy + wellEnergy;
+}
+
+
+void domain_energy_randWellGradient(gsl_vector *grad, unsigned n,
+	const gsl_vector *z, double c, unsigned ord,  const gsl_vector *k,
+	const gsl_vector *a, const gsl_vector *phi, double w, double s) {
+
+	domain_energy_nakedGradient(grad, n, z, c);
+	domain_energy_randGradient(grad, n, z, ord, k, a, phi);
+	domain_energy_wellGradient(grad, n, z, w, s);
+}
+
+
+void domain_energy_randWellHessian(gsl_matrix *hess, unsigned n,
+	const gsl_vector *z, double c, unsigned ord,  const gsl_vector *k,
+	const gsl_vector *a, const gsl_vector *phi, double w, double s) {
+
+	domain_energy_nakedHalfHessian(hess, n, z, c);
+	domain_energy_randHalfHessian(hess, n, z, ord, k, a, phi);
+	domain_energy_wellHalfHessian(hess, n, z, w, s);
+
+	#pragma omp parallel for
+	for (unsigned i = 1; i < 3 * n + 2; i++) {
+		for (unsigned j = 0; j < i; j++) {
+			gsl_matrix_set(hess, j, i, gsl_matrix_get(hess, i, j));
+		}
+	}
+}
diff --git a/src/domain_energy.h b/src/domain_energy.h
new file mode 100644
index 0000000..6cdf882
--- /dev/null
+++ b/src/domain_energy.h
@@ -0,0 +1,41 @@
+#ifndef DOMAIN_ENERGY_H
+#define DOMAIN_ENERGY_H
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+
+double domain_energy_nakedEnergy(unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_nakedLagrangian(unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_nakedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_nakedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_fixedEnergy(unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_fixedLagrangian(unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_fixedGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_fixedHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c);
+
+double domain_energy_nakedRandLagrangian(unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi);
+
+double domain_energy_nakedRandGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi);
+
+double domain_energy_nakedRandHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi);
+
+double domain_energy_nakedWellLagrangian(unsigned n, const gsl_vector *z, double c, double w, double s);
+
+double domain_energy_nakedWellGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, double w, double s);
+
+double domain_energy_nakedWellHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, double w, double s);
+
+double domain_energy_randWellLagrangian(unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s);
+
+double domain_energy_randWellGradient(gsl_vector *grad, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s);
+
+double domain_energy_randWellHessian(gsl_matrix *hess, unsigned n, const gsl_vector *z, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s);
+
+#endif
diff --git a/src/domain_improve.cpp b/src/domain_improve.cpp
new file mode 100644
index 0000000..f69e85f
--- /dev/null
+++ b/src/domain_improve.cpp
@@ -0,0 +1,145 @@
+/* domain_improve.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* A program which facilitates automated mapping of bifurcation points in the
+ * energy of a system where the Hessian is available.  Currently, only a one
+ * dimensional parameter space is supported.
+ */
+
+#include "domain_energy.h"
+#include "domain_minimize.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt, min_fails;
+	unsigned n, N, num;
+	double c, g0, g, eps, energy;
+	char *filename;
+	bool eigenpres = true;
+
+	// Setting default values.
+	eps = 0;
+	num = 25;
+
+	gsl_vector *z, *old_z, *eigenvalues, *trueEigenvalues;
+	gsl_permutation *eigenorder, *trueEigenorder;
+
+	while ((opt = getopt(argc, argv, "n:c:d:g:h:i:N:p:m:j:e:t:s")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'N':
+				N = atoi(optarg);
+				break;
+			case 'g':
+				g0 = atof(optarg);
+				break;
+			case 'i':
+				filename = optarg;
+				break;
+			case 'e':
+				eps = atof(optarg);
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	z = gsl_vector_alloc(3 * n + 3);
+	old_z = gsl_vector_alloc(3 * n + 3);
+	eigenvalues = gsl_vector_alloc(3 * n + 3);
+	trueEigenvalues = gsl_vector_alloc(3 * n + 4);
+	eigenorder = gsl_permutation_alloc(3 * n + 3);
+	trueEigenorder = gsl_permutation_alloc(3 * n + 4);
+
+	g = g0;
+
+	char ch;
+	double throwaway;
+
+	FILE *f = fopen(filename, "r+");
+	while (ch != '\n') ch = fgetc(f);
+	ch = 'a';
+	while (ch != '\n' && ch != '\t') ch = fgetc(f);
+	if (ch == '\n') eigenpres = false;
+
+	rewind(f);
+
+	fscanf(f, "%le\t", &c);
+	fscanf(f, "%le\n", &energy);
+
+	if (eigenpres) {
+		ch = 'a';
+		while (ch != '\n') ch = fgetc(f);
+	}
+	gsl_vector_fscanf(f, z);
+	fclose(f);
+
+	min_fails = domain_minimize_fixed(z, n, c, eps, g, N, 4, 2);
+
+	if (min_fails) {
+		printf("BIFUR: Initial relaxation failed, exiting.\n");
+		return 1;
+	}
+
+//	domain_eigen_values(eigenvalues, n, z, c);
+//	domain_eigen_sort(eigenorder, n, num, eigenvalues);
+
+	energy = domain_energy_fixedEnergy(n, z, c);
+	unsigned ii;
+
+	FILE *newf = fopen(filename, "w");
+	fprintf(newf, "%.12le\t%.12le\n", c, energy);
+	for (unsigned i = 0; i < num; i++) {
+		ii = gsl_permutation_get(eigenorder, i);
+		fprintf(newf, "%.12le\t", gsl_vector_get(eigenvalues, ii));
+	}
+	fprintf(newf, "\n");
+	for (unsigned i = 0; i < num; i++) {
+		ii = gsl_permutation_get(trueEigenorder, i);
+		fprintf(newf, "%.12le\t", gsl_vector_get(trueEigenvalues, ii));
+	}
+	fprintf(newf, "\n");
+	for (unsigned i = 0; i < 3 * n + 3; i++) {
+		fprintf(newf, "%.12le\t", gsl_vector_get(z, i));
+	}
+	fclose(newf);
+}
diff --git a/src/domain_increase.cpp b/src/domain_increase.cpp
new file mode 100644
index 0000000..9c95d8d
--- /dev/null
+++ b/src/domain_increase.cpp
@@ -0,0 +1,120 @@
+/* domain_init.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// Initializes a circular domain, or a circular domain with a perturbation.
+
+#include <unistd.h>
+#include <iostream>
+#include <string>
+#include <stdlib.h>
+#include <math.h>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+
+#include "domain_minimize.h"
+#include "domain_energy.h"
+
+
+
+
+int main(int argc, char *argv[]) {
+
+	// Declaring variables.
+	gsl_vector *z, *new_z;
+	int opt;
+	unsigned n;
+	double c;
+	char out_file[20];
+
+	// Default values.
+	sprintf(out_file, "%s", "out.dat");
+
+	// GNU getopt in action.
+	while ((opt = getopt(argc, argv, "n:c:i:")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'i':
+				sprintf(out_file, "%s", optarg);
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	z = gsl_vector_alloc(3 * n + 3);
+	new_z = gsl_vector_alloc(6 * n + 3);
+
+	FILE *f = fopen(out_file, "r");
+	gsl_vector_fscanf(f, z);
+	fclose(f);
+
+	gsl_vector_set(new_z, 2 * n - 1, (gsl_vector_get(z, 0) + gsl_vector_get(z, n - 1)) / 2);
+		gsl_vector_set(new_z, 2 * n - 2, gsl_vector_get(z, n - 2));
+
+	for (unsigned i = 0; i < n - 1; i ++) {
+		gsl_vector_set(new_z, 2 * i, gsl_vector_get(z, i));
+		gsl_vector_set(new_z, 2 * i + 1,
+				(gsl_vector_get(z, i) + gsl_vector_get(z, i + 1)) / 2);
+	}
+
+	gsl_vector_set(new_z, 2 * n, gsl_vector_get(z, n) / 2);
+	gsl_vector_set(new_z, 2 * n + 1, gsl_vector_get(z, n));
+	gsl_vector_set(new_z, 4 * n - 2, gsl_vector_get(z, 2 * n - 2) / 2);
+	gsl_vector_set(new_z, 4 * n - 3, gsl_vector_get(z, 2 * n - 2));
+
+	for (unsigned i = 1; i < n - 1; i ++) {
+		gsl_vector_set(new_z, 2 * (n + i) + 1, gsl_vector_get(z, n + i));
+		gsl_vector_set(new_z, 2 * (n + i) ,
+				(gsl_vector_get(z, n + i - 1) + gsl_vector_get(z, n + i)) / 2);
+	}
+
+	gsl_vector_set(new_z, 4 * n - 1, gsl_vector_get(z, 2 * n -1));
+
+	gsl_vector_set(new_z, 4 * n, gsl_vector_get(z, 2 * n));
+
+	for (unsigned i = 0; i < n; i++) {
+		gsl_vector_set(new_z, 4 * n + 1 + 2 * i, gsl_vector_get(z, 2 * n + 1 + i));
+		gsl_vector_set(new_z, 4 * n + 2 + 2 * i, gsl_vector_get(z, 2 * n + 1 + i));
+	}
+
+	gsl_vector_set(new_z, 6 * n + 1, gsl_vector_get(z, 3 * n + 1));
+	gsl_vector_set(new_z, 6 * n + 2, gsl_vector_get(z, 3 * n + 2));
+
+	int result;
+	result = domain_minimize_fixed(new_z, 2 * n, c, 1e-8, 0.00001, 5000, 4, 2);
+
+	if (result) printf("Converging failed.");
+	else {
+		FILE *fout = fopen(out_file, "w");
+		gsl_vector_fprintf(fout, new_z, "%.10g");
+		fclose(fout);
+	}
+
+	gsl_vector_free(z);
+	gsl_vector_free(new_z);
+
+}
+
diff --git a/src/domain_minimize.cpp b/src/domain_minimize.cpp
new file mode 100644
index 0000000..f316e4b
--- /dev/null
+++ b/src/domain_minimize.cpp
@@ -0,0 +1,300 @@
+/* domain_minimize.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// A Newton's method solver for modulated domains.
+
+// GSL includes.
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+// Gives the necessary functions for the Lagrangian, gradient, and Hessian.
+#include "domain_energy.h"
+#include "domain_newton.h"
+
+struct nakedgetgrad {
+	nakedgetgrad(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedGradient(grad, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct nakedgethess {
+	nakedgethess(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedHessian(hess, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct nakedgetenergy {
+	nakedgetenergy(double c, unsigned n): c(c), n(n) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedLagrangian(n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+// Carries out Newton's method.
+int domain_minimize_naked(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) {
+
+	unsigned size = 3 * n + 2;
+	unsigned params = 2 * n + 1;
+	nakedgetgrad grad(c, n);
+	nakedgethess hess(c, n);
+	nakedgetenergy energy(c, n);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb);
+}
+
+struct fixedgetgrad {
+	fixedgetgrad(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_fixedGradient(grad, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct fixedgethess {
+	fixedgethess(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_fixedHessian(hess, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct fixedgetenergy {
+	fixedgetenergy(double c, unsigned n): c(c), n(n) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_fixedLagrangian(n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+// Carries out Newton's method.
+int domain_minimize_fixed(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma) {
+
+	unsigned size = 3 * n + 3;
+	unsigned params = 2 * n;
+	fixedgetgrad grad(c, n);
+	fixedgethess hess(c, n);
+	fixedgetenergy energy(c, n);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, 0, 0, 0.1, 10, true);
+}
+
+struct randgetgrad {
+	randgetgrad(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedRandGradient(grad, n, state, c, ord, k, a, phi);}
+
+	private:
+		double c;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+struct randgethess {
+	randgethess(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedRandHessian(hess, n, state, c, ord, k, a, phi);}
+
+	private:
+		double c;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+struct randgetenergy {
+	randgetenergy(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi): c(c), n(n), ord(ord), k(k), a(a), phi(phi) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedRandLagrangian(n, state, c, ord, k, a, phi);}
+
+	private:
+		double c;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+// Carries out Newton's method.
+int domain_minimize_rand(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb) {
+
+	unsigned size = 3 * n + 2;
+	unsigned params = 2 * n + 1;
+	randgetgrad grad(c, n, ord, k, a, phi);
+	randgethess hess(c, n, ord, k, a, phi);
+	randgetenergy energy(c, n, ord, k, a, phi);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, bound, 0.1, 2, verb);
+}
+
+struct nakedwellgetgrad {
+	nakedwellgetgrad(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_nakedWellGradient(grad, n, state, c, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+};
+
+struct nakedwellgethess {
+	nakedwellgethess(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_nakedWellHessian(hess, n, state, c, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+};
+
+struct nakedwellgetenergy {
+	nakedwellgetenergy(double c, unsigned n, double w, double s): c(c), n(n), w(w), s(s) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_nakedWellLagrangian(n, state, c, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+};
+
+// Carries out Newton's method.
+int domain_minimize_nakedWell(gsl_vector *z, unsigned n, double c, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) {
+
+	unsigned size = 3 * n + 2;
+	unsigned params = 2 * n + 1;
+	nakedwellgetgrad grad(c, n, w, ss);
+	nakedwellgethess hess(c, n, w, ss);
+	nakedwellgetenergy energy(c, n, w, ss);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb);
+}
+
+
+struct randwellgetgrad {
+	randwellgetgrad(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_randWellGradient(grad, n, state, c, ord, k, a, phi, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+struct randwellgethess {
+	randwellgethess(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_randWellHessian(hess, n, state, c, ord, k, a, phi, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+struct randwellgetenergy {
+	randwellgetenergy(double c, unsigned n, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double s): c(c), n(n), ord(ord), k(k), a(a), phi(phi), w(w), s(s) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_randWellLagrangian(n, state, c, ord, k, a, phi, w, s);}
+
+	private:
+		double c;
+		double s;
+		double w;
+		unsigned n;
+		unsigned ord;
+		const gsl_vector *k;
+		const gsl_vector *a;
+		const gsl_vector *phi;
+};
+
+// Carries out Newton's method.
+int domain_minimize_randWell(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb) {
+
+	unsigned size = 3 * n + 2;
+	unsigned params = 2 * n + 1;
+	randwellgetgrad grad(c, n, ord, k, a, phi, w, ss);
+	randwellgethess hess(c, n, ord, k, a, phi, w, ss);
+	randwellgetenergy energy(c, n, ord, k, a, phi, w, ss);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, eta0, 0.1, 100, verb);
+}
+
+
+struct fixedmingetgrad {
+	fixedmingetgrad(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_vector* grad, unsigned size, gsl_vector *state) {domain_energy_fixedGradient(grad, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct fixedmingethess {
+	fixedmingethess(double c, unsigned n): c(c), n(n) {}
+	void operator()(gsl_matrix* hess, unsigned size, gsl_vector *state) {domain_energy_fixedHessian(hess, n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+struct fixedmingetenergy {
+	fixedmingetenergy(double c, unsigned n): c(c), n(n) {}
+	double operator()(unsigned size, gsl_vector *state) {return domain_energy_fixedLagrangian(n, state, c);}
+
+	private:
+		double c;
+		unsigned n;
+};
+
+// Carries out Newton's method.
+int domain_minimize_fixedmin(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb) {
+
+	unsigned size = 3 * n + 3;
+	unsigned params = 2 * n;
+	fixedmingetgrad grad(c, n);
+	fixedmingethess hess(c, n);
+	fixedmingetenergy energy(c, n);
+
+	return domain_newton(z, size, params, energy, grad, hess, eps, N, beta, s, sigma, gamma, bound, 0.1, 10, verb);
+}
diff --git a/src/domain_minimize.h b/src/domain_minimize.h
new file mode 100644
index 0000000..3754e6a
--- /dev/null
+++ b/src/domain_minimize.h
@@ -0,0 +1,24 @@
+#ifndef DOMAIN_MINIMIZE_H
+#define DOMAIN_MINIMIZE_H
+
+// GSL includes.
+#include <gsl/gsl_vector.h>
+
+// Gives the necessary functions for the Lagrangian, gradient, and Hessian.
+#include "domain_energy.h"
+#include "domain_newton.h"
+
+int domain_minimize_naked(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double eta0, bool verb);
+
+int domain_minimize_fixed(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma);
+
+int domain_minimize_rand(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb);
+
+int domain_minimize_nakedWell(gsl_vector *z, unsigned n, double c, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb);
+
+int domain_minimize_randWell(gsl_vector *z, unsigned n, double c, unsigned ord, const gsl_vector *k, const gsl_vector *a, const gsl_vector *phi, double w, double ss, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb);
+
+int domain_minimize_fixedmin(gsl_vector *z, unsigned n, double c, double eps, unsigned N, double beta, double s, double sigma, double gamma, double bound, bool verb);
+
+#endif
+
diff --git a/src/domain_newton.h b/src/domain_newton.h
new file mode 100644
index 0000000..abb7ea3
--- /dev/null
+++ b/src/domain_newton.h
@@ -0,0 +1,224 @@
+#ifndef DOMAIN_NEWTON_H
+#define DOMAIN_NEWTON_H
+
+#include <math.h>
+
+// GSL includes.
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_sf.h>
+
+// Eigen's linear solving uses cheap parallelization.
+#include <eigen3/Eigen/Dense>
+
+/* This function is templated so that any set of functions which return an
+ * energy, gradient, and Hessian given an empty object, the size of the state
+ * vector, and the state vector can be used.  This allows many such sets of
+ * functions, e.g., that for a fixed domain or a domain on a random background,
+ * to be used.  See the file domain_minimize.cpp for examples of construction
+ * of these functions.
+ */
+template <class energy_func, class grad_func, class hess_func>
+
+int domain_newton(gsl_vector *state, unsigned size, unsigned params,
+		energy_func get_energy, grad_func get_grad, hess_func get_hess, double
+		epsilon, unsigned max_iterations, double beta, double s, double sigma,
+		double gamma, double eta_0, double delta, double bound, bool verbose) {
+/* The function domain_newton carries out a modified version of Newton's
+ * method.  On success, 0 is returned.  On failure, 1 is returned.
+ * 
+ * state          - GSL_VECTOR
+ *                On entry, state gives the system's initial condition.  On
+ *                exit, state contains the result Newton's method.
+ *
+ * size           - UNSIGNED INTEGER
+ *                On entry, size gives the size of the vector state.  Unchanged
+ *                on exit.
+ *
+ * params         - UNSIGNED INTEGER
+ *                On entry, params gives the number of non-multiplier elements
+ *                in state, which are assumed by the function to be the first
+ *                elements of state.  Unchanged on exit.
+ *
+ * get_energy     - ENERGY_FUNC
+ *                On entry, get_energy is a function that returns a double
+ *                float.  The first argument of get_energy is an unsigned
+ *                integer and the second argument is a gsl_vector object.  This
+ *                function is expected to take size and state, respectively,
+ *                and return the energy of that state.  Unchanged on exit.
+ *
+ * get_grad       - GRAD_FUNC
+ *                On entry, get_grad is a function that returns void.  The
+ *                first argument of get_grad is a gsl_vector object, the second
+ *                argument of get_grad is an unsigned integer, and the third
+ *                argument of get_grad is a gsl_vector object.  This function
+ *                is expected to take a vector of size size, size, and state,
+ *                respectively.  It leaves the gradient of the energy function
+ *                in the first argument.  Unchanged on exit.
+ *
+ * get_hess       - HESS_FUNC
+ *                On entry, get_hess is a function that returns void.  The
+ *                first argument of get_hess is a gsl_matrix object, the second
+ *                argument of get_hess is an unsigned integer, and the third
+ *                argument of get_hess is a gsl_vector object.  This function
+ *                is expected to take a matrix of size size by size, size, and
+ *                state, respectively.  It leaves the Hessian of the energy
+ *                function in the first argument.  Unchanged on exit.
+ *
+ * epsilon        - DOUBLE FLOAT
+ *                On entry, epsilon gives the number that is used to judge
+ *                convergence.  When the norm of the gradient is less than
+ *                epsilon * size, the process is deemed complete and the
+ *                iterations are stopped.  Unchanged on exit.
+ *
+ * max_iterations - UNSIGNED INTEGER
+ *                On entry, max_iterations gives the maximum number of times
+ *                the algorithm will repeat before failing.  Unchanged on exit.
+ *
+ * beta           - DOUBLE FLOAT
+ *                On entry, beta gives the number which is exponentiated to
+ *                scale the step size in Newton's method.  Unchanged on exit.
+ *
+ * s              - DOUBLE FLOAT
+ *                On entry, s gives a constant scaling of the step size in
+ *                Newton's method.  Unchanged on exit.
+ *
+ * sigma          - DOUBLE FLOAT
+ *                On entry, sigma gives a scaling to the condition on the step
+ *                size in Newton's method.  Unchanged on exit.
+ *
+ * gamma          - DOUBLE FLOAT
+ *                On entry, gamma gives the amount by which the norm of the
+ *                gradient must change for eta to decrement by a factor delta.
+ *                Unchanged on exit.
+ *
+ * eta_0          - DOUBLE FLOAT
+ *                On entry, eta_0 gives the starting value of eta.  Unchanged
+ *                on exit.
+ *
+ * delta          - DOUBLE FLOAT
+ *                On entry, delta gives the factor by which eta is decremented.
+ *                Unchanged on exit.
+ *
+ * bound          - DOUBLE FLOAT
+ *                On entry, delta gives an upper bound to the gradient norm.
+ *                If surpassed, the execution is halted and the program returns
+ *                failure.  Unchanged on exit.
+ *
+ * verbose        - BOOLEAN
+ *                On entry, verbose indicates whether verbose output will be
+ *                printed to stdout by this program.  Unchanged on exit.
+ */
+
+	// Declaring variables.
+	double ratio, norm, old_norm, old_energy, energy, grad_dz_prod, alpha, eta;
+	unsigned iterations, m;
+	bool converged, bound_exceeded;
+
+	// Declaring GSL variables.
+	gsl_vector *grad, *dz;
+	gsl_matrix *hess;
+
+	// Allocating memory for GSL objects
+	grad = gsl_vector_alloc(size);
+	dz = gsl_vector_alloc(size);
+	hess = gsl_matrix_alloc(size, size);
+
+	// Declaring Eigen map objects to wrap the GSL ones.
+	Eigen::Map<Eigen::VectorXd> grad_eigen(grad->data, size);
+	Eigen::Map<Eigen::VectorXd> dz_eigen(dz->data, size);
+	Eigen::Map<Eigen::MatrixXd> hess_eigen(hess->data, size, size);
+
+	// If epsilon > 0, use its value.  Otherwise, set to machine precision.
+	if (epsilon == 0) epsilon = DBL_EPSILON;
+
+	// Initializes the starting value of old_norm at effectively infinity.
+	old_norm = 1 / DBL_EPSILON;
+
+	// Start the loop parameter at zero.
+	iterations = 0;
+
+	/* If the loop ends and this boolean has not been flipped, the program will
+	 * know it has not converged.
+	 */
+	converged = false;
+
+	// Initializes the value of eta.
+	eta = eta_0;
+
+	// Begins the algorithm's loop.
+	while (iterations < max_iterations) {
+
+		// Gets the energy, gradient and Hessian for this iteration.
+		old_energy = get_energy(size, state);
+		get_grad(grad, size, state);
+		get_hess(hess, size, state);
+
+		// Adds eta along the diagonal of the Hessian for non-multiplier entries.
+		for (unsigned i = 0; i < params; i++) {
+			gsl_matrix_set(hess, i, i, gsl_matrix_get(hess, i, i) + eta);
+		}
+
+		// Use LU decomposition to solve for the next step in Newton's method.
+		dz_eigen = hess_eigen.lu().solve(grad_eigen);
+
+		// Dots the gradient into the step in order to judge the step size.
+		gsl_blas_ddot(grad, dz, &grad_dz_prod);
+
+		// Initializes the Armijo counter.
+		m = 0;
+
+		// This loop determines the Armijo step size.
+		while (true) {
+			alpha = gsl_sf_pow_int(beta, m) * s;
+			gsl_vector_scale(dz, alpha);
+			gsl_vector_sub(state, dz);
+
+			energy = get_energy(size, state);
+
+			if (fabs(old_energy - energy) >= sigma * alpha * grad_dz_prod) break;
+			else {
+				gsl_vector_add(state, dz);
+				gsl_vector_scale(dz, 1 / alpha);
+				m++;
+			}
+		}
+
+		// Gets the new norm of the gradient for comparison.
+		norm = gsl_blas_dnrm2(grad) / size;
+
+		// Judges if the norm has changed sufficiently little to decrement eta.
+		if (fabs(norm - old_norm) < gamma * eta) eta *= delta;
+
+		// Prints several useful statistics for debugging purposes.
+		if (verbose) printf("m was %i, grad_norm was %e, eta was %e, energy was %e\n",
+				m, norm, eta, energy);
+
+		// Determines if the process has converged to acceptable precision.
+		if (norm < epsilon) {
+			converged = true;
+			break;
+		}
+
+		// Causes the program to fail if norm has diverged to a large number.
+		if (norm > bound) break;
+
+		// Reset the norm for the next iteration.
+		old_norm = norm;
+
+		// Increment the counter.
+		iterations++;
+	}
+
+	// Gotta live free, die hard.  No one likes memory leaks.
+	gsl_vector_free(grad);
+	gsl_vector_free(dz);
+	gsl_matrix_free(hess);
+
+	// Return conditions to indicate success or failure.
+	if (converged) return 0;
+	else return 1;
+}
+
+#endif
diff --git a/src/eigenvalues.cpp b/src/eigenvalues.cpp
new file mode 100644
index 0000000..c812d0f
--- /dev/null
+++ b/src/eigenvalues.cpp
@@ -0,0 +1,156 @@
+/* eigenvalues.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* This program allows for the generalized eigenvalues of a modulated domain to
+ * be computed and returned.
+ */
+
+#include "domain_energy.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt, min_fails, eigen_follow, eigen_num, examining;
+	unsigned n, N, ord, size, params, j, M;
+	double d, c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, eps2, state, old_state, h, bound;
+	char *in_filename, *out_filename, *k_filename, *a_filename, *phi_filename, str[19], in;
+	bool subcrit, reset, rand, verbose, fixed;
+
+	// Setting default values.
+
+	gsl_vector *z, *k, *a, *phi, *old_z, *eigenvalues;
+	gsl_permutation *eigenorder;
+	gsl_matrix *hess;
+	rand = false;
+	fixed = false;
+	verbose = false;
+	N=25;
+
+	while ((opt = getopt(argc, argv, "n:c:i:o:O:K:A:P:e:g:N:b:rvd:M:f")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'i':
+				in_filename = optarg;
+				break;
+			case 'O':
+				ord = atoi(optarg);
+				break;
+			case 'K':
+				k_filename = optarg;
+				break;
+			case 'A':
+				a_filename = optarg;
+				break;
+			case 'P':
+				phi_filename = optarg;
+				break;
+			case 'N':
+				N = atoi(optarg);
+				break;
+			case 'r':
+				rand = true;
+				break;
+			case 'f':
+				fixed = true;
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (rand || !fixed) {
+		size = 3 * n + 2;
+		params = 2 * n + 1;
+	} else {
+		size = 3 * n + 3;
+		params = 2 * n;
+	}
+
+	z = gsl_vector_alloc(size);
+	eigenvalues = gsl_vector_alloc(size);
+	eigenorder = gsl_permutation_alloc(size);
+	hess = gsl_matrix_alloc(size, size);
+	if (rand) {
+		k = gsl_vector_alloc(2 * ord);
+		a = gsl_vector_alloc(ord);
+		phi = gsl_vector_alloc(ord);
+	}
+
+	FILE *in_file = fopen(in_filename, "r");
+	gsl_vector_fscanf(in_file, z);
+	fclose(in_file);
+
+	if (rand) {
+		FILE *k_file = fopen(k_filename, "r");
+		gsl_vector_fscanf(k_file, k);
+		fclose(k_file);
+
+		FILE *a_file = fopen(a_filename, "r");
+		gsl_vector_fscanf(a_file, a);
+		fclose(a_file);
+
+		FILE *phi_file = fopen(phi_filename, "r");
+		gsl_vector_fscanf(phi_file, phi);
+		fclose(phi_file);
+	}
+
+	if (rand) domain_energy_nakedRandHessian(hess, n, z, c, ord, k, a, phi);
+	else {
+		if (fixed) domain_energy_fixedHessian(hess, n, z, c);
+		else domain_energy_nakedHessian(hess, n, z, c);
+	}
+
+	domain_eigen_values(eigenvalues, size, params, hess);
+	domain_eigen_sort(eigenorder, size, N, eigenvalues);
+
+	for (unsigned i = 0; i < N; i++) {
+		printf("%e\t", gsl_vector_get(eigenvalues, gsl_permutation_get(eigenorder, i)));
+	}
+	printf("\n");
+
+	gsl_vector_free(z);
+
+
+}
diff --git a/src/evolve.cpp b/src/evolve.cpp
new file mode 100644
index 0000000..f23e15e
--- /dev/null
+++ b/src/evolve.cpp
@@ -0,0 +1,237 @@
+/* evolve.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* This program allows for the iterated minimization of modulated domains as
+ * the dimensionless parameter Lambda is varied.
+ */
+
+#include "domain_energy.h"
+#include "domain_minimize.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt, min_fails, eigen_follow, eigen_num, examining;
+	unsigned n, N, ord, size, params, j, M;
+	double d, c, dc0, dc, g0, g, eigen_thres, approach_thres, eps, eps2, state, old_state, h, bound, da, w, ss;
+	char *in_filename, *out_filename, *k_filename, *a_filename, *phi_filename, str[19], in;
+	bool subcrit, reset, rand, verbose, fixed, well;
+
+	// Setting default values.
+
+	gsl_vector *z, *k, *a, *phi, *old_z;
+	rand = false;
+	fixed = false;
+	well = false;
+	verbose = false;
+	j=0;
+	ss=1;
+
+	while ((opt = getopt(argc, argv, "n:c:i:o:O:K:A:P:e:g:N:b:rvd:M:a:fws:W:j:")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'b':
+				bound = atof(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'i':
+				in_filename = optarg;
+				break;
+			case 'o':
+				out_filename = optarg;
+				break;
+			case 'O':
+				ord = atoi(optarg);
+				break;
+			case 'K':
+				k_filename = optarg;
+				break;
+			case 'A':
+				a_filename = optarg;
+				break;
+			case 'P':
+				phi_filename = optarg;
+				break;
+			case 'g':
+				g0 = atof(optarg);
+				break;
+			case 'N':
+				N = atoi(optarg);
+				break;
+			case 'j':
+				j = atoi(optarg);
+				break;
+			case 'M':
+				M = atoi(optarg);
+				break;
+			case 'e':
+				eps = atof(optarg);
+				break;
+			case 'd':
+				dc0 = atof(optarg);
+				break;
+			case 'a':
+				da = atof(optarg);
+				break;
+			case 'r':
+				rand = true;
+				break;
+			case 'f':
+				fixed = true;
+				break;
+			case 'w':
+				well = true;
+				break;
+			case 'W':
+				w = atof(optarg);
+				break;
+			case 's':
+				ss = atof(optarg);
+				break;
+			case 'v':
+				verbose = true;
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (rand || !fixed) {
+		size = 3 * n + 2;
+		params = 2 * n + 1;
+	} else {
+		size = 3 * n + 3;
+		params = 2 * n;
+	}
+
+	z = gsl_vector_alloc(size);
+	old_z = gsl_vector_alloc(size);
+	if (rand) {
+		k = gsl_vector_alloc(2 * ord);
+		a = gsl_vector_alloc(ord);
+		phi = gsl_vector_alloc(ord);
+	}
+
+	FILE *in_file = fopen(in_filename, "r");
+	gsl_vector_fscanf(in_file, z);
+	fclose(in_file);
+
+	if (rand) {
+		FILE *k_file = fopen(k_filename, "r");
+		gsl_vector_fscanf(k_file, k);
+		fclose(k_file);
+
+		FILE *a_file = fopen(a_filename, "r");
+		gsl_vector_fscanf(a_file, a);
+		fclose(a_file);
+
+		FILE *phi_file = fopen(phi_filename, "r");
+		gsl_vector_fscanf(phi_file, phi);
+		fclose(phi_file);
+	}
+
+	g = g0;
+	dc = dc0;
+
+	double beta = 0.9;
+	double s = 1;
+	double sigma = 0.5;
+
+	if (rand && well) min_fails = domain_minimize_randWell(z, n, c, ord, k, a, phi, w, ss, eps, N, beta, s, sigma, g, bound, verbose);
+	else if (rand) min_fails = domain_minimize_rand(z, n, c, ord, k, a, phi, eps, N, beta, s, sigma, g, bound, verbose);
+	else {
+		if (fixed) min_fails = domain_minimize_fixedmin(z, n, c, eps, N, beta, ss, sigma, g, bound, verbose);
+		else {
+			if (well) min_fails = domain_minimize_nakedWell(z, n, c, w, ss, eps, N, beta, s, sigma, g, bound, verbose);
+			else min_fails = domain_minimize_naked(z, n, c, eps, N, beta, s, sigma, g, bound, verbose);
+		}
+	}
+
+	if (min_fails) {
+		printf("BIFUR: Initial relaxation failed, exiting.\n");
+		FILE *out_file = fopen(out_filename, "w");
+		gsl_vector_fprintf(out_file, z, "%.10e");
+		fclose(out_file);
+		return 1;
+	}
+
+
+	while (j < M) {
+		j += 1;
+		c += dc;
+		g = g0;
+		if (rand) gsl_vector_scale(a, da);
+
+		gsl_vector_memcpy(old_z, z);
+
+		printf("EVOLVE: Step %05d, starting with c = %f.\n", j, c);
+
+		while (true) {
+			if (rand && well) min_fails = domain_minimize_randWell(z, n, c, ord, k, a, phi, w, ss, eps, N, beta, s, sigma, g, bound, verbose);
+			else if (rand) min_fails = domain_minimize_rand(z, n, c, ord, k, a, phi, eps, N, beta, s, sigma, g, bound, verbose);
+			else if (fixed) min_fails = domain_minimize_fixedmin(z, n, c, eps, N, beta, s, sigma, g, bound, verbose);
+			else if (well) min_fails = domain_minimize_nakedWell(z, n, c, w, ss, eps, N, beta, s, sigma, g, bound, verbose);
+			else min_fails = domain_minimize_naked(z, n, c, eps, N, beta, s, sigma, g, bound, verbose);
+
+			if (!min_fails) break;
+			printf("EVOLVE: Newton's method failed to converge, reducing gamma.\n");
+			gsl_vector_memcpy(z, old_z);
+			g *= 0.1;
+		}
+
+		sprintf(str, "output/out-%05d.dat", j);
+		FILE *fout = fopen(str, "w");
+		fprintf(fout, "%.10e\n", c);
+		gsl_vector_fprintf(fout, z, "%.10e");
+		fclose(fout);
+	}
+
+	FILE *out_file = fopen(out_filename, "w");
+	gsl_vector_fprintf(out_file, z, "%.10e");
+	fclose(out_file);
+
+	gsl_vector_free(z);
+
+	return 0;
+
+}
diff --git a/src/gradget.cpp b/src/gradget.cpp
new file mode 100644
index 0000000..2f8128f
--- /dev/null
+++ b/src/gradget.cpp
@@ -0,0 +1,107 @@
+/* eigenget.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// A utility which returns a given generalized eigenvector to the user.
+
+#include "domain_energy.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt;
+	unsigned n, m, nn;
+	double c, d;
+	char *in_filename, *out_filename;
+	bool truehessian;
+
+	gsl_vector *z, *grad, *testgrad;
+
+	truehessian = false;
+
+	while ((opt = getopt(argc, argv, "n:m:c:i:o:td:")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'm':
+				m = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'd':
+				d = atof(optarg);
+				break;
+			case 'i':
+				in_filename = optarg;
+				break;
+			case 'o':
+				out_filename = optarg;
+				break;
+			case 't':
+				truehessian = true;
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (truehessian) nn = 3 * n + 3;
+	else nn = 3 * n + 3;
+
+	z = gsl_vector_alloc(3 * n + 3);
+	grad = gsl_vector_alloc(nn);
+	testgrad = gsl_vector_alloc(nn);
+
+	FILE *in_file = fopen(in_filename, "r");
+	gsl_vector_fscanf(in_file, z);
+	fclose(in_file);
+
+	if (truehessian) domain_energy_fixedGradient(grad, n, z, c);
+	else domain_energy_fixedGradient(grad, n, z, c);
+	domain_energy_fixedGradient(testgrad, n, z, c);
+
+	FILE *out_file = fopen(out_filename, "w");
+	for (unsigned i = 0; i < nn; i++) {
+		fprintf(out_file, "%e\t", gsl_vector_get(grad, i) - gsl_vector_get(testgrad, i));
+	}
+	fclose(out_file);
+}
+
+
diff --git a/src/hessget.cpp b/src/hessget.cpp
new file mode 100644
index 0000000..ad4ae2d
--- /dev/null
+++ b/src/hessget.cpp
@@ -0,0 +1,107 @@
+/* eigenget.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// A utility which returns a given generalized eigenvector to the user.
+
+#include "domain_energy.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt;
+	unsigned n, m, nn;
+	double c;
+	char *in_filename, *out_filename;
+	bool truehessian;
+
+	gsl_vector *z;
+	gsl_matrix *hess, *testhess;
+
+	truehessian = false;
+
+	while ((opt = getopt(argc, argv, "n:m:c:i:o:t")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'm':
+				m = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'i':
+				in_filename = optarg;
+				break;
+			case 'o':
+				out_filename = optarg;
+				break;
+			case 't':
+				truehessian = true;
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (truehessian) nn = 3 * n + 4;
+	else nn = 3 * n + 3;
+
+	z = gsl_vector_alloc(3 * n + 3);
+	hess = gsl_matrix_alloc(nn, nn);
+	testhess = gsl_matrix_alloc(nn, nn);
+
+	FILE *in_file = fopen(in_filename, "r");
+	gsl_vector_fscanf(in_file, z);
+	fclose(in_file);
+
+	domain_energy_fixedHessian(hess, n, z, c);
+	domain_energy_fixedHessian(testhess, n, z, c);
+
+	FILE *out_file = fopen(out_filename, "w");
+	for (unsigned i = 0; i < nn; i++) {
+		for (unsigned j = 0; j < nn; j++) {
+			fprintf(out_file, "%e\t", gsl_matrix_get(hess, i, j) - gsl_matrix_get(testhess, i, j));
+		}
+		fprintf(out_file, "\n");
+	}
+	fclose(out_file);
+}
+
+
diff --git a/src/initialize.cpp b/src/initialize.cpp
new file mode 100644
index 0000000..3c77fc2
--- /dev/null
+++ b/src/initialize.cpp
@@ -0,0 +1,201 @@
+/* initialize.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// Initializes modulated domains.
+
+#include <unistd.h>
+#include <iostream>
+#include <string>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_vector.h>
+
+
+int main(int argc, char *argv[]) {
+
+	// Declaring variables.
+	gsl_vector *z, *k, *a, *phi;
+	int opt;
+	unsigned n, m, size, params, word, ord;
+	double p, k0, a0, R, w0, slope;
+	char *out_filename, *k_filename, *a_filename, *phi_filename;
+	bool rand;
+
+	// Default values.
+	n = 100;
+	p = 0;
+	m = 0;
+	rand = false;
+	k0 = 1;
+	a0 = 1;
+	ord = 0;
+	word = 0;
+	w0=1;
+	slope=1;
+
+	// GNU getopt in action.
+	while ((opt = getopt(argc, argv, "n:p:m:o:O:rK:A:P:k:a:wW:u:s:")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'p':
+				p = atof(optarg);
+				break;
+			case 'a':
+				a0 = atof(optarg);
+				break;
+			case 'k':
+				k0 = atof(optarg);
+				break;
+			case 'm':
+				m = atoi(optarg);
+				break;
+			case 'O':
+				ord = atoi(optarg);
+				break;
+			case 'o':
+				out_filename = optarg;
+				break;
+			case 'K':
+				k_filename = optarg;
+				break;
+			case 'A':
+				a_filename = optarg;
+				break;
+			case 'P':
+				phi_filename = optarg;
+				break;
+			case 'r':
+				rand = true;
+				break;
+			case 'W':
+				word = atoi(optarg);
+				break;
+			case 'u':
+				w0 = atof(optarg);
+				break;
+			case 's':
+				slope = atof(optarg);
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (rand) {
+		size = 3 * n + 2;
+		params = 2 * n + 1;
+	} else {
+		size = 3 * n + 3;
+		params = 2 * n;
+	}
+
+	z = gsl_vector_alloc(size);
+	if (rand) {
+		k = gsl_vector_alloc(2 * (ord + 2 * word));
+		a = gsl_vector_alloc(ord + 2 * word);
+		phi = gsl_vector_alloc(ord + 2* word);
+	}
+
+	R = sqrt(2 * M_PI / (n * sin(2 * M_PI / n)));
+
+	// setting x[0..n], y[1..n]
+	for (unsigned i = 0; i < n; i++) {
+		gsl_vector_set(z, i, R * cos(2 * M_PI * i / n));
+		if (rand) gsl_vector_set(z, n + i, R * sin(2 * M_PI * i / n));
+		else if (i != 0) gsl_vector_set(z, n + i - 1, R * sin(2 * M_PI * i / n));
+	}
+
+	// setting L
+	gsl_vector_set(z, params - 1, 2 * R * n * sin(M_PI / n));
+
+	for (unsigned i = 0; i < size - params; i++) {
+		gsl_vector_set(z, params + i, 1);
+	}
+
+	if (p > 0) {
+		for (unsigned i = 0; i < n; i++) {
+			gsl_vector_set(z, i, gsl_vector_get(z, i) *
+					(1 + p * cos(2 * M_PI * m * i / n)));
+
+			if (rand) gsl_vector_set(z, n + i, gsl_vector_get(z, n + i) *
+					(1 + p * cos(2 * M_PI * m * i / n)));
+			else if (i != 0) gsl_vector_set(z, n + i - 1, gsl_vector_get(z, n + i - 1) *
+					(1 + p * cos(2 * M_PI * m * i / n)));
+		}
+	}
+	if (rand) {
+		// The GSL random number generator.  Don't try to make more than one random
+		// background in the same second, or you'll get identical results.
+		gsl_rng *rand;
+		rand = gsl_rng_alloc(gsl_rng_ranlux);
+		gsl_rng_set(rand, time(NULL));
+
+		double kx, ky, kk;
+
+		for (unsigned i = 0; i < ord; i++) {
+			while (true) {
+				kx = gsl_rng_uniform(rand) * 2 - 1;
+				ky = gsl_rng_uniform(rand) * 2 - 1;
+				kk = gsl_pow_2(kx)+gsl_pow_2(ky);
+
+				if (kk < 1) {
+					gsl_vector_set(k, i, k0 * kx);
+					gsl_vector_set(k, i + ord + 2 * word, k0 * ky);
+					break;
+				}
+			}
+			gsl_vector_set(a, i, 2 * a0 / ord * gsl_rng_uniform(rand));
+			gsl_vector_set(phi, i, 2 * M_PI * gsl_rng_uniform(rand));
+		}
+	}
+
+	if (rand) {
+		FILE *k_file = fopen(k_filename, "w");
+		gsl_vector_fprintf(k_file, k, "%.10g");
+		fclose(k_file);
+
+		FILE *a_file = fopen(a_filename, "w");
+		gsl_vector_fprintf(a_file, a, "%.10g");
+		fclose(a_file);
+
+		FILE *phi_file = fopen(phi_filename, "w");
+		gsl_vector_fprintf(phi_file, phi, "%.10g");
+		fclose(phi_file);
+	}
+
+	FILE *out_file = fopen(out_filename, "w");
+	gsl_vector_fprintf(out_file, z, "%.10g");
+	fclose(out_file);
+
+
+	gsl_vector_free(z);
+	if (rand) {
+		gsl_vector_free(phi);
+		gsl_vector_free(k);
+		gsl_vector_free(a);
+	}
+
+}
+
diff --git a/src/perturb.cpp b/src/perturb.cpp
new file mode 100644
index 0000000..9eaf434
--- /dev/null
+++ b/src/perturb.cpp
@@ -0,0 +1,113 @@
+/* eigenget.cpp
+ *
+ * Copyright (C) 2013 Jaron Kent-Dobias
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// A utility which returns a given generalized eigenvector to the user.
+
+#include "domain_energy.h"
+#include "domain_min.h"
+#include "domain_eigen.h"
+
+#include <unistd.h>
+#include <stdio.h>
+#include <iostream>
+#include <stdlib.h>
+#include <math.h>
+#include <string>
+
+// GSL includes.
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_permute_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_sort_vector.h>
+
+
+// Initializes the program.
+int main(int argc, char *argv[]) {
+
+	int opt;
+	unsigned n, m, nn;
+	double c;
+	char *in_filename, *out_filename;
+	bool truehessian;
+
+	truehessian = false;
+
+	gsl_vector *z, *eigenvalues, *eigenvector;
+	gsl_permutation *eigenorder;
+
+	while ((opt = getopt(argc, argv, "n:m:c:i:o:t")) != -1) {
+		switch (opt) {
+			case 'n':
+				n = atoi(optarg);
+				break;
+			case 'm':
+				m = atoi(optarg);
+				break;
+			case 'c':
+				c = atof(optarg);
+				break;
+			case 'i':
+				in_filename = optarg;
+				break;
+			case 'o':
+				out_filename = optarg;
+				break;
+			case 't':
+				truehessian = true;
+				break;
+			default:
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	if (truehessian) nn = 3 * n + 4;
+	else nn = 3 * n + 3;
+
+	z = gsl_vector_alloc(3 * n + 3);
+	eigenvalues = gsl_vector_alloc(nn);
+	eigenvector = gsl_vector_alloc(nn);
+	eigenorder = gsl_permutation_alloc(nn);
+
+	FILE *in_file = fopen(in_filename, "r");
+	gsl_vector_fscanf(in_file, z);
+	fclose(in_file);
+
+	if (truehessian) {
+		domain_eigen_truevalues(eigenvalues, n, z, c);
+		domain_eigen_truesort(eigenorder, n, m, eigenvalues);
+		domain_eigen_truevector(eigenvector, gsl_permutation_get(eigenorder, m), n, z, c);
+	} else {
+		domain_eigen_values(eigenvalues, n, z, c);
+		domain_eigen_sort(eigenorder, n, m, eigenvalues);
+		domain_eigen_vector(eigenvector, gsl_permutation_get(eigenorder, m), n, z, c);
+	}
+
+	FILE *out_file = fopen(out_filename, "w");
+	for (unsigned i = 0; i < nn; i++) {
+		fprintf(out_file, "%e\t", gsl_vector_get(eigenvector, i));
+	}
+	fclose(out_file);
+}
+
+