123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347 |
- /* eigen/sort.c
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Patrick Alken
- *
- * 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, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
- /* Author: G. Jungman, Modified: B. Gough. */
- #include "gsl__config.h"
- #include <stdlib.h>
- #include "gsl_math.h"
- #include "gsl_eigen.h"
- #include "gsl_complex.h"
- #include "gsl_complex_math.h"
- /* The eigen_sort below is not very good, but it is simple and
- * self-contained. We can always implement an improved sort later. */
- int
- gsl_eigen_symmv_sort (gsl_vector * eval, gsl_matrix * evec,
- gsl_eigen_sort_t sort_type)
- {
- if (evec->size1 != evec->size2)
- {
- GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
- }
- else if (eval->size != evec->size1)
- {
- GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
- }
- else
- {
- const size_t N = eval->size;
- size_t i;
- for (i = 0; i < N - 1; i++)
- {
- size_t j;
- size_t k = i;
- double ek = gsl_vector_get (eval, i);
- /* search for something to swap */
- for (j = i + 1; j < N; j++)
- {
- int test;
- const double ej = gsl_vector_get (eval, j);
- switch (sort_type)
- {
- case GSL_EIGEN_SORT_VAL_ASC:
- test = (ej < ek);
- break;
- case GSL_EIGEN_SORT_VAL_DESC:
- test = (ej > ek);
- break;
- case GSL_EIGEN_SORT_ABS_ASC:
- test = (fabs (ej) < fabs (ek));
- break;
- case GSL_EIGEN_SORT_ABS_DESC:
- test = (fabs (ej) > fabs (ek));
- break;
- default:
- GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
- }
- if (test)
- {
- k = j;
- ek = ej;
- }
- }
- if (k != i)
- {
- /* swap eigenvalues */
- gsl_vector_swap_elements (eval, i, k);
- /* swap eigenvectors */
- gsl_matrix_swap_columns (evec, i, k);
- }
- }
- return GSL_SUCCESS;
- }
- }
- int
- gsl_eigen_hermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
- gsl_eigen_sort_t sort_type)
- {
- if (evec->size1 != evec->size2)
- {
- GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
- }
- else if (eval->size != evec->size1)
- {
- GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
- }
- else
- {
- const size_t N = eval->size;
- size_t i;
- for (i = 0; i < N - 1; i++)
- {
- size_t j;
- size_t k = i;
- double ek = gsl_vector_get (eval, i);
- /* search for something to swap */
- for (j = i + 1; j < N; j++)
- {
- int test;
- const double ej = gsl_vector_get (eval, j);
- switch (sort_type)
- {
- case GSL_EIGEN_SORT_VAL_ASC:
- test = (ej < ek);
- break;
- case GSL_EIGEN_SORT_VAL_DESC:
- test = (ej > ek);
- break;
- case GSL_EIGEN_SORT_ABS_ASC:
- test = (fabs (ej) < fabs (ek));
- break;
- case GSL_EIGEN_SORT_ABS_DESC:
- test = (fabs (ej) > fabs (ek));
- break;
- default:
- GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
- }
- if (test)
- {
- k = j;
- ek = ej;
- }
- }
- if (k != i)
- {
- /* swap eigenvalues */
- gsl_vector_swap_elements (eval, i, k);
- /* swap eigenvectors */
- gsl_matrix_complex_swap_columns (evec, i, k);
- }
- }
- return GSL_SUCCESS;
- }
- }
- int
- gsl_eigen_nonsymmv_sort (gsl_vector_complex * eval,
- gsl_matrix_complex * evec,
- gsl_eigen_sort_t sort_type)
- {
- if (evec && (evec->size1 != evec->size2))
- {
- GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
- }
- else if (evec && (eval->size != evec->size1))
- {
- GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
- }
- else
- {
- const size_t N = eval->size;
- size_t i;
- for (i = 0; i < N - 1; i++)
- {
- size_t j;
- size_t k = i;
- gsl_complex ek = gsl_vector_complex_get (eval, i);
- /* search for something to swap */
- for (j = i + 1; j < N; j++)
- {
- int test;
- const gsl_complex ej = gsl_vector_complex_get (eval, j);
- switch (sort_type)
- {
- case GSL_EIGEN_SORT_ABS_ASC:
- test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
- break;
- case GSL_EIGEN_SORT_ABS_DESC:
- test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
- break;
- case GSL_EIGEN_SORT_VAL_ASC:
- case GSL_EIGEN_SORT_VAL_DESC:
- default:
- GSL_ERROR ("invalid sort type", GSL_EINVAL);
- }
- if (test)
- {
- k = j;
- ek = ej;
- }
- }
- if (k != i)
- {
- /* swap eigenvalues */
- gsl_vector_complex_swap_elements (eval, i, k);
- /* swap eigenvectors */
- if (evec)
- gsl_matrix_complex_swap_columns (evec, i, k);
- }
- }
- return GSL_SUCCESS;
- }
- }
- int
- gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
- gsl_eigen_sort_t sort_type)
- {
- int s;
- s = gsl_eigen_symmv_sort(eval, evec, sort_type);
- return s;
- }
- int
- gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
- gsl_eigen_sort_t sort_type)
- {
- int s;
- s = gsl_eigen_hermv_sort(eval, evec, sort_type);
- return s;
- }
- int
- gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
- gsl_matrix_complex * evec, gsl_eigen_sort_t sort_type)
- {
- if (evec->size1 != evec->size2)
- {
- GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
- }
- else if (alpha->size != evec->size1 || beta->size != evec->size1)
- {
- GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
- }
- else
- {
- const size_t N = alpha->size;
- size_t i;
- for (i = 0; i < N - 1; i++)
- {
- size_t j;
- size_t k = i;
- gsl_complex ak = gsl_vector_complex_get (alpha, i);
- double bk = gsl_vector_get(beta, i);
- gsl_complex ek;
- if (bk < GSL_DBL_EPSILON)
- {
- GSL_SET_COMPLEX(&ek,
- GSL_SIGN(GSL_REAL(ak)) ? GSL_POSINF : GSL_NEGINF,
- GSL_SIGN(GSL_IMAG(ak)) ? GSL_POSINF : GSL_NEGINF);
- }
- else
- ek = gsl_complex_div_real(ak, bk);
- /* search for something to swap */
- for (j = i + 1; j < N; j++)
- {
- int test;
- const gsl_complex aj = gsl_vector_complex_get (alpha, j);
- double bj = gsl_vector_get(beta, j);
- gsl_complex ej;
- if (bj < GSL_DBL_EPSILON)
- {
- GSL_SET_COMPLEX(&ej,
- GSL_SIGN(GSL_REAL(aj)) ? GSL_POSINF : GSL_NEGINF,
- GSL_SIGN(GSL_IMAG(aj)) ? GSL_POSINF : GSL_NEGINF);
- }
- else
- ej = gsl_complex_div_real(aj, bj);
- switch (sort_type)
- {
- case GSL_EIGEN_SORT_ABS_ASC:
- test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
- break;
- case GSL_EIGEN_SORT_ABS_DESC:
- test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
- break;
- case GSL_EIGEN_SORT_VAL_ASC:
- case GSL_EIGEN_SORT_VAL_DESC:
- default:
- GSL_ERROR ("invalid sort type", GSL_EINVAL);
- }
- if (test)
- {
- k = j;
- ek = ej;
- }
- }
- if (k != i)
- {
- /* swap eigenvalues */
- gsl_vector_complex_swap_elements (alpha, i, k);
- gsl_vector_swap_elements (beta, i, k);
- /* swap eigenvectors */
- gsl_matrix_complex_swap_columns (evec, i, k);
- }
- }
- return GSL_SUCCESS;
- }
- }
|