/* 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 <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
/* Compares real parts of a and b and, if they are not approximately equal,
* returns Re(a) < Re(b); otherwise returns Im(a) < Im(b). */
static int
complex_less(gsl_complex a, gsl_complex b)
{
return gsl_fcmp(GSL_REAL(a), GSL_REAL(b), GSL_DBL_EPSILON) == 0 ?
GSL_IMAG(a) < GSL_IMAG(b) : GSL_REAL(a) < GSL_REAL(b);
}
/* 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:
test = complex_less(ej, ek);
break;
case GSL_EIGEN_SORT_VAL_DESC:
test = complex_less(ek, ej);
break;
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;
}
}