From a79c2f04682c44275d2415d39e6996802c4a83c4 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 1 May 2014 16:01:48 -0700 Subject: created a git repository for my thesis code. --- src/domain_eigen.cpp | 187 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 187 insertions(+) create mode 100644 src/domain_eigen.cpp (limited to 'src/domain_eigen.cpp') 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 . + */ + +/* A set of utilities for find the generalized eigenvalues and eigenvectors of + * modulated domains. + */ + +// GSL includes. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +// 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; +} + + -- cgit v1.2.3-54-g00ecf