test_complex_source.c

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

void
FUNCTION (test, trap) (void)
{
  TYPE (gsl_matrix) * mc = FUNCTION (gsl_matrix, alloc) (M, N);
  size_t i = 0, j = 0;

  BASE z = { {(ATOMIC) 1.2, (ATOMIC) 3.4} };
  BASE z1;

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, i - 1, j, z);
  gsl_test (!status,
            NAME (gsl_matrix) "_set traps 1st index below lower bound");

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, i, j - 1, z);
  gsl_test (!status,
            NAME (gsl_matrix) "_set traps 2nd index below lower bound");

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, M + 1, 0, z);
  gsl_test (!status,
            NAME (gsl_matrix) "_set traps 1st index above upper bound");

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, 0, N + 1, z);
  gsl_test (!status,
            NAME (gsl_matrix) "_set traps 2nd index above upper bound");

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, M, 0, z);
  gsl_test (!status, NAME (gsl_matrix) "_set traps 1st index at upper bound");

  status = 0;
  FUNCTION (gsl_matrix, set) (mc, 0, N, z);
  gsl_test (!status, NAME (gsl_matrix) "_set traps 2nd index at upper bound");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, i - 1, 0);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 1st index below lower bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 1st index below l.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 1st index below l.b.");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, 0, j - 1);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 2nd index below lower bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 2nd index below l.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 2nd index below l.b.");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, M + 1, 0);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 1st index above upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 1st index above u.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 1st index above u.b.");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, 0, N + 1);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 2nd index above upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 2nd index above u.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 2nd index above u.b.");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, M, 0);
  gsl_test (!status, NAME (gsl_matrix) "_get traps 1st index at upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 1st index at u.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 1st index at u.b.");

  status = 0;
  z1 = FUNCTION (gsl_matrix, get) (mc, 0, N);
  gsl_test (!status, NAME (gsl_matrix) "_get traps 2nd index at upper bound");
  gsl_test (GSL_REAL (z1) != 0,
            NAME (gsl_matrix) "_get, zero real for 2nd index at u.b.");
  gsl_test (GSL_IMAG (z1) != 0,
            NAME (gsl_matrix) "_get, zero imag for 2nd index at u.b.");

  FUNCTION (gsl_matrix, free) (mc);
}


void
FUNCTION (test, arith) (void)
{

#define P 8
#define Q 12
/* Must use smaller dimensions to prevent approximation of floats in float_mul_elements test*/

  TYPE (gsl_matrix) * a = FUNCTION (gsl_matrix, alloc) (P, Q);
  TYPE (gsl_matrix) * b = FUNCTION (gsl_matrix, alloc) (P, Q);
  TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (P, Q);
  size_t i, j;
  size_t k = 0;

  size_t status = 0;

  for (i = 0; i < P; i++)
    {
      for (j = 0; j < Q; j++)
        {
          BASE z, z1;
          GSL_REAL (z) = (ATOMIC) k;
          GSL_IMAG (z) = (ATOMIC) (k + 10);
          GSL_REAL (z1) = (ATOMIC) (k + 5);
          GSL_IMAG (z1) = (ATOMIC) (k + 20);

          FUNCTION (gsl_matrix, set) (a, i, j, z);
          FUNCTION (gsl_matrix, set) (b, i, j, z1);
          k++;
        }
    }

  {
    FUNCTION (gsl_matrix, memcpy) (m, a);

    FUNCTION (gsl_matrix, add) (m, b);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (GSL_REAL (z) != (ATOMIC) (2 * k + 5) ||
                GSL_IMAG (z) != (ATOMIC) (2 * k + 30))
              status = 1;
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_add matrix addition");
  }

  {
    FUNCTION (gsl_matrix, memcpy) (m, a);

    FUNCTION (gsl_matrix, sub) (m, b);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (GSL_REAL (z) != (ATOMIC) (-5)
                || GSL_IMAG (z) != (ATOMIC) (-10))
              status = 1;
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_sub matrix subtraction");
  }

  {
    FUNCTION (gsl_matrix, memcpy) (m, a);

    FUNCTION (gsl_matrix, mul_elements) (m, b);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            ATOMIC real = -(ATOMIC) (25 * k + 200);
            ATOMIC imag = (ATOMIC) (2 * k * k + 35 * k + 50);
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (fabs (GSL_REAL (z) - real) > 100 * BASE_EPSILON ||
                fabs (GSL_IMAG (z) - imag) > 100 * BASE_EPSILON)
              {
                status = 1;
#ifdef DEBUG
                printf ("%d\t%d\n", i, j);
                printf (OUT_FORMAT "\n",
                        GSL_REAL (z) + (ATOMIC) (25 * (ATOMIC) k + 200));
                printf (OUT_FORMAT "\n",
                        GSL_IMAG (z) - (ATOMIC) (2 * k * k + 35 * k + 50));
                printf ("\n");
#endif
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_mul_elements multiplication");
  }


  {
    FUNCTION (gsl_matrix, memcpy) (m, a);

    FUNCTION (gsl_matrix, div_elements) (m, b);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            ATOMIC denom = (2 * k * k + 50 * k + 425);
            ATOMIC real = (ATOMIC) (2 * k * k + 35 * k + 200) / denom;
            ATOMIC imag = ((ATOMIC) (50) - (ATOMIC) (5 * k)) / denom;
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (fabs (GSL_REAL (z) - real) > 100 * BASE_EPSILON ||
                fabs (GSL_IMAG (z) - imag) > 100 * BASE_EPSILON)
              {
#ifdef DEBUG
                printf (OUT_FORMAT "\t",
                        GSL_REAL (z) - (ATOMIC) (2 * k * k + 35 * k +
                                                 200) / denom);
                printf (OUT_FORMAT "\n",
                        GSL_IMAG (z) - ((ATOMIC) (50) -
                                        (ATOMIC) (5 * k)) / denom);
#endif
                status = 1;
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_div_elements division");
  }

  {
    BASE s;
    GSL_SET_COMPLEX(&s, 2.0, 3.0);

    FUNCTION (gsl_matrix, memcpy) (m, a);
    FUNCTION (gsl_matrix, scale) (m, s);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            ATOMIC real = (ATOMIC) (-(ATOMIC)k - 30);
            ATOMIC imag = (ATOMIC) (5 * (ATOMIC)k + 20);
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (GSL_REAL (z) != real || GSL_IMAG (z) != imag)
              {
                status = 1;
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_scale");
  }

  {
    BASE s;
    GSL_SET_COMPLEX(&s, 2.0, 3.0);

    FUNCTION (gsl_matrix, memcpy) (m, a);
    FUNCTION (gsl_matrix, add_constant) (m, s);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            ATOMIC real = (ATOMIC) ((ATOMIC)k + 2);
            ATOMIC imag = (ATOMIC) ((ATOMIC)k + 10 + 3);
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (GSL_REAL (z) != real || GSL_IMAG (z) != imag)
              {
                status = 1;
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_add_constant");
  }

  {
    BASE s;
    GSL_SET_COMPLEX(&s, 2.0, 3.0);

    FUNCTION (gsl_matrix, memcpy) (m, a);
    FUNCTION (gsl_matrix, add_diagonal) (m, s);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            ATOMIC real = (ATOMIC) ((ATOMIC)k + ((i==j) ? 2 : 0));
            ATOMIC imag = (ATOMIC) ((ATOMIC)k + 10 +((i==j) ? 3 : 0));
            BASE z = FUNCTION (gsl_matrix, get) (m, i, j);
            if (GSL_REAL (z) != real || GSL_IMAG (z) != imag)
              {
                status = 1;
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_add_diagonal");
  }

  {
    FUNCTION (gsl_matrix, swap) (a, b);

    k = 0;
    status = 0;

    for (i = 0; i < P; i++)
      {
        for (j = 0; j < Q; j++)
          {
            BASE x = FUNCTION (gsl_matrix, get) (a, i, j);
            BASE y = FUNCTION (gsl_matrix, get) (b, i, j);
            if (GSL_REAL (x) != (ATOMIC) (k + 5) || GSL_IMAG (x) != (ATOMIC) (k + 20) ||
                GSL_REAL (y) != (ATOMIC) (k) || GSL_IMAG (y) != (ATOMIC) (k + 10))
              {
                status = 1;
              }
            k++;
          }
      }
    gsl_test (status, NAME (gsl_matrix) "_swap");
  }



  FUNCTION (gsl_matrix, free) (a);
  FUNCTION (gsl_matrix, free) (b);
  FUNCTION (gsl_matrix, free) (m);

}