diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/bifurchaser.cpp | 359 | ||||
| -rw-r--r-- | src/domain_check.cpp | 274 | ||||
| -rw-r--r-- | src/domain_eigen.cpp | 187 | ||||
| -rw-r--r-- | src/domain_eigen.h | 28 | ||||
| -rw-r--r-- | src/domain_energy.cpp | 2103 | ||||
| -rw-r--r-- | src/domain_energy.h | 41 | ||||
| -rw-r--r-- | src/domain_improve.cpp | 145 | ||||
| -rw-r--r-- | src/domain_increase.cpp | 120 | ||||
| -rw-r--r-- | src/domain_minimize.cpp | 300 | ||||
| -rw-r--r-- | src/domain_minimize.h | 24 | ||||
| -rw-r--r-- | src/domain_newton.h | 224 | ||||
| -rw-r--r-- | src/eigenvalues.cpp | 156 | ||||
| -rw-r--r-- | src/evolve.cpp | 237 | ||||
| -rw-r--r-- | src/gradget.cpp | 107 | ||||
| -rw-r--r-- | src/hessget.cpp | 107 | ||||
| -rw-r--r-- | src/initialize.cpp | 201 | ||||
| -rw-r--r-- | src/perturb.cpp | 113 |
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); +} + + |