diff options
Diffstat (limited to 'gsl-1.9/eigen/jacobi.c')
| -rw-r--r-- | gsl-1.9/eigen/jacobi.c | 259 |
1 files changed, 259 insertions, 0 deletions
diff --git a/gsl-1.9/eigen/jacobi.c b/gsl-1.9/eigen/jacobi.c new file mode 100644 index 0000000..7f4a4dc --- /dev/null +++ b/gsl-1.9/eigen/jacobi.c @@ -0,0 +1,259 @@ +/* eigen/jacobi.c + * + * Copyright (C) 2004 Brian Gough, Gerard Jungman + * + * 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 2 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, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_eigen.h> + +/* Algorithm 8.4.3 - Cyclic Jacobi. Golub & Van Loan, Matrix Computations */ + +static inline double +symschur2 (gsl_matrix * A, size_t p, size_t q, double *c, double *s) +{ + double Apq = gsl_matrix_get (A, p, q); + + if (Apq != 0.0) + { + double App = gsl_matrix_get (A, p, p); + double Aqq = gsl_matrix_get (A, q, q); + double tau = (Aqq - App) / (2.0 * Apq); + double t, c1; + + if (tau >= 0.0) + { + t = 1.0 / (tau + hypot (1.0, tau)); + } + else + { + t = -1.0 / (-tau + hypot (1.0, tau)); + } + + c1 = 1.0 / hypot (1.0, t); + + *c = c1; + *s = t * c1; + } + else + { + *c = 1.0; + *s = 0.0; + } + + /* reduction in off(A) is 2*(A_pq)^2 */ + + return fabs (Apq); +} + +inline static void +apply_jacobi_L (gsl_matrix * A, size_t p, size_t q, double c, double s) +{ + size_t j; + const size_t N = A->size2; + + /* Apply rotation to matrix A, A' = J^T A */ + + for (j = 0; j < N; j++) + { + double Apj = gsl_matrix_get (A, p, j); + double Aqj = gsl_matrix_get (A, q, j); + gsl_matrix_set (A, p, j, Apj * c - Aqj * s); + gsl_matrix_set (A, q, j, Apj * s + Aqj * c); + } +} + +inline static void +apply_jacobi_R (gsl_matrix * A, size_t p, size_t q, double c, double s) +{ + size_t i; + const size_t M = A->size1; + + /* Apply rotation to matrix A, A' = A J */ + + for (i = 0; i < M; i++) + { + double Aip = gsl_matrix_get (A, i, p); + double Aiq = gsl_matrix_get (A, i, q); + gsl_matrix_set (A, i, p, Aip * c - Aiq * s); + gsl_matrix_set (A, i, q, Aip * s + Aiq * c); + } +} + +inline static double +norm (gsl_matrix * A) +{ + size_t i, j, M = A->size1, N = A->size2; + double sum = 0.0, scale = 0.0, ssq = 1.0; + + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + double Aij = gsl_matrix_get (A, i, j); + + if (Aij != 0.0) + { + double ax = fabs (Aij); + + if (scale < ax) + { + ssq = 1.0 + ssq * (scale / ax) * (scale / ax); + scale = ax; + } + else + { + ssq += (ax / scale) * (ax / scale); + } + } + + } + } + + sum = scale * sqrt (ssq); + + return sum; +} + +int +gsl_eigen_jacobi (gsl_matrix * a, + gsl_vector * eval, + gsl_matrix * evec, unsigned int max_rot, unsigned int *nrot) +{ + size_t i, p, q; + const size_t M = a->size1, N = a->size2; + double red, redsum = 0.0; + + if (M != N) + { + GSL_ERROR ("eigenproblem requires square matrix", GSL_ENOTSQR); + } + else if (M != evec->size1 || M != evec->size2) + { + GSL_ERROR ("eigenvector matrix must match input matrix", GSL_EBADLEN); + } + else if (M != eval->size) + { + GSL_ERROR ("eigenvalue vector must match input matrix", GSL_EBADLEN); + } + + gsl_vector_set_zero (eval); + gsl_matrix_set_identity (evec); + + for (i = 0; i < max_rot; i++) + { + double nrm = norm (a); + + if (nrm == 0.0) + break; + + for (p = 0; p < N; p++) + { + for (q = p + 1; q < N; q++) + { + double c, s; + + red = symschur2 (a, p, q, &c, &s); + redsum += red; + + /* Compute A <- J^T A J */ + apply_jacobi_L (a, p, q, c, s); + apply_jacobi_R (a, p, q, c, s); + + /* Compute V <- V J */ + apply_jacobi_R (evec, p, q, c, s); + } + } + } + + *nrot = i; + + for (p = 0; p < N; p++) + { + double ep = gsl_matrix_get (a, p, p); + gsl_vector_set (eval, p, ep); + } + + if (i == max_rot) + { + return GSL_EMAXITER; + } + + return GSL_SUCCESS; +} + +int +gsl_eigen_invert_jacobi (const gsl_matrix * a, + gsl_matrix * ainv, unsigned int max_rot) +{ + if (a->size1 != a->size2 || ainv->size1 != ainv->size2) + { + GSL_ERROR("jacobi method requires square matrix", GSL_ENOTSQR); + } + else if (a->size1 != ainv->size2) + { + GSL_ERROR ("inverse matrix must match input matrix", GSL_EBADLEN); + } + + { + const size_t n = a->size2; + size_t i,j,k; + unsigned int nrot = 0; + int status; + + gsl_vector * eval = gsl_vector_alloc(n); + gsl_matrix * evec = gsl_matrix_alloc(n, n); + gsl_matrix * tmp = gsl_matrix_alloc(n, n); + + gsl_matrix_memcpy (tmp, a); + + status = gsl_eigen_jacobi(tmp, eval, evec, max_rot, &nrot); + + for(i=0; i<n; i++) + { + for(j=0; j<n; j++) + { + double ainv_ij = 0.0; + + for(k = 0; k<n; k++) + { + double f = 1.0 / gsl_vector_get(eval, k); + double vik = gsl_matrix_get (evec, i, k); + double vjk = gsl_matrix_get (evec, j, k); + ainv_ij += vik * vjk * f; + } + gsl_matrix_set (ainv, i, j, ainv_ij); + } + } + + gsl_vector_free(eval); + gsl_matrix_free(evec); + gsl_matrix_free(tmp); + + if (status) + { + return status; + } + else + { + return GSL_SUCCESS; + } + } +} |