test.c

来自「math library from gnu」· C语言 代码 · 共 1,306 行 · 第 1/3 页

	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,
	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,
	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,
	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,
	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
    };
        
    gsl_matrix_view m;
    m = gsl_matrix_view_array (dat, 27, 27);
    test_eigen_symm_matrix(&m.matrix, 0, "symm(27)");
  };

} /* test_eigen_symm() */

/******************************************
 * herm test code                         *
 ******************************************/

void
test_eigen_herm_results (const gsl_matrix_complex * A, 
                         const gsl_vector * eval, 
                         const gsl_matrix_complex * evec, 
                         size_t count,
                         const char * desc,
                         const char * desc2)
{
  const size_t N = A->size1;
  size_t i, j;

  gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  gsl_vector_complex * y = gsl_vector_complex_alloc(N);

  /* check eigenvalues */

  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      gsl_vector_complex_const_view vi =
        gsl_matrix_complex_const_column(evec, i);
      gsl_vector_complex_memcpy(x, &vi.vector);
      /* compute y = m x (should = lambda v) */
      gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x, 
                      GSL_COMPLEX_ZERO, y);
      for (j = 0; j < N; j++)
        {
          gsl_complex xj = gsl_vector_complex_get (x, j);
          gsl_complex yj = gsl_vector_complex_get (y, j);
          gsl_test_rel(GSL_REAL(yj), ei * GSL_REAL(xj), 1e8*GSL_DBL_EPSILON, 
                       "%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2);
          gsl_test_rel(GSL_IMAG(yj), ei * GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON, 
                       "%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2);
        }
    }

  /* check eigenvectors are orthonormal */

  for (i = 0; i < N; i++)
    {
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      double nrm_v = gsl_blas_dznrm2(&vi.vector);
      gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s", 
                    desc, i, desc2);
    }

  for (i = 0; i < N; i++)
    {
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      for (j = i + 1; j < N; j++)
        {
          gsl_vector_complex_const_view vj 
            = gsl_matrix_complex_const_column(evec, j);
          gsl_complex vivj;
          gsl_blas_zdotc (&vi.vector, &vj.vector, &vivj);
          gsl_test_abs (gsl_complex_abs(vivj), 0.0, 10.0 * N * GSL_DBL_EPSILON, 
                        "%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
        }
    }

  gsl_vector_complex_free(x);
  gsl_vector_complex_free(y);
} /* test_eigen_herm_results() */

void
test_eigen_herm_matrix(const gsl_matrix_complex * m, size_t count,
                       const char * desc)
{
  const size_t N = m->size1;
  gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
  gsl_vector * eval = gsl_vector_alloc(N);
  gsl_vector * evalv = gsl_vector_alloc(N);
  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * y = gsl_vector_alloc(N);
  gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
  gsl_eigen_herm_workspace * w = gsl_eigen_herm_alloc(N);
  gsl_eigen_hermv_workspace * wv = gsl_eigen_hermv_alloc(N);

  gsl_matrix_complex_memcpy(A, m);

  gsl_eigen_hermv(A, evalv, evec, wv);
  test_eigen_herm_results(m, evalv, evec, count, desc, "unsorted");

  gsl_matrix_complex_memcpy(A, m);

  gsl_eigen_herm(A, eval, w);

  /* sort eval and evalv */
  gsl_vector_memcpy(x, eval);
  gsl_vector_memcpy(y, evalv);
  gsl_sort_vector(x);
  gsl_sort_vector(y);
  test_eigenvalues_real(y, x, desc, "unsorted");

  gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
  test_eigen_herm_results(m, evalv, evec, count, desc, "val/asc");

  gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
  test_eigen_herm_results(m, evalv, evec, count, desc, "val/desc");

  gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
  test_eigen_herm_results(m, evalv, evec, count, desc, "abs/asc");

  gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
  test_eigen_herm_results(m, evalv, evec, count, desc, "abs/desc");

  gsl_matrix_complex_free(A);
  gsl_vector_free(eval);
  gsl_vector_free(evalv);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_complex_free(evec);
  gsl_eigen_herm_free(w);
  gsl_eigen_hermv_free(wv);
} /* test_eigen_herm_matrix() */

void
test_eigen_herm(void)
{
  size_t N_max = 20;
  size_t n, i;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  for (n = 1; n <= N_max; ++n)
    {
      gsl_matrix_complex * A = gsl_matrix_complex_alloc(n, n);

      for (i = 0; i < 5; ++i)
        {
          create_random_herm_matrix(A, r, -10, 10);
          test_eigen_herm_matrix(A, i, "herm random");
        }

      gsl_matrix_complex_free(A);
    }

  gsl_rng_free(r);

  {
    double dat1[] =  { 0,0,  0,0, -1,0,  0,0, 
                       0,0,  1,0,  0,0,  1,0,
                       -1,0,  0,0,  0,0,  0,0,
                       0,0,  1,0,  0,0,  0,0 };
    double dat2[] =  { 1,0,  0,0, 0,0,  0,0, 
                       0,0,  2,0, 0,0,  0,0,
                       0,0,  0,0, 3,0,  0,0,
                       0,0,  0,0, 0,0,  4,0 };
    gsl_matrix_complex_view m;
    
    m = gsl_matrix_complex_view_array (dat1, 4, 4);
    test_eigen_herm_matrix(&m.matrix, 0, "herm(4)");

    m = gsl_matrix_complex_view_array (dat2, 4, 4);
    test_eigen_herm_matrix(&m.matrix, 1, "herm(4) diag");
  }
} /* test_eigen_herm() */

/******************************************
 * nonsymm test code                      *
 ******************************************/

void
test_eigen_nonsymm_results (const gsl_matrix * m, 
                            const gsl_vector_complex * eval, 
                            const gsl_matrix_complex * evec, 
                            size_t count,
                            const char * desc,
                            const char * desc2)
{
  size_t i,j;
  size_t N = m->size1;

  gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  gsl_vector_complex * y = gsl_vector_complex_alloc(N);
  gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);

  /* we need a complex matrix for the blas routines, so copy m into A */
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex z;
          GSL_SET_COMPLEX(&z, gsl_matrix_get(m, i, j), 0.0);
          gsl_matrix_complex_set(A, i, j, z);
        }
    }

  for (i = 0; i < N; i++)
    {
      gsl_complex ei = gsl_vector_complex_get (eval, i);
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      double norm = gsl_blas_dznrm2(&vi.vector);

      /* check that eigenvector is normalized */
      gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
                   "nonsymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count, desc, i, desc2);

      gsl_vector_complex_memcpy(x, &vi.vector);

      /* compute y = m x (should = lambda v) */
      gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x, 
                      GSL_COMPLEX_ZERO, y);

      /* compute x = lambda v */
      gsl_blas_zscal(ei, x);

      /* now test if y = x */
      for (j = 0; j < N; j++)
        {
          gsl_complex xj = gsl_vector_complex_get (x, j);
          gsl_complex yj = gsl_vector_complex_get (y, j);

          /* use abs here in case the values are close to 0 */
          gsl_test_abs(GSL_REAL(yj), GSL_REAL(xj), 1e8*GSL_DBL_EPSILON, 
                       "nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
          gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON, 
                       "nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
        }
    }

  gsl_matrix_complex_free(A);
  gsl_vector_complex_free(x);
  gsl_vector_complex_free(y);
}

void
test_eigen_nonsymm_matrix(const gsl_matrix * m, size_t count,
                          const char * desc,
                          gsl_eigen_nonsymmv_workspace *w)
{
  const size_t N = m->size1;
  gsl_matrix * A = gsl_matrix_alloc(N, N);
  gsl_matrix * Z = gsl_matrix_alloc(N, N);
  gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
  gsl_vector_complex * eval = gsl_vector_complex_alloc(N);

  /*
   * calculate eigenvalues and eigenvectors - it is sufficient to
   * test gsl_eigen_nonsymmv() since that function calls
   * gsl_eigen_nonsymm() for the eigenvalues
   */ 
  gsl_matrix_memcpy(A, m);
  gsl_eigen_nonsymmv(A, eval, evec, w);
  test_eigen_nonsymm_results (m, eval, evec, count, desc, "unsorted");

  /* test sort routines */
  gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
  test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/asc");

  gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);
  test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/desc");

  /* test Schur vectors */
  gsl_matrix_memcpy(A, m);
  gsl_eigen_nonsymmv_Z(A, eval, evec, Z, w);
  gsl_linalg_hessenberg_set_zero(A);
  test_eigen_schur(m, A, Z, Z, count, "nonsymm", desc);

  gsl_matrix_free(A);
  gsl_matrix_free(Z);
  gsl_matrix_complex_free(evec);
  gsl_vector_complex_free(eval);
}

void
test_eigen_nonsymm(void)
{
  size_t N_max = 20;
  size_t n, i;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  for (n = 1; n <= N_max; ++n)
    {
      gsl_matrix * m = gsl_matrix_alloc(n, n);
      gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(n);

      for (i = 0; i < 5; ++i)
        {
          create_random_nonsymm_matrix(m, r, -10, 10);

          gsl_eigen_nonsymm_params(1, 0, w->nonsymm_workspace_p);
          test_eigen_nonsymm_matrix(m, i, "random, unbalanced", w);

          gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);
          test_eigen_nonsymm_matrix(m, i, "random, balanced", w);
        }

      gsl_matrix_free(m);
      gsl_eigen_nonsymmv_free(w);
    }

  gsl_rng_free(r);

  {
    double dat1[] = { 0, 1, 1, 1,
                      1, 1, 1, 1,
                      0, 0, 0, 0,
                      0, 0, 0, 0 };
    double dat2[] = { 1, 1, 0, 1,
                      1, 1, 1, 1,
                      1, 1, 1, 1,
                      0, 1, 0, 0 };
    gsl_matrix_view v;
    gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(4);
    
    v = gsl_matrix_view_array (dat1, 4, 4);
    test_eigen_nonsymm_matrix(&v.matrix, 0, "integer", w);

    v = gsl_matrix_view_array (dat2, 4, 4);
    test_eigen_nonsymm_matrix(&v.matrix, 1, "integer", w);

    gsl_eigen_nonsymmv_free(w);
  }
} /* test_eigen_nonsymm() */

/******************************************
 * gensymm test code                      *
 ******************************************/

void
test_eigen_gensymm_results (const gsl_matrix * A, 
                            const gsl_matrix * B,
                            const gsl_vector * eval, 
                            const gsl_matrix * evec, 
                            size_t count,
                            const char * desc,
                            const char * desc2)
{
  const size_t N = A->size1;
  size_t i, j;

  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * y = gsl_vector_alloc(N);
  gsl_vector * z = gsl_vector_alloc(N);

  /* check A v = lambda B v */
  for (i = 0; i < N; i++)
    {
      double ei = gsl_vector_get (eval, i);
      gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
      double norm = gsl_blas_dnrm2(&vi.vector);

      /* check that eigenvector is normalized */
      gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
                   "gensymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
                   desc, i, desc2);

      gsl_vector_memcpy(z, &vi.vector);

      /* compute y = A z */
      gsl_blas_dgemv (CblasNoTrans, 1.0, A, z, 0.0, y);

      /* compute x = B z */
      gsl_blas_dgemv (CblasNoTrans, 1.0, B, z, 0.0, x);

      /* compute x = lambda B z */
      gsl_blas_dscal(ei, x);

      /* now test if y = x */
      for (j = 0; j < N; j++)
        {
          double xj = gsl_vector_get (x, j);
          double yj = gsl_vector_get (y, j);

          gsl_test_rel(yj, xj, 1e9 * GSL_DBL_EPSILON, 
                       "gensymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
        }
    }

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
}

void
test_eigen_gensymm(void)
{
  size_t N_max = 20;
  size_t n, i;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  for (n = 1; n <= N_max; ++n)
    {
      gsl_matrix * A = gsl_matrix_alloc(n, n);
      gsl_matrix * B = gsl_matrix_alloc(n, n);
      gsl_matrix * ma = gsl_matrix_alloc(n, n);
      gsl_matrix * mb = gsl_matrix_alloc(n, n);
      gsl_vector * eval = gsl_vector_alloc(n);
      gsl_vector * evalv = gsl_vector_alloc(n);
      gsl_vector * x = gsl_vector_alloc(n);
      gsl_vector * y = gsl_vector_alloc(n);
      gsl_matrix * evec = gsl_matrix_alloc(n, n);
      gsl_eigen_gensymm_workspace * w = gsl_eigen_gensymm_alloc(n);
      gsl_eigen_gensymmv_workspace * wv = gsl_eigen_gensymmv_alloc(n);

      for (i = 0; i < 5; ++i)
        {
          create_random_symm_matrix(A, r, -10, 10);
          create_random_posdef_matrix(B, r);

          gsl_matrix_memcpy(ma, A);
          gsl_matrix_memcpy(mb, B);

          gsl_eigen_gensymmv(ma, mb, evalv, evec, wv);
          test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "unsorted");

          gsl_matrix_memcpy(ma, A);
          gsl_matrix_memcpy(mb, B);

          gsl_eigen_gensymm(ma, mb, eval, w);

          /* eval and evalv have to be sorted? not sure why */
          gsl_vector_memcpy(x, eval);
          gsl_vector_memcpy(y, evalv);
          gsl_sort_vector(x);
          gsl_sort_vector(y);
          test_eigenvalues_real(y, x, "gensymm, random", "unsorted");