summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jkentdobias@g.hmc.edu>2014-05-01 16:01:48 -0700
committerJaron Kent-Dobias <jkentdobias@g.hmc.edu>2014-05-01 16:01:48 -0700
commita79c2f04682c44275d2415d39e6996802c4a83c4 (patch)
treed811bfc2f9e781c6b3d262c758223179f9b4f511 /src
downloadcode-a79c2f04682c44275d2415d39e6996802c4a83c4.tar.gz
code-a79c2f04682c44275d2415d39e6996802c4a83c4.tar.bz2
code-a79c2f04682c44275d2415d39e6996802c4a83c4.zip
created a git repository for my thesis code.
Diffstat (limited to 'src')
-rw-r--r--src/bifurchaser.cpp359
-rw-r--r--src/domain_check.cpp274
-rw-r--r--src/domain_eigen.cpp187
-rw-r--r--src/domain_eigen.h28
-rw-r--r--src/domain_energy.cpp2103
-rw-r--r--src/domain_energy.h41
-rw-r--r--src/domain_improve.cpp145
-rw-r--r--src/domain_increase.cpp120
-rw-r--r--src/domain_minimize.cpp300
-rw-r--r--src/domain_minimize.h24
-rw-r--r--src/domain_newton.h224
-rw-r--r--src/eigenvalues.cpp156
-rw-r--r--src/evolve.cpp237
-rw-r--r--src/gradget.cpp107
-rw-r--r--src/hessget.cpp107
-rw-r--r--src/initialize.cpp201
-rw-r--r--src/perturb.cpp113
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);
+}
+
+