htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/sort.c

   1 /* eigen/sort.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman, Modified: B. Gough. */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_eigen.h>
  26 #include <gsl/gsl_complex.h>
  27 #include <gsl/gsl_complex_math.h>
  28
  29 /* The eigen_sort below is not very good, but it is simple and
  30  * self-contained. We can always implement an improved sort later.  */
  31
  32 int
  33 gsl_eigen_symmv_sort (gsl_vector * eval, gsl_matrix * evec,
  34                       gsl_eigen_sort_t sort_type)
  35 {
  36   if (evec->size1 != evec->size2)
  37     {
  38       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
  39     }
  40   else if (eval->size != evec->size1)
  41     {
  42       GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
  43     }
  44   else
  45     {
  46       const size_t N = eval->size;
  47       size_t i;
  48
  49       for (i = 0; i < N - 1; i++)
  50         {
  51           size_t j;
  52           size_t k = i;
  53
  54           double ek = gsl_vector_get (eval, i);
  55
  56           /* search for something to swap */
  57           for (j = i + 1; j < N; j++)
  58             {
  59               int test;
  60               const double ej = gsl_vector_get (eval, j);
  61
  62               switch (sort_type)
  63                 {
  64                 case GSL_EIGEN_SORT_VAL_ASC:
  65                   test = (ej < ek);
  66                   break;
  67                 case GSL_EIGEN_SORT_VAL_DESC:
  68                   test = (ej > ek);
  69                   break;
  70                 case GSL_EIGEN_SORT_ABS_ASC:
  71                   test = (fabs (ej) < fabs (ek));
  72                   break;
  73                 case GSL_EIGEN_SORT_ABS_DESC:
  74                   test = (fabs (ej) > fabs (ek));
  75                   break;
  76                 default:
  77                   GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
  78                 }
  79
  80               if (test)
  81                 {
  82                   k = j;
  83                   ek = ej;
  84                 }
  85             }
  86
  87           if (k != i)
  88             {
  89               /* swap eigenvalues */
  90               gsl_vector_swap_elements (eval, i, k);
  91
  92               /* swap eigenvectors */
  93               gsl_matrix_swap_columns (evec, i, k);
  94             }
  95         }
  96
  97       return GSL_SUCCESS;
  98     }
  99 }
 100
 101
 102 int
 103 gsl_eigen_hermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
 104                       gsl_eigen_sort_t sort_type)
 105 {
 106   if (evec->size1 != evec->size2)
 107     {
 108       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 109     }
 110   else if (eval->size != evec->size1)
 111     {
 112       GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
 113     }
 114   else
 115     {
 116       const size_t N = eval->size;
 117       size_t i;
 118
 119       for (i = 0; i < N - 1; i++)
 120         {
 121           size_t j;
 122           size_t k = i;
 123
 124           double ek = gsl_vector_get (eval, i);
 125
 126           /* search for something to swap */
 127           for (j = i + 1; j < N; j++)
 128             {
 129               int test;
 130               const double ej = gsl_vector_get (eval, j);
 131
 132               switch (sort_type)
 133                 {
 134                 case GSL_EIGEN_SORT_VAL_ASC:
 135                   test = (ej < ek);
 136                   break;
 137                 case GSL_EIGEN_SORT_VAL_DESC:
 138                   test = (ej > ek);
 139                   break;
 140                 case GSL_EIGEN_SORT_ABS_ASC:
 141                   test = (fabs (ej) < fabs (ek));
 142                   break;
 143                 case GSL_EIGEN_SORT_ABS_DESC:
 144                   test = (fabs (ej) > fabs (ek));
 145                   break;
 146                 default:
 147                   GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
 148                 }
 149
 150               if (test)
 151                 {
 152                   k = j;
 153                   ek = ej;
 154                 }
 155             }
 156
 157           if (k != i)
 158             {
 159               /* swap eigenvalues */
 160               gsl_vector_swap_elements (eval, i, k);
 161
 162               /* swap eigenvectors */
 163               gsl_matrix_complex_swap_columns (evec, i, k);
 164             }
 165         }
 166
 167       return GSL_SUCCESS;
 168     }
 169 }
 170
 171 int
 172 gsl_eigen_nonsymmv_sort (gsl_vector_complex * eval,
 173                          gsl_matrix_complex * evec,
 174                          gsl_eigen_sort_t sort_type)
 175 {
 176   if (evec && (evec->size1 != evec->size2))
 177     {
 178       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 179     }
 180   else if (evec && (eval->size != evec->size1))
 181     {
 182       GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
 183     }
 184   else
 185     {
 186       const size_t N = eval->size;
 187       size_t i;
 188
 189       for (i = 0; i < N - 1; i++)
 190         {
 191           size_t j;
 192           size_t k = i;
 193
 194           gsl_complex ek = gsl_vector_complex_get (eval, i);
 195
 196           /* search for something to swap */
 197           for (j = i + 1; j < N; j++)
 198             {
 199               int test;
 200               const gsl_complex ej = gsl_vector_complex_get (eval, j);
 201
 202               switch (sort_type)
 203                 {
 204                 case GSL_EIGEN_SORT_ABS_ASC:
 205                   test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
 206                   break;
 207                 case GSL_EIGEN_SORT_ABS_DESC:
 208                   test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
 209                   break;
 210                 case GSL_EIGEN_SORT_VAL_ASC:
 211                 case GSL_EIGEN_SORT_VAL_DESC:
 212                 default:
 213                   GSL_ERROR ("invalid sort type", GSL_EINVAL);
 214                 }
 215
 216               if (test)
 217                 {
 218                   k = j;
 219                   ek = ej;
 220                 }
 221             }
 222
 223           if (k != i)
 224             {
 225               /* swap eigenvalues */
 226               gsl_vector_complex_swap_elements (eval, i, k);
 227
 228               /* swap eigenvectors */
 229               if (evec)
 230                 gsl_matrix_complex_swap_columns (evec, i, k);
 231             }
 232         }
 233
 234       return GSL_SUCCESS;
 235     }
 236 }
 237
 238 int
 239 gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
 240                          gsl_eigen_sort_t sort_type)
 241 {
 242   int s;
 243
 244   s = gsl_eigen_symmv_sort(eval, evec, sort_type);
 245
 246   return s;
 247 }
 248
 249 int
 250 gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
 251                          gsl_eigen_sort_t sort_type)
 252 {
 253   int s;
 254
 255   s = gsl_eigen_hermv_sort(eval, evec, sort_type);
 256
 257   return s;
 258 }
 259
 260 int
 261 gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
 262                      gsl_matrix_complex * evec, gsl_eigen_sort_t sort_type)
 263 {
 264   if (evec->size1 != evec->size2)
 265     {
 266       GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
 267     }
 268   else if (alpha->size != evec->size1 || beta->size != evec->size1)
 269     {
 270       GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
 271     }
 272   else
 273     {
 274       const size_t N = alpha->size;
 275       size_t i;
 276
 277       for (i = 0; i < N - 1; i++)
 278         {
 279           size_t j;
 280           size_t k = i;
 281
 282           gsl_complex ak = gsl_vector_complex_get (alpha, i);
 283           double bk = gsl_vector_get(beta, i);
 284           gsl_complex ek;
 285
 286           if (bk < GSL_DBL_EPSILON)
 287             {
 288               GSL_SET_COMPLEX(&ek,
 289                               GSL_SIGN(GSL_REAL(ak)) ? GSL_POSINF : GSL_NEGINF,
 290                               GSL_SIGN(GSL_IMAG(ak)) ? GSL_POSINF : GSL_NEGINF);
 291             }
 292           else
 293             ek = gsl_complex_div_real(ak, bk);
 294
 295           /* search for something to swap */
 296           for (j = i + 1; j < N; j++)
 297             {
 298               int test;
 299               const gsl_complex aj = gsl_vector_complex_get (alpha, j);
 300               double bj = gsl_vector_get(beta, j);
 301               gsl_complex ej;
 302
 303               if (bj < GSL_DBL_EPSILON)
 304                 {
 305                   GSL_SET_COMPLEX(&ej,
 306                                   GSL_SIGN(GSL_REAL(aj)) ? GSL_POSINF : GSL_NEGINF,
 307                                   GSL_SIGN(GSL_IMAG(aj)) ? GSL_POSINF : GSL_NEGINF);
 308                 }
 309               else
 310                 ej = gsl_complex_div_real(aj, bj);
 311
 312               switch (sort_type)
 313                 {
 314                 case GSL_EIGEN_SORT_ABS_ASC:
 315                   test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
 316                   break;
 317                 case GSL_EIGEN_SORT_ABS_DESC:
 318                   test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
 319                   break;
 320                 case GSL_EIGEN_SORT_VAL_ASC:
 321                 case GSL_EIGEN_SORT_VAL_DESC:
 322                 default:
 323                   GSL_ERROR ("invalid sort type", GSL_EINVAL);
 324                 }
 325
 326               if (test)
 327                 {
 328                   k = j;
 329                   ek = ej;
 330                 }
 331             }
 332
 333           if (k != i)
 334             {
 335               /* swap eigenvalues */
 336               gsl_vector_complex_swap_elements (alpha, i, k);
 337               gsl_vector_swap_elements (beta, i, k);
 338
 339               /* swap eigenvectors */
 340               gsl_matrix_complex_swap_columns (evec, i, k);
 341             }
 342         }
 343
 344       return GSL_SUCCESS;
 345     }
 346 }