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

   1 /* eigen/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Patrick Alken, Brian Gough
   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
  21  */
  22
  23 #include <config.h>
  24 #include <stdlib.h>
  25 #include <gsl/gsl_test.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_blas.h>
  28 #include <gsl/gsl_ieee_utils.h>
  29 #include <gsl/gsl_complex_math.h>
  30 #include <gsl/gsl_eigen.h>
  31 #include <gsl/gsl_linalg.h>
  32 #include <gsl/gsl_rng.h>
  33 #include <gsl/gsl_randist.h>
  34 #include <gsl/gsl_errno.h>
  35 #include <gsl/gsl_sort.h>
  36 #include <gsl/gsl_sort_vector.h>
  37 #include <gsl/gsl_matrix.h>
  38 #include <gsl/gsl_vector.h>
  39
  40 /******************************************
  41  * common test code                       *
  42  ******************************************/
  43
  44 double
  45 chop_subnormals (double x)
  46 {
  47   /* Chop any subnormal values */
  48   return fabs(x) < GSL_DBL_MIN ? 0 : x;
  49 }
  50
  51 void
  52 create_random_symm_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
  53 {
  54   size_t i, j;
  55
  56   for (i = 0; i < m->size1; ++i)
  57     {
  58       for (j = i; j < m->size2; ++j)
  59       {
  60         double x = gsl_rng_uniform(r) * (upper - lower) + lower;
  61         gsl_matrix_set(m, i, j, x);
  62         gsl_matrix_set(m, j, i, x);
  63       }
  64     }
  65 } /* create_random_symm_matrix() */
  66
  67 void
  68 create_random_herm_matrix(gsl_matrix_complex *m, gsl_rng *r, int lower,
  69                           int upper)
  70 {
  71   size_t i, j;
  72
  73   for (i = 0; i < m->size1; ++i)
  74     {
  75       for (j = i; j < m->size2; ++j)
  76       {
  77         gsl_complex z;
  78
  79         GSL_REAL(z) = gsl_rng_uniform(r) * (upper - lower) + lower;
  80
  81         if (i == j)
  82           GSL_IMAG(z) = 0.0;
  83         else
  84           GSL_IMAG(z) = gsl_rng_uniform(r) * (upper - lower) + lower;
  85
  86         gsl_matrix_complex_set(m, i, j, z);
  87         gsl_matrix_complex_set(m, j, i, gsl_complex_conjugate(z));
  88       }
  89     }
  90 } /* create_random_herm_matrix() */
  91
  92 /* with r \in (0,1) if m_{ij} = r^{|i - j|} then m is positive definite */
  93 void
  94 create_random_posdef_matrix(gsl_matrix *m, gsl_rng *r)
  95 {
  96   size_t i, j;
  97   double x = gsl_rng_uniform(r);
  98
  99   for (i = 0; i < m->size1; ++i)
 100     {
 101       for (j = i; j < m->size2; ++j)
 102       {
 103         double a = pow(x, (double) (j - i));
 104
 105         gsl_matrix_set(m, i, j, a);
 106         gsl_matrix_set(m, j, i, a);
 107       }
 108     }
 109 } /* create_random_posdef_matrix() */
 110
 111 void
 112 create_random_complex_posdef_matrix(gsl_matrix_complex *m, gsl_rng *r,
 113                                     gsl_vector_complex *work)
 114 {
 115   const size_t N = m->size1;
 116   size_t i, j;
 117   double x, y;
 118   gsl_complex z;
 119   gsl_complex tau;
 120
 121   GSL_SET_IMAG(&z, 0.0);
 122
 123   /* make a positive diagonal matrix */
 124   gsl_matrix_complex_set_zero(m);
 125   for (i = 0; i < N; ++i)
 126     {
 127       x = gsl_rng_uniform(r);
 128       GSL_SET_REAL(&z, x);
 129       gsl_matrix_complex_set(m, i, i, z);
 130     }
 131
 132   /* now generate random householder reflections and form P D P^H */
 133   for (i = 0; i < N; ++i)
 134     {
 135       /* form complex vector */
 136       for (j = 0; j < N; ++j)
 137         {
 138           x = 2.0 * gsl_rng_uniform(r) - 1.0;
 139           y = 2.0 * gsl_rng_uniform(r) - 1.0;
 140           GSL_SET_COMPLEX(&z, x, y);
 141           gsl_vector_complex_set(work, j, z);
 142         }
 143
 144       tau = gsl_linalg_complex_householder_transform(work);
 145       gsl_linalg_complex_householder_hm(tau, work, m);
 146       gsl_linalg_complex_householder_mh(gsl_complex_conjugate(tau), work, m);
 147     }
 148 } /* create_random_complex_posdef_matrix() */
 149
 150 void
 151 create_random_nonsymm_matrix(gsl_matrix *m, gsl_rng *r, int lower,
 152                              int upper)
 153 {
 154   size_t i, j;
 155
 156   for (i = 0; i < m->size1; ++i)
 157     {
 158       for (j = 0; j < m->size2; ++j)
 159       {
 160         gsl_matrix_set(m,
 161                        i,
 162                        j,
 163                        gsl_rng_uniform(r) * (upper - lower) + lower);
 164       }
 165     }
 166 } /* create_random_nonsymm_matrix() */
 167
 168 /* test if A Z = Q S */
 169 void
 170 test_eigen_schur(const gsl_matrix * A, const gsl_matrix * S,
 171                  const gsl_matrix * Q, const gsl_matrix * Z,
 172                  size_t count, const char * desc,
 173                  const char * desc2)
 174 {
 175   const size_t N = A->size1;
 176   size_t i, j;
 177
 178   gsl_matrix * T1 = gsl_matrix_alloc(N, N);
 179   gsl_matrix * T2 = gsl_matrix_alloc(N, N);
 180
 181   /* compute T1 = A Z */
 182   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, Z, 0.0, T1);
 183
 184   /* compute T2 = Q S */
 185   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Q, S, 0.0, T2);
 186
 187   for (i = 0; i < N; ++i)
 188     {
 189       for (j = 0; j < N; ++j)
 190         {
 191           double x = gsl_matrix_get(T1, i, j);
 192           double y = gsl_matrix_get(T2, i, j);
 193
 194           gsl_test_abs(x, y, 1.0e8 * GSL_DBL_EPSILON,
 195                        "%s(N=%u,cnt=%u), %s, schur(%d,%d)", desc, N, count, desc2, i, j);
 196         }
 197     }
 198
 199   gsl_matrix_free (T1);
 200   gsl_matrix_free (T2);
 201 } /* test_eigen_schur() */
 202
 203 void
 204 test_eigenvalues_real (const gsl_vector *eval, const gsl_vector * eval2,
 205                        const char * desc, const char * desc2)
 206 {
 207   const size_t N = eval->size;
 208   size_t i;
 209
 210   double emax = 0;
 211
 212   /* check eigenvalues */
 213   for (i = 0; i < N; i++)
 214     {
 215       double ei = gsl_vector_get (eval, i);
 216       if (fabs(ei) > emax) emax = fabs(ei);
 217     }
 218
 219   for (i = 0; i < N; i++)
 220     {
 221       double ei = gsl_vector_get (eval, i);
 222       double e2i = gsl_vector_get (eval2, i);
 223       e2i = chop_subnormals(e2i);
 224       gsl_test_abs(ei, e2i, emax * 1e8 * GSL_DBL_EPSILON,
 225                    "%s, direct eigenvalue(%d), %s",
 226                    desc, i, desc2);
 227     }
 228 }
 229
 230 void
 231 test_eigenvalues_complex (const gsl_vector_complex * eval,
 232                           const gsl_vector_complex * eval2,
 233                           const char * desc, const char * desc2)
 234 {
 235   const size_t N = eval->size;
 236   size_t i;
 237
 238   for (i = 0; i < N; i++)
 239     {
 240       gsl_complex ei = gsl_vector_complex_get (eval, i);
 241       gsl_complex e2i = gsl_vector_complex_get (eval2, i);
 242       gsl_test_rel(GSL_REAL(ei), GSL_REAL(e2i), 10*N*GSL_DBL_EPSILON,
 243                    "%s, direct eigenvalue(%d) real, %s",
 244                    desc, i, desc2);
 245       gsl_test_rel(GSL_IMAG(ei), GSL_IMAG(e2i), 10*N*GSL_DBL_EPSILON,
 246                    "%s, direct eigenvalue(%d) imag, %s",
 247                    desc, i, desc2);
 248     }
 249 }
 250
 251 /******************************************
 252  * symm test code                         *
 253  ******************************************/
 254
 255 void
 256 test_eigen_symm_results (const gsl_matrix * A,
 257                          const gsl_vector * eval,
 258                          const gsl_matrix * evec,
 259                          size_t count,
 260                          const char * desc,
 261                          const char * desc2)
 262 {
 263   const size_t N = A->size1;
 264   size_t i, j;
 265   double emax = 0;
 266
 267   gsl_vector * x = gsl_vector_alloc(N);
 268   gsl_vector * y = gsl_vector_alloc(N);
 269
 270   /* check eigenvalues */
 271   for (i = 0; i < N; i++)
 272     {
 273       double ei = gsl_vector_get (eval, i);
 274       if (fabs(ei) > emax) emax = fabs(ei);
 275     }
 276
 277   for (i = 0; i < N; i++)
 278     {
 279       double ei = gsl_vector_get (eval, i);
 280       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
 281       gsl_vector_memcpy(x, &vi.vector);
 282       /* compute y = A x (should = lambda v) */
 283       gsl_blas_dgemv (CblasNoTrans, 1.0, A, x, 0.0, y);
 284       for (j = 0; j < N; j++)
 285         {
 286           double xj = gsl_vector_get (x, j);
 287           double yj = gsl_vector_get (y, j);
 288           double eixj = chop_subnormals(ei * xj);
 289           gsl_test_abs(yj, eixj,  emax * 1e8 * GSL_DBL_EPSILON,
 290                        "%s, eigenvalue(%d,%d), %s", desc, i, j, desc2);
 291         }
 292     }
 293
 294   /* check eigenvectors are orthonormal */
 295
 296   for (i = 0; i < N; i++)
 297     {
 298       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
 299       double nrm_v = gsl_blas_dnrm2(&vi.vector);
 300       gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
 301                     desc, i, desc2);
 302     }
 303
 304   for (i = 0; i < N; i++)
 305     {
 306       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
 307       for (j = i + 1; j < N; j++)
 308         {
 309           gsl_vector_const_view vj = gsl_matrix_const_column(evec, j);
 310           double vivj;
 311           gsl_blas_ddot (&vi.vector, &vj.vector, &vivj);
 312           gsl_test_abs (vivj, 0.0, N * GSL_DBL_EPSILON,
 313                         "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
 314         }
 315     }
 316
 317   gsl_vector_free(x);
 318   gsl_vector_free(y);
 319 }
 320
 321 void
 322 test_eigen_symm_matrix(const gsl_matrix * m, size_t count,
 323                        const char * desc)
 324 {
 325   const size_t N = m->size1;
 326   gsl_matrix * A = gsl_matrix_alloc(N, N);
 327   gsl_vector * eval = gsl_vector_alloc(N);
 328   gsl_vector * evalv = gsl_vector_alloc(N);
 329   gsl_vector * x = gsl_vector_alloc(N);
 330   gsl_vector * y = gsl_vector_alloc(N);
 331   gsl_matrix * evec = gsl_matrix_alloc(N, N);
 332   gsl_eigen_symm_workspace * w = gsl_eigen_symm_alloc(N);
 333   gsl_eigen_symmv_workspace * wv = gsl_eigen_symmv_alloc(N);
 334
 335   gsl_matrix_memcpy(A, m);
 336
 337   gsl_eigen_symmv(A, evalv, evec, wv);
 338   test_eigen_symm_results(m, evalv, evec, count, desc, "unsorted");
 339
 340   gsl_matrix_memcpy(A, m);
 341
 342   gsl_eigen_symm(A, eval, w);
 343
 344   /* sort eval and evalv */
 345   gsl_vector_memcpy(x, eval);
 346   gsl_vector_memcpy(y, evalv);
 347   gsl_sort_vector(x);
 348   gsl_sort_vector(y);
 349   test_eigenvalues_real(y, x, desc, "unsorted");
 350
 351   gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
 352   test_eigen_symm_results(m, evalv, evec, count, desc, "val/asc");
 353
 354   gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
 355   test_eigen_symm_results(m, evalv, evec, count, desc, "val/desc");
 356
 357   gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
 358   test_eigen_symm_results(m, evalv, evec, count, desc, "abs/asc");
 359
 360   gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
 361   test_eigen_symm_results(m, evalv, evec, count, desc, "abs/desc");
 362
 363   gsl_matrix_free(A);
 364   gsl_vector_free(eval);
 365   gsl_vector_free(evalv);
 366   gsl_vector_free(x);
 367   gsl_vector_free(y);
 368   gsl_matrix_free(evec);
 369   gsl_eigen_symm_free(w);
 370   gsl_eigen_symmv_free(wv);
 371 } /* test_eigen_symm_matrix() */
 372
 373 void
 374 test_eigen_symm(void)
 375 {
 376   size_t N_max = 20;
 377   size_t n, i;
 378   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
 379
 380   for (n = 1; n <= N_max; ++n)
 381     {
 382       gsl_matrix * A = gsl_matrix_alloc(n, n);
 383
 384       for (i = 0; i < 5; ++i)
 385         {
 386           create_random_symm_matrix(A, r, -10, 10);
 387           test_eigen_symm_matrix(A, i, "symm random");
 388         }
 389
 390       gsl_matrix_free(A);
 391     }
 392
 393   gsl_rng_free(r);
 394
 395   {
 396     double dat1[] =  { 0,  0, -1,  0,
 397                        0,  1,  0,  1,
 398                        -1,  0,  0,  0,
 399                        0,  1,  0,  0 };
 400     double dat2[] =  { 1,  0,  0,  0,
 401                        0,  2,  0,  0,
 402                        0,  0,  3,  0,
 403                        0,  0,  0,  4 };
 404     gsl_matrix_view m;
 405
 406     m = gsl_matrix_view_array (dat1, 4, 4);
 407     test_eigen_symm_matrix(&m.matrix, 0, "symm(4)");
 408
 409     m = gsl_matrix_view_array (dat2, 4, 4);
 410     test_eigen_symm_matrix(&m.matrix, 0, "symm(4) diag");
 411   }
 412
 413   {
 414     double dat[27*27] = {
 415         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 416         0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,
 417         0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,
 418         0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 419         0,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
 420         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 421         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 422         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 423         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 424         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 425         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 426         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 427         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 428         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 429         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 430         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 431         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 432         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 433         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 434         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 435         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 436         0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 437         0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 438         0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 439         0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 440         0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 441         0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
 442     };
 443
 444     gsl_matrix_view m;
 445     m = gsl_matrix_view_array (dat, 27, 27);
 446     test_eigen_symm_matrix(&m.matrix, 0, "symm(27)");
 447   };
 448
 449 } /* test_eigen_symm() */
 450
 451 /******************************************
 452  * herm test code                         *
 453  ******************************************/
 454
 455 void
 456 test_eigen_herm_results (const gsl_matrix_complex * A,
 457                          const gsl_vector * eval,
 458                          const gsl_matrix_complex * evec,
 459                          size_t count,
 460                          const char * desc,
 461                          const char * desc2)
 462 {
 463   const size_t N = A->size1;
 464   size_t i, j;
 465
 466   gsl_vector_complex * x = gsl_vector_complex_alloc(N);
 467   gsl_vector_complex * y = gsl_vector_complex_alloc(N);
 468
 469   /* check eigenvalues */
 470
 471   for (i = 0; i < N; i++)
 472     {
 473       double ei = gsl_vector_get (eval, i);
 474       gsl_vector_complex_const_view vi =
 475         gsl_matrix_complex_const_column(evec, i);
 476       gsl_vector_complex_memcpy(x, &vi.vector);
 477       /* compute y = m x (should = lambda v) */
 478       gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x,
 479                       GSL_COMPLEX_ZERO, y);
 480       for (j = 0; j < N; j++)
 481         {
 482           gsl_complex xj = gsl_vector_complex_get (x, j);
 483           gsl_complex yj = gsl_vector_complex_get (y, j);
 484           gsl_test_rel(GSL_REAL(yj), ei * GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
 485                        "%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2);
 486           gsl_test_rel(GSL_IMAG(yj), ei * GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
 487                        "%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2);
 488         }
 489     }
 490
 491   /* check eigenvectors are orthonormal */
 492
 493   for (i = 0; i < N; i++)
 494     {
 495       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
 496       double nrm_v = gsl_blas_dznrm2(&vi.vector);
 497       gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
 498                     desc, i, desc2);
 499     }
 500
 501   for (i = 0; i < N; i++)
 502     {
 503       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
 504       for (j = i + 1; j < N; j++)
 505         {
 506           gsl_vector_complex_const_view vj
 507             = gsl_matrix_complex_const_column(evec, j);
 508           gsl_complex vivj;
 509           gsl_blas_zdotc (&vi.vector, &vj.vector, &vivj);
 510           gsl_test_abs (gsl_complex_abs(vivj), 0.0, N * GSL_DBL_EPSILON,
 511                         "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
 512         }
 513     }
 514
 515   gsl_vector_complex_free(x);
 516   gsl_vector_complex_free(y);
 517 } /* test_eigen_herm_results() */
 518
 519 void
 520 test_eigen_herm_matrix(const gsl_matrix_complex * m, size_t count,
 521                        const char * desc)
 522 {
 523   const size_t N = m->size1;
 524   gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
 525   gsl_vector * eval = gsl_vector_alloc(N);
 526   gsl_vector * evalv = gsl_vector_alloc(N);
 527   gsl_vector * x = gsl_vector_alloc(N);
 528   gsl_vector * y = gsl_vector_alloc(N);
 529   gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
 530   gsl_eigen_herm_workspace * w = gsl_eigen_herm_alloc(N);
 531   gsl_eigen_hermv_workspace * wv = gsl_eigen_hermv_alloc(N);
 532
 533   gsl_matrix_complex_memcpy(A, m);
 534
 535   gsl_eigen_hermv(A, evalv, evec, wv);
 536   test_eigen_herm_results(m, evalv, evec, count, desc, "unsorted");
 537
 538   gsl_matrix_complex_memcpy(A, m);
 539
 540   gsl_eigen_herm(A, eval, w);
 541
 542   /* sort eval and evalv */
 543   gsl_vector_memcpy(x, eval);
 544   gsl_vector_memcpy(y, evalv);
 545   gsl_sort_vector(x);
 546   gsl_sort_vector(y);
 547   test_eigenvalues_real(y, x, desc, "unsorted");
 548
 549   gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
 550   test_eigen_herm_results(m, evalv, evec, count, desc, "val/asc");
 551
 552   gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
 553   test_eigen_herm_results(m, evalv, evec, count, desc, "val/desc");
 554
 555   gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
 556   test_eigen_herm_results(m, evalv, evec, count, desc, "abs/asc");
 557
 558   gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
 559   test_eigen_herm_results(m, evalv, evec, count, desc, "abs/desc");
 560
 561   gsl_matrix_complex_free(A);
 562   gsl_vector_free(eval);
 563   gsl_vector_free(evalv);
 564   gsl_vector_free(x);
 565   gsl_vector_free(y);
 566   gsl_matrix_complex_free(evec);
 567   gsl_eigen_herm_free(w);
 568   gsl_eigen_hermv_free(wv);
 569 } /* test_eigen_herm_matrix() */
 570
 571 void
 572 test_eigen_herm(void)
 573 {
 574   size_t N_max = 20;
 575   size_t n, i;
 576   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
 577
 578   for (n = 1; n <= N_max; ++n)
 579     {
 580       gsl_matrix_complex * A = gsl_matrix_complex_alloc(n, n);
 581
 582       for (i = 0; i < 5; ++i)
 583         {
 584           create_random_herm_matrix(A, r, -10, 10);
 585           test_eigen_herm_matrix(A, i, "herm random");
 586         }
 587
 588       gsl_matrix_complex_free(A);
 589     }
 590
 591   gsl_rng_free(r);
 592
 593   {
 594     double dat1[] =  { 0,0,  0,0, -1,0,  0,0,
 595                        0,0,  1,0,  0,0,  1,0,
 596                        -1,0,  0,0,  0,0,  0,0,
 597                        0,0,  1,0,  0,0,  0,0 };
 598     double dat2[] =  { 1,0,  0,0, 0,0,  0,0,
 599                        0,0,  2,0, 0,0,  0,0,
 600                        0,0,  0,0, 3,0,  0,0,
 601                        0,0,  0,0, 0,0,  4,0 };
 602     gsl_matrix_complex_view m;
 603
 604     m = gsl_matrix_complex_view_array (dat1, 4, 4);
 605     test_eigen_herm_matrix(&m.matrix, 0, "herm(4)");
 606
 607     m = gsl_matrix_complex_view_array (dat2, 4, 4);
 608     test_eigen_herm_matrix(&m.matrix, 1, "herm(4) diag");
 609   }
 610 } /* test_eigen_herm() */
 611
 612 /******************************************
 613  * nonsymm test code                      *
 614  ******************************************/
 615
 616 void
 617 test_eigen_nonsymm_results (const gsl_matrix * m,
 618                             const gsl_vector_complex * eval,
 619                             const gsl_matrix_complex * evec,
 620                             size_t count,
 621                             const char * desc,
 622                             const char * desc2)
 623 {
 624   size_t i,j;
 625   size_t N = m->size1;
 626
 627   gsl_vector_complex * x = gsl_vector_complex_alloc(N);
 628   gsl_vector_complex * y = gsl_vector_complex_alloc(N);
 629   gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
 630
 631   /* we need a complex matrix for the blas routines, so copy m into A */
 632   for (i = 0; i < N; ++i)
 633     {
 634       for (j = 0; j < N; ++j)
 635         {
 636           gsl_complex z;
 637           GSL_SET_COMPLEX(&z, gsl_matrix_get(m, i, j), 0.0);
 638           gsl_matrix_complex_set(A, i, j, z);
 639         }
 640     }
 641
 642   for (i = 0; i < N; i++)
 643     {
 644       gsl_complex ei = gsl_vector_complex_get (eval, i);
 645       gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
 646       double norm = gsl_blas_dznrm2(&vi.vector);
 647
 648       /* check that eigenvector is normalized */
 649       gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
 650                    "nonsymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count, desc, i, desc2);
 651
 652       gsl_vector_complex_memcpy(x, &vi.vector);
 653
 654       /* compute y = m x (should = lambda v) */
 655       gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x,
 656                       GSL_COMPLEX_ZERO, y);
 657
 658       /* compute x = lambda v */
 659       gsl_blas_zscal(ei, x);
 660
 661       /* now test if y = x */
 662       for (j = 0; j < N; j++)
 663         {
 664           gsl_complex xj = gsl_vector_complex_get (x, j);
 665           gsl_complex yj = gsl_vector_complex_get (y, j);
 666
 667           /* use abs here in case the values are close to 0 */
 668           gsl_test_abs(GSL_REAL(yj), GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
 669                        "nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
 670           gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
 671                        "nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
 672         }
 673     }
 674
 675   gsl_matrix_complex_free(A);
 676   gsl_vector_complex_free(x);
 677   gsl_vector_complex_free(y);
 678 }
 679
 680 void
 681 test_eigen_nonsymm_matrix(const gsl_matrix * m, size_t count,
 682                           const char * desc,
 683                           gsl_eigen_nonsymmv_workspace *w)
 684 {
 685   const size_t N = m->size1;
 686   gsl_matrix * A = gsl_matrix_alloc(N, N);
 687   gsl_matrix * Z = gsl_matrix_alloc(N, N);
 688   gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
 689   gsl_vector_complex * eval = gsl_vector_complex_alloc(N);
 690
 691   /*
 692    * calculate eigenvalues and eigenvectors - it is sufficient to
 693    * test gsl_eigen_nonsymmv() since that function calls
 694    * gsl_eigen_nonsymm() for the eigenvalues
 695    */
 696   gsl_matrix_memcpy(A, m);
 697   gsl_eigen_nonsymmv(A, eval, evec, w);
 698   test_eigen_nonsymm_results (m, eval, evec, count, desc, "unsorted");
 699
 700   /* test sort routines */
 701   gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
 702   test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/asc");
 703
 704   gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);
 705   test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/desc");
 706
 707   /* test Schur vectors */
 708   gsl_matrix_memcpy(A, m);
 709   gsl_eigen_nonsymmv_Z(A, eval, evec, Z, w);
 710   gsl_linalg_hessenberg_set_zero(A);
 711   test_eigen_schur(m, A, Z, Z, count, "nonsymm", desc);
 712
 713   gsl_matrix_free(A);
 714   gsl_matrix_free(Z);
 715   gsl_matrix_complex_free(evec);
 716   gsl_vector_complex_free(eval);
 717 }
 718
 719 void
 720 test_eigen_nonsymm(void)
 721 {
 722   size_t N_max = 20;
 723   size_t n, i;
 724   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
 725
 726   for (n = 1; n <= N_max; ++n)
 727     {
 728       gsl_matrix * m = gsl_matrix_alloc(n, n);
 729       gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(n);
 730
 731       for (i = 0; i < 5; ++i)
 732         {
 733           create_random_nonsymm_matrix(m, r, -10, 10);
 734
 735           gsl_eigen_nonsymm_params(1, 0, w->nonsymm_workspace_p);
 736           test_eigen_nonsymm_matrix(m, i, "random, unbalanced", w);
 737
 738           gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);
 739           test_eigen_nonsymm_matrix(m, i, "random, balanced", w);
 740         }
 741
 742       gsl_matrix_free(m);
 743       gsl_eigen_nonsymmv_free(w);
 744     }
 745
 746   gsl_rng_free(r);
 747
 748   {
 749     double dat1[] = { 0, 1, 1, 1,
 750                       1, 1, 1, 1,
 751                       0, 0, 0, 0,
 752                       0, 0, 0, 0 };
 753     double dat2[] = { 1, 1, 0, 1,
 754                       1, 1, 1, 1,
 755                       1, 1, 1, 1,
 756                       0, 1, 0, 0 };
 757     gsl_matrix_view v;
 758     gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(4);
 759
 760     v = gsl_matrix_view_array (dat1, 4, 4);
 761     test_eigen_nonsymm_matrix(&v.matrix, 0, "integer", w);
 762
 763     v = gsl_matrix_view_array (dat2, 4, 4);
 764     test_eigen_nonsymm_matrix(&v.matrix, 1, "integer", w);
 765
 766     gsl_eigen_nonsymmv_free(w);
 767   }
 768 } /* test_eigen_nonsymm() */
 769
 770 /******************************************
 771  * gensymm test code                      *
 772  ******************************************/
 773
 774 void
 775 test_eigen_gensymm_results (const gsl_matrix * A,
 776                             const gsl_matrix * B,
 777                             const gsl_vector * eval,
 778                             const gsl_matrix * evec,
 779                             size_t count,
 780                             const char * desc,
 781                             const char * desc2)
 782 {
 783   const size_t N = A->size1;
 784   size_t i, j;
 785
 786   gsl_vector * x = gsl_vector_alloc(N);
 787   gsl_vector * y = gsl_vector_alloc(N);
 788   gsl_vector * z = gsl_vector_alloc(N);
 789
 790   /* check A v = lambda B v */
 791   for (i = 0; i < N; i++)
 792     {
 793       double ei = gsl_vector_get (eval, i);
 794       gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
 795       double norm = gsl_blas_dnrm2(&vi.vector);
 796
 797       /* check that eigenvector is normalized */
 798       gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
 799                    "gensymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
 800                    desc, i, desc2);
 801
 802       gsl_vector_memcpy(z, &vi.vector);
 803
 804       /* compute y = A z */
 805       gsl_blas_dgemv (CblasNoTrans, 1.0, A, z, 0.0, y);
 806
 807       /* compute x = B z */
 808       gsl_blas_dgemv (CblasNoTrans, 1.0, B, z, 0.0, x);
 809
 810       /* compute x = lambda B z */
 811       gsl_blas_dscal(ei, x);
 812
 813       /* now test if y = x */
 814       for (j = 0; j < N; j++)
 815         {
 816           double xj = gsl_vector_get (x, j);
 817           double yj = gsl_vector_get (y, j);
 818
 819           gsl_test_rel(yj, xj, 1e9 * GSL_DBL_EPSILON,
 820                        "gensymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
 821         }
 822     }
 823
 824   gsl_vector_free(x);
 825   gsl_vector_free(y);
 826   gsl_vector_free(z);
 827 }
 828
 829 void
 830 test_eigen_gensymm(void)
 831 {
 832   size_t N_max = 20;
 833   size_t n, i;
 834   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
 835
 836   for (n = 1; n <= N_max; ++n)
 837     {
 838       gsl_matrix * A = gsl_matrix_alloc(n, n);
 839       gsl_matrix * B = gsl_matrix_alloc(n, n);
 840       gsl_matrix * ma = gsl_matrix_alloc(n, n);
 841       gsl_matrix * mb = gsl_matrix_alloc(n, n);
 842       gsl_vector * eval = gsl_vector_alloc(n);
 843       gsl_vector * evalv = gsl_vector_alloc(n);
 844       gsl_vector * x = gsl_vector_alloc(n);
 845       gsl_vector * y = gsl_vector_alloc(n);
 846       gsl_matrix * evec = gsl_matrix_alloc(n, n);
 847       gsl_eigen_gensymm_workspace * w = gsl_eigen_gensymm_alloc(n);
 848       gsl_eigen_gensymmv_workspace * wv = gsl_eigen_gensymmv_alloc(n);
 849
 850       for (i = 0; i < 5; ++i)
 851         {
 852           create_random_symm_matrix(A, r, -10, 10);
 853           create_random_posdef_matrix(B, r);
 854
 855           gsl_matrix_memcpy(ma, A);
 856           gsl_matrix_memcpy(mb, B);
 857
 858           gsl_eigen_gensymmv(ma, mb, evalv, evec, wv);
 859           test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "unsorted");
 860
 861           gsl_matrix_memcpy(ma, A);
 862           gsl_matrix_memcpy(mb, B);
 863
 864           gsl_eigen_gensymm(ma, mb, eval, w);
 865
 866           /* eval and evalv have to be sorted? not sure why */
 867           gsl_vector_memcpy(x, eval);
 868           gsl_vector_memcpy(y, evalv);
 869           gsl_sort_vector(x);
 870           gsl_sort_vector(y);
 871           test_eigenvalues_real(y, x, "gensymm, random", "unsorted");
 872
 873           gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
 874           test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "val/asc");
 875
 876           gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
 877           test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "val/desc");
 878
 879           gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
 880           test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "abs/asc");
 881           gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
 882           test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "abs/desc");
 883         }
 884
 885       gsl_matrix_free(A);
 886       gsl_matrix_free(B);
 887       gsl_matrix_free(ma);
 888       gsl_matrix_free(mb);
 889       gsl_vector_free(eval);
 890       gsl_vector_free(evalv);
 891       gsl_vector_free(x);
 892       gsl_vector_free(y);
 893       gsl_matrix_free(evec);
 894       gsl_eigen_gensymm_free(w);
 895       gsl_eigen_gensymmv_free(wv);
 896     }
 897
 898   gsl_rng_free(r);
 899 } /* test_eigen_gensymm() */
 900
 901 /******************************************
 902  * genherm test code                      *
 903  ******************************************/
 904
 905 void
 906 test_eigen_genherm_results (const gsl_matrix_complex * A,
 907                             const gsl_matrix_complex * B,
 908                             const gsl_vector * eval,
 909                             const gsl_matrix_complex * evec,
 910                             size_t count,
 911                             const char * desc,
 912                             const char * desc2)
 913 {
 914   const size_t N = A->size1;
 915   size_t i, j;
 916
 917   gsl_vector_complex * x = gsl_vector_complex_alloc(N);
 918   gsl_vector_complex * y = gsl_vector_complex_alloc(N);
 919
 920   /* check A v = lambda B v */
 921   for (i = 0; i < N; i++)
 922     {
 923       double ei = gsl_vector_get (eval, i);
 924       gsl_vector_complex_const_view vi =
 925         gsl_matrix_complex_const_column(evec, i);
 926       double norm = gsl_blas_dznrm2(&vi.vector);
 927
 928       /* check that eigenvector is normalized */
 929       gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
 930                    "genherm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
 931                    desc, i, desc2);
 932
 933       /* compute y = A z */
 934       gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, &vi.vector, GSL_COMPLEX_ZERO, y);
 935
 936       /* compute x = B z */
 937       gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, B, &vi.vector, GSL_COMPLEX_ZERO, x);
 938
 939       /* compute x = lambda B z */
 940       gsl_blas_zdscal(ei, x);
 941
 942       /* now test if y = x */
 943       for (j = 0; j < N; j++)
 944         {
 945           gsl_complex xj = gsl_vector_complex_get (x, j);
 946           gsl_complex yj = gsl_vector_complex_get (y, j);
 947
 948           gsl_test_rel(GSL_REAL(yj), GSL_REAL(xj), 1e9 * GSL_DBL_EPSILON,
 949                        "genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
 950           gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e9 * GSL_DBL_EPSILON,
 951                        "genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
 952         }
 953     }
 954
 955   gsl_vector_complex_free(x);
 956   gsl_vector_complex_free(y);
 957 }
 958
 959 void
 960 test_eigen_genherm(void)
 961 {
 962   size_t N_max = 20;
 963   size_t n, i;
 964   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
 965
 966   for (n = 1; n <= N_max; ++n)
 967     {
 968       gsl_matrix_complex * A = gsl_matrix_complex_alloc(n, n);
 969       gsl_matrix_complex * B = gsl_matrix_complex_alloc(n, n);
 970       gsl_matrix_complex * ma = gsl_matrix_complex_alloc(n, n);
 971       gsl_matrix_complex * mb = gsl_matrix_complex_alloc(n, n);
 972       gsl_vector * eval = gsl_vector_alloc(n);
 973       gsl_vector * evalv = gsl_vector_alloc(n);
 974       gsl_vector * x = gsl_vector_alloc(n);
 975       gsl_vector * y = gsl_vector_alloc(n);
 976       gsl_vector_complex * work = gsl_vector_complex_alloc(n);
 977       gsl_matrix_complex * evec = gsl_matrix_complex_alloc(n, n);
 978       gsl_eigen_genherm_workspace * w = gsl_eigen_genherm_alloc(n);
 979       gsl_eigen_genhermv_workspace * wv = gsl_eigen_genhermv_alloc(n);
 980
 981       for (i = 0; i < 5; ++i)
 982         {
 983           create_random_herm_matrix(A, r, -10, 10);
 984           create_random_complex_posdef_matrix(B, r, work);
 985
 986           gsl_matrix_complex_memcpy(ma, A);
 987           gsl_matrix_complex_memcpy(mb, B);
 988
 989           gsl_eigen_genhermv(ma, mb, evalv, evec, wv);
 990           test_eigen_genherm_results(A, B, evalv, evec, i, "random", "unsorted");
 991
 992           gsl_matrix_complex_memcpy(ma, A);
 993           gsl_matrix_complex_memcpy(mb, B);
 994
 995           gsl_eigen_genherm(ma, mb, eval, w);
 996
 997           /* eval and evalv have to be sorted? not sure why */
 998           gsl_vector_memcpy(x, eval);
 999           gsl_vector_memcpy(y, evalv);
1000           gsl_sort_vector(x);
1001           gsl_sort_vector(y);
1002           test_eigenvalues_real(y, x, "genherm, random", "unsorted");
1003
1004           gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
1005           test_eigen_genherm_results(A, B, evalv, evec, i, "random", "val/asc");
1006
1007           gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
1008           test_eigen_genherm_results(A, B, evalv, evec, i, "random", "val/desc");
1009
1010           gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
1011           test_eigen_genherm_results(A, B, evalv, evec, i, "random", "abs/asc");
1012           gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
1013           test_eigen_genherm_results(A, B, evalv, evec, i, "random", "abs/desc");
1014         }
1015
1016       gsl_matrix_complex_free(A);
1017       gsl_matrix_complex_free(B);
1018       gsl_matrix_complex_free(ma);
1019       gsl_matrix_complex_free(mb);
1020       gsl_vector_free(eval);
1021       gsl_vector_free(evalv);
1022       gsl_vector_free(x);
1023       gsl_vector_free(y);
1024       gsl_vector_complex_free(work);
1025       gsl_matrix_complex_free(evec);
1026       gsl_eigen_genherm_free(w);
1027       gsl_eigen_genhermv_free(wv);
1028     }
1029
1030   gsl_rng_free(r);
1031 } /* test_eigen_genherm() */
1032
1033 /******************************************
1034  * gen test code                          *
1035  ******************************************/
1036
1037 typedef struct
1038 {
1039   gsl_matrix *A;
1040   gsl_matrix *B;
1041   gsl_vector_complex *alpha;
1042   gsl_vector *beta;
1043   gsl_vector_complex *alphav;
1044   gsl_vector *betav;
1045   gsl_vector_complex *eval;
1046   gsl_vector_complex *evalv;
1047   gsl_vector *x;
1048   gsl_vector *y;
1049   gsl_matrix *Q;
1050   gsl_matrix *Z;
1051   gsl_matrix_complex *evec;
1052   gsl_eigen_gen_workspace *gen_p;
1053   gsl_eigen_genv_workspace *genv_p;
1054 } test_eigen_gen_workspace;
1055
1056 test_eigen_gen_workspace *
1057 test_eigen_gen_alloc(const size_t n)
1058 {
1059   test_eigen_gen_workspace *w;
1060
1061   w = calloc(1, sizeof(test_eigen_gen_workspace));
1062
1063   w->A = gsl_matrix_alloc(n, n);
1064   w->B = gsl_matrix_alloc(n, n);
1065   w->alpha = gsl_vector_complex_alloc(n);
1066   w->beta = gsl_vector_alloc(n);
1067   w->alphav = gsl_vector_complex_alloc(n);
1068   w->betav = gsl_vector_alloc(n);
1069   w->eval = gsl_vector_complex_alloc(n);
1070   w->evalv = gsl_vector_complex_alloc(n);
1071   w->x = gsl_vector_alloc(n);
1072   w->y = gsl_vector_alloc(n);
1073   w->Q = gsl_matrix_alloc(n, n);
1074   w->Z = gsl_matrix_alloc(n, n);
1075   w->evec = gsl_matrix_complex_alloc(n, n);
1076   w->gen_p = gsl_eigen_gen_alloc(n);
1077   w->genv_p = gsl_eigen_genv_alloc(n);
1078
1079   return (w);
1080 } /* test_eigen_gen_alloc() */
1081
1082 void
1083 test_eigen_gen_free(test_eigen_gen_workspace *w)
1084 {
1085   gsl_matrix_free(w->A);
1086   gsl_matrix_free(w->B);
1087   gsl_vector_complex_free(w->alpha);
1088   gsl_vector_free(w->beta);
1089   gsl_vector_complex_free(w->alphav);
1090   gsl_vector_free(w->betav);
1091   gsl_vector_complex_free(w->eval);
1092   gsl_vector_complex_free(w->evalv);
1093   gsl_vector_free(w->x);
1094   gsl_vector_free(w->y);
1095   gsl_matrix_free(w->Q);
1096   gsl_matrix_free(w->Z);
1097   gsl_matrix_complex_free(w->evec);
1098   gsl_eigen_gen_free(w->gen_p);
1099   gsl_eigen_genv_free(w->genv_p);
1100   free(w);
1101 } /* test_eigen_gen_free() */
1102
1103 void
1104 test_eigen_gen_results (const gsl_matrix * A, const gsl_matrix * B,
1105                         const gsl_vector_complex * alpha,
1106                         const gsl_vector * beta,
1107                         const gsl_matrix_complex * evec,
1108                         size_t count, const char * desc,
1109                         const char * desc2)
1110 {
1111   const size_t N = A->size1;
1112   size_t i, j;
1113   gsl_matrix_complex *ma, *mb;
1114   gsl_vector_complex *x, *y;
1115   gsl_complex z_one, z_zero;
1116
1117   ma = gsl_matrix_complex_alloc(N, N);
1118   mb = gsl_matrix_complex_alloc(N, N);
1119   y = gsl_vector_complex_alloc(N);
1120   x = gsl_vector_complex_alloc(N);
1121
1122   /* ma <- A, mb <- B */
1123   for (i = 0; i < N; ++i)
1124     {
1125       for (j = 0; j < N; ++j)
1126         {
1127           gsl_complex z;
1128
1129           GSL_SET_COMPLEX(&z, gsl_matrix_get(A, i, j), 0.0);
1130           gsl_matrix_complex_set(ma, i, j, z);
1131
1132           GSL_SET_COMPLEX(&z, gsl_matrix_get(B, i, j), 0.0);
1133           gsl_matrix_complex_set(mb, i, j, z);
1134         }
1135     }
1136
1137   GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
1138   GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
1139
1140   /* check eigenvalues */
1141   for (i = 0; i < N; ++i)
1142     {
1143       gsl_vector_complex_const_view vi =
1144         gsl_matrix_complex_const_column(evec, i);
1145       gsl_complex ai = gsl_vector_complex_get(alpha, i);
1146       double bi = gsl_vector_get(beta, i);
1147
1148       /* compute x = alpha * B * v */
1149       gsl_blas_zgemv(CblasNoTrans, z_one, mb, &vi.vector, z_zero, x);
1150       gsl_blas_zscal(ai, x);
1151
1152       /* compute y = beta * A v */
1153       gsl_blas_zgemv(CblasNoTrans, z_one, ma, &vi.vector, z_zero, y);
1154       gsl_blas_zdscal(bi, y);
1155
1156       /* now test if y = x */
1157       for (j = 0; j < N; ++j)
1158         {
1159           gsl_complex xj = gsl_vector_complex_get(x, j);
1160           gsl_complex yj = gsl_vector_complex_get(y, j);
1161
1162           gsl_test_abs(GSL_REAL(yj), GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
1163                        "gen(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s",
1164                        N, count, desc, i, j, desc2);
1165           gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
1166                        "gen(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s",
1167                        N, count, desc, i, j, desc2);
1168         }
1169     }
1170
1171   gsl_matrix_complex_free(ma);
1172   gsl_matrix_complex_free(mb);
1173   gsl_vector_complex_free(y);
1174   gsl_vector_complex_free(x);
1175 } /* test_eigen_gen_results() */
1176
1177 void
1178 test_eigen_gen_pencil(const gsl_matrix * A, const gsl_matrix * B,
1179                       size_t count, const char * desc, int test_schur,
1180                       test_eigen_gen_workspace *w)
1181 {
1182   const size_t N = A->size1;
1183   size_t i;
1184
1185   gsl_matrix_memcpy(w->A, A);
1186   gsl_matrix_memcpy(w->B, B);
1187
1188   if (test_schur)
1189     {
1190       gsl_eigen_genv_QZ(w->A, w->B, w->alphav, w->betav, w->evec, w->Q, w->Z, w->genv_p);
1191       test_eigen_schur(A, w->A, w->Q, w->Z, count, "genv/A", desc);
1192       test_eigen_schur(B, w->B, w->Q, w->Z, count, "genv/B", desc);
1193     }
1194   else
1195     gsl_eigen_genv(w->A, w->B, w->alphav, w->betav, w->evec, w->genv_p);
1196
1197   test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "unsorted");
1198
1199   gsl_matrix_memcpy(w->A, A);
1200   gsl_matrix_memcpy(w->B, B);
1201
1202   if (test_schur)
1203     {
1204       gsl_eigen_gen_params(1, 1, 0, w->gen_p);
1205       gsl_eigen_gen_QZ(w->A, w->B, w->alpha, w->beta, w->Q, w->Z, w->gen_p);
1206       test_eigen_schur(A, w->A, w->Q, w->Z, count, "gen/A", desc);
1207       test_eigen_schur(B, w->B, w->Q, w->Z, count, "gen/B", desc);
1208     }
1209   else
1210     {
1211       gsl_eigen_gen_params(0, 0, 0, w->gen_p);
1212       gsl_eigen_gen(w->A, w->B, w->alpha, w->beta, w->gen_p);
1213     }
1214
1215   /* compute eval = alpha / beta values */
1216   for (i = 0; i < N; ++i)
1217     {
1218       gsl_complex z, ai;
1219       double bi;
1220
1221       ai = gsl_vector_complex_get(w->alpha, i);
1222       bi = gsl_vector_get(w->beta, i);
1223       GSL_SET_COMPLEX(&z, GSL_REAL(ai) / bi, GSL_IMAG(ai) / bi);
1224       gsl_vector_complex_set(w->eval, i, z);
1225
1226       ai = gsl_vector_complex_get(w->alphav, i);
1227       bi = gsl_vector_get(w->betav, i);
1228       GSL_SET_COMPLEX(&z, GSL_REAL(ai) / bi, GSL_IMAG(ai) / bi);
1229       gsl_vector_complex_set(w->evalv, i, z);
1230     }
1231
1232   /* sort eval and evalv and test them */
1233   gsl_eigen_nonsymmv_sort(w->eval, NULL, GSL_EIGEN_SORT_ABS_ASC);
1234   gsl_eigen_nonsymmv_sort(w->evalv, NULL, GSL_EIGEN_SORT_ABS_ASC);
1235   test_eigenvalues_complex(w->evalv, w->eval, "gen", desc);
1236
1237   gsl_eigen_genv_sort(w->alphav, w->betav, w->evec, GSL_EIGEN_SORT_ABS_ASC);
1238   test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "abs/asc");
1239   gsl_eigen_genv_sort(w->alphav, w->betav, w->evec, GSL_EIGEN_SORT_ABS_DESC);
1240   test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "abs/desc");
1241 } /* test_eigen_gen_pencil() */
1242
1243 void
1244 test_eigen_gen(void)
1245 {
1246   size_t N_max = 20;
1247   size_t n, i;
1248   gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
1249
1250   for (n = 1; n <= N_max; ++n)
1251     {
1252       gsl_matrix * A = gsl_matrix_alloc(n, n);
1253       gsl_matrix * B = gsl_matrix_alloc(n, n);
1254       test_eigen_gen_workspace * w = test_eigen_gen_alloc(n);
1255
1256       for (i = 0; i < 5; ++i)
1257         {
1258           create_random_nonsymm_matrix(A, r, -10, 10);
1259           create_random_nonsymm_matrix(B, r, -10, 10);
1260
1261           test_eigen_gen_pencil(A, B, i, "random", 0, w);
1262           test_eigen_gen_pencil(A, B, i, "random", 1, w);
1263         }
1264
1265       gsl_matrix_free(A);
1266       gsl_matrix_free(B);
1267       test_eigen_gen_free(w);
1268     }
1269
1270   gsl_rng_free(r);
1271
1272   /* this system will test the exceptional shift code */
1273   {
1274     double datA[] = { 1, 1, 0,
1275                       0, 0, -1,
1276                       1, 0, 0 };
1277     double datB[] = { -1, 0, -1,
1278                       0, -1, 0,
1279                       0, 0, -1 };
1280     gsl_matrix_view va = gsl_matrix_view_array (datA, 3, 3);
1281     gsl_matrix_view vb = gsl_matrix_view_array (datB, 3, 3);
1282     test_eigen_gen_workspace * w = test_eigen_gen_alloc(3);
1283
1284     test_eigen_gen_pencil(&va.matrix, &vb.matrix, 0, "integer", 0, w);
1285     test_eigen_gen_pencil(&va.matrix, &vb.matrix, 0, "integer", 1, w);
1286
1287     test_eigen_gen_free(w);
1288   }
1289 } /* test_eigen_gen() */
1290
1291 int
1292 main()
1293 {
1294   gsl_ieee_env_setup ();
1295   gsl_rng_env_setup ();
1296
1297   test_eigen_symm();
1298   test_eigen_herm();
1299   test_eigen_nonsymm();
1300   test_eigen_gensymm();
1301   test_eigen_genherm();
1302   test_eigen_gen();
1303
1304   exit (gsl_test_summary());
1305 }