test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

    double ai = actual[i];
    int foo = check(si, ai, eps);
    if(foo) {
      printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
    }
    s += foo;
  }

  gsl_vector_free(x);
  gsl_vector_free(rhs);
  gsl_vector_free(diag);
  gsl_vector_free(offdiag);

  return s;
}


int test_TD_solve(void)
{
  int f;
  int s = 0;
  double actual[16];

  actual[0] =  0.0;
  actual[1] =  2.0;
  f = test_TD_solve_dim(2, 1.0, 0.5, actual, 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TD dim=2 A");
  s += f;

  actual[0] =  3.0/8.0;
  actual[1] =  15.0/8.0;
  f = test_TD_solve_dim(2, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TD dim=2 B");
  s += f;

  actual[0] =  5.0/8.0;
  actual[1] =  9.0/8.0;
  actual[2] =  2.0;
  actual[3] =  15.0/8.0;
  actual[4] =  35.0/8.0;
  f = test_TD_solve_dim(5, 1.0, 1.0/3.0, actual, 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TD dim=5");
  s += f;

  return s;
}

int
test_TD_cyc_solve_one(const size_t dim, const double * d, const double * od,
                      const double * r, const double * actual, double eps)
{
  int s = 0;
  size_t i;

  gsl_vector * offdiag = gsl_vector_alloc(dim);
  gsl_vector * diag = gsl_vector_alloc(dim);
  gsl_vector * rhs = gsl_vector_alloc(dim);
  gsl_vector * x = gsl_vector_alloc(dim);

  for(i=0; i<dim; i++) {
    gsl_vector_set(diag, i, d[i]);
    gsl_vector_set(rhs,  i, r[i]);
    gsl_vector_set(offdiag, i, od[i]);
  }

  s += gsl_linalg_solve_symm_cyc_tridiag(diag, offdiag, rhs, x);

  for(i=0; i<dim; i++) {
    double si = gsl_vector_get(x, i);
    double ai = actual[i];
    int foo = check(si, ai, eps);
    if(foo) {
      printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
    }
    s += foo;
  }

  gsl_vector_free(x);
  gsl_vector_free(rhs);
  gsl_vector_free(diag);
  gsl_vector_free(offdiag);

  return s;
}

int test_TD_cyc_solve(void)
{
  int f;
  int s = 0;

  double diag_1[] = { 1, 1, 1 };
  double offdiag_1[] = { 3, 3, 3 };
  double rhs_1[] = { 7, -7, 7 };
  double actual_1[] = { -2, 5, -2 };
  size_t dim_1 = sizeof(diag_1)/sizeof(double);

  double diag_2[] = { 4, 2, 1, 2, 4 };
  double offdiag_2[] = { 1, 1, 1, 1, 1 };
  double rhs_2[] = { 30, -24, 3, 21, -30 };
  double actual_2[] = { 12, 3, -42, 42, -21 };
  size_t dim_2 = sizeof(diag_2)/sizeof(double);

  f = test_TD_cyc_solve_one(dim_1, diag_1, offdiag_1, rhs_1, actual_1, 28.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TD_cyc dim=%u A", dim_1);
  s += f;

/*  f = test_TD_cyc_solve_one(dim_2, diag_2, offdiag_2, rhs_2, actual_2, 7.0 * GSL_DBL_EPSILON);
  FIXME: bad accuracy */
  f = test_TD_cyc_solve_one(dim_2, diag_2, offdiag_2, rhs_2, actual_2, 35.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TD_cyc dim=%u B", dim_2);
  s += f;

  return s;
}

int
test_TDN_solve_dim(size_t dim, double d, double a, double b, const double * actual, double eps)
{
  int s = 0;
  size_t i;

  gsl_vector * abovediag = gsl_vector_alloc(dim-1);
  gsl_vector * belowdiag = gsl_vector_alloc(dim-1);
  gsl_vector * diag = gsl_vector_alloc(dim);
  gsl_vector * rhs = gsl_vector_alloc(dim);
  gsl_vector * x = gsl_vector_alloc(dim);

  for(i=0; i<dim; i++) {
    gsl_vector_set(diag, i, d);
    gsl_vector_set(rhs,  i, i + 1.0);
  }
  for(i=0; i<dim-1; i++) {
    gsl_vector_set(abovediag, i, a);
    gsl_vector_set(belowdiag, i, b);
  }

  s += gsl_linalg_solve_tridiag(diag, abovediag, belowdiag, rhs, x);

  for(i=0; i<dim; i++) {
    double si = gsl_vector_get(x, i);
    double ai = actual[i];
    int foo = check(si, ai, eps);
    if(foo) {
      printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
    }
    s += foo;
  }

  gsl_vector_free(x);
  gsl_vector_free(rhs);
  gsl_vector_free(diag);
  gsl_vector_free(abovediag);
  gsl_vector_free(belowdiag);

  return s;
}


int test_TDN_solve(void)
{
  int f;
  int s = 0;
  double actual[16];

  actual[0] =  3.0;
  actual[1] = -1.0;
  f = test_TDN_solve_dim(2, 1.0, 2.0, 1.0, actual, 1.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TDN dim=2 A");
  s += f;

  actual[0] =  2.0/5.0;
  actual[1] =  9.0/5.0;
  f = test_TDN_solve_dim(2, 1.0, 1.0/3.0, 1.0/2.0, actual, 1.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TDN dim=2 B");
  s += f;

  actual[0] =  99.0/140.0;
  actual[1] =  41.0/35.0;
  actual[2] =  19.0/10.0;
  actual[3] =  72.0/35.0;
  actual[4] =  139.0/35.0;
  f = test_TDN_solve_dim(5, 1.0, 1.0/4.0, 1.0/2.0, actual, 35.0/8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TDN dim=5");
  s += f;

  return s;
}

int
test_TDN_cyc_solve_dim(size_t dim, double d, double a, double b, const double * actual, double eps)
{
  int s = 0;
  size_t i;

  gsl_vector * abovediag = gsl_vector_alloc(dim);
  gsl_vector * belowdiag = gsl_vector_alloc(dim);
  gsl_vector * diag = gsl_vector_alloc(dim);
  gsl_vector * rhs = gsl_vector_alloc(dim);
  gsl_vector * x = gsl_vector_alloc(dim);

  for(i=0; i<dim; i++) {
    gsl_vector_set(diag, i, d);
    gsl_vector_set(rhs,  i, i + 1.0);
  }
  for(i=0; i<dim; i++) {
    gsl_vector_set(abovediag, i, a);
    gsl_vector_set(belowdiag, i, b);
  }

  s += gsl_linalg_solve_cyc_tridiag(diag, abovediag, belowdiag, rhs, x);

  for(i=0; i<dim; i++) {
    double si = gsl_vector_get(x, i);
    double ai = actual[i];
    int foo = check(si, ai, eps);
    if(foo) {
      printf("%3d[%d]: %22.18g   %22.18g\n", dim, i, gsl_vector_get(x, i), actual[i]);
    }
    s += foo;
  }

  gsl_vector_free(x);
  gsl_vector_free(rhs);
  gsl_vector_free(diag);
  gsl_vector_free(abovediag);
  gsl_vector_free(belowdiag);

  return s;
}


int test_TDN_cyc_solve(void)
{
  int f;
  int s = 0;
  double actual[16];

  actual[0] =  3.0/2.0;
  actual[1] = -1.0/2.0;
  actual[2] =  1.0/2.0;
  f = test_TDN_cyc_solve_dim(3, 1.0, 2.0, 1.0, actual, 4.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TDN_cyc dim=2 A");
  s += f;

  actual[0] = -5.0/22.0;
  actual[1] = -3.0/22.0;
  actual[2] =  29.0/22.0;
  actual[3] = -9.0/22.0;
  actual[4] =  43.0/22.0;
  f = test_TDN_cyc_solve_dim(5, 3.0, 2.0, 1.0, actual, 66.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  solve_TDN_cyc dim=5");
  s += f;

  return s;
}

int
test_bidiag_decomp_dim(const gsl_matrix * m, double eps)
{
  int s = 0;
  size_t i,j,k,r, M = m->size1, N = m->size2;

  gsl_matrix * A  = gsl_matrix_alloc(M,N);
  gsl_matrix * a  = gsl_matrix_alloc(M,N);
  gsl_matrix * b  = gsl_matrix_alloc(N,N);

  gsl_matrix * u  = gsl_matrix_alloc(M,N);
  gsl_matrix * v  = gsl_matrix_alloc(N,N);

  gsl_vector * tau1  = gsl_vector_alloc(N);
  gsl_vector * tau2  = gsl_vector_alloc(N-1);
  gsl_vector * d  = gsl_vector_alloc(N);
  gsl_vector * sd  = gsl_vector_alloc(N-1);

  gsl_matrix_memcpy(A,m);

  s += gsl_linalg_bidiag_decomp(A, tau1, tau2);
  s += gsl_linalg_bidiag_unpack(A, tau1, u, tau2, v, d, sd);

  gsl_matrix_set_zero(b);
  for (i = 0; i < N; i++) gsl_matrix_set(b, i,i, gsl_vector_get(d,i));
  for (i = 0; i < N-1; i++) gsl_matrix_set(b, i,i+1, gsl_vector_get(sd,i));
  
  /* Compute A = U B V^T */
  
  for (i = 0; i < M ; i++)
    {
      for (j = 0; j < N; j++)
        {
          double sum = 0;

          for (k = 0; k < N; k++)
            {
              for (r = 0; r < N; r++)
                {
                  sum += gsl_matrix_get(u, i, k) * gsl_matrix_get (b, k, r)
                    * gsl_matrix_get(v, j, r);
                }
            }
          gsl_matrix_set (a, i, j, sum);
        }
    }

  for(i=0; i<M; i++) {
    for(j=0; j<N; j++) {
      double aij = gsl_matrix_get(a, i, j);
      double mij = gsl_matrix_get(m, i, j);
      int foo = check(aij, mij, eps);
      if(foo) {
        printf("(%3d,%3d)[%d,%d]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
      }
      s += foo;
    }
  }

  gsl_matrix_free(A);
  gsl_matrix_free(a);
  gsl_matrix_free(u);
  gsl_matrix_free(v);
  gsl_matrix_free(b);
  gsl_vector_free(tau1);
  gsl_vector_free(tau2);
  gsl_vector_free(d);
  gsl_vector_free(sd);

  return s;
}

int test_bidiag_decomp(void)
{
  int f;
  int s = 0;

  f = test_bidiag_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp m(5,3)");
  s += f;

  f = test_bidiag_decomp_dim(m97, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp m(9,7)");
  s += f;

  f = test_bidiag_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp hilbert(2)");
  s += f;

  f = test_bidiag_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp hilbert(3)");
  s += f;

  f = test_bidiag_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp hilbert(4)");
  s += f;

  f = test_bidiag_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  bidiag_decomp hilbert(12)");
  s += f;

  return s;
}


int main(void)
{
  gsl_ieee_env_setup ();

  m35 = create_general_matrix(3,5);
  m53 = create_general_matrix(5,3);
  m97 = create_general_matrix(9,7);

  s35 = create_singular_matrix(3,5);
  s53 = create_singular_matrix(5,3);

  hilb2 = create_hilbert_matrix(2);
  hilb3 = create_hilbert_matrix(3);
  hilb4 = create_hilbert_matrix(4);
  hilb12 = create_hilbert_matrix(12);

  vander2 = create_vandermonde_matrix(2);
  vander3 = create_vandermonde_matrix(3);
  vander4 = create_vandermonde_matrix(4);
  vander12 = create_vandermonde_matrix(12);

  moler10 = create_moler_matrix(10);

  c7 = create_complex_matrix(7);

  row3 = create_row_matrix(3,3);
  row5 = create_row_matrix(5,5);
  row12 = create_row_matrix(12,12);

  A22 = create_2x2_matrix (0.0, 0.0, 0.0, 0.0);
  A33 = gsl_matrix_alloc(3,3);
  A44 = gsl_matrix_alloc(4,4);

  /* Matmult now obsolete */
#ifdef MATMULT
  gsl_test(test_matmult(),        "Matrix Multiply"); 
  gsl_test(test_matmult_mod(),    "Matrix Multiply with Modification"); 
#endif
  gsl_test(test_bidiag_decomp(),  "Bidiagonal Decomposition");
  gsl_test(test_LU_solve(),       "LU Decomposition and Solve");
  gsl_test(test_LUc_solve(),      "Complex LU Decomposition and Solve");
  gsl_test(test_QR_decomp(),      "QR Decomposition");
  gsl_test(test_QR_solve(),       "QR Solve");
  gsl_test(test_QR_QRsolve(),     "QR QR Solve");
  gsl_test(test_QR_lssolve(),     "QR LS Solve");
  gsl_test(test_QR_update(),      "QR Rank-1 Update");
  gsl_test(test_QRPT_decomp(),    "QRPT Decomposition");
  gsl_test(test_QRPT_solve(),     "QRPT Solve");
  gsl_test(test_QRPT_QRsolve(),   "QRPT QR Solve");
  gsl_test(test_SV_decomp(),      "Singular Value Decomposition");
  gsl_test(test_SV_solve(),       "SVD Solve");
  gsl_test(test_cholesky_decomp(),"Cholesky Decomposition");
  gsl_test(test_cholesky_solve(), "Cholesky Solve");
  gsl_test(test_HH_solve(),       "Householder solve");
  gsl_test(test_TD_solve(),       "Tridiagonal symmetric solve");
  gsl_test(test_TD_cyc_solve(),   "Tridiagonal symmetric cyclic solve");

#ifdef NONSYM
  gsl_test(test_TDN_solve(),      "Tridiagonal nonsymmetric solve");
  gsl_test(test_TDN_cyc_solve(),  "Tridiagonal nonsymmetric cyclic solve");
#endif

  gsl_matrix_free(m35);
  gsl_matrix_free(m53);
  gsl_matrix_free(m97);
  gsl_matrix_free(s35);
  gsl_matrix_free(s53);

  gsl_matrix_free(hilb2);
  gsl_matrix_free(hilb3);
  gsl_matrix_free(hilb4);
  gsl_matrix_free(hilb12);

  gsl_matrix_free(vander2);
  gsl_matrix_free(vander3);
  gsl_matrix_free(vander4);
  gsl_matrix_free(vander12);

  gsl_matrix_free(moler10);

  gsl_matrix_complex_free(c7);
  gsl_matrix_free(row3);
  gsl_matrix_free(row5);
  gsl_matrix_free(row12);

  gsl_matrix_free(A22);
  gsl_matrix_free(A33);
  gsl_matrix_free(A44);

  exit (gsl_test_summary());
}