summaryrefslogtreecommitdiff
path: root/src/domain_eigen.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/domain_eigen.cpp')
-rw-r--r--src/domain_eigen.cpp187
1 files changed, 187 insertions, 0 deletions
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;
+}
+
+