test_source.c

math library from gnu

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          k++;
          FUNCTION (gsl_matrix, set) (m, i, j, -(ATOMIC) (k % 10 + 1));
        }
    }

  {
    status = (FUNCTION(gsl_matrix,isnull)(m) != 0);
    TEST (status, "_isnull" DESC " on negative matrix") ;

    status = (FUNCTION(gsl_matrix,ispos)(m) != 0);
    TEST (status, "_ispos" DESC " on negative matrix") ;

    status = (FUNCTION(gsl_matrix,isneg)(m) != 1);
    TEST (status, "_isneg" DESC " on negative matrix") ;

    status = (FUNCTION(gsl_matrix,isnonneg)(m) != 0);
    TEST (status, "_isnonneg" DESC " on negative matrix") ;
  }
#endif

  {
    TYPE (gsl_matrix) * a = FUNCTION (gsl_matrix, calloc) (M, N);
    TYPE (gsl_matrix) * b = FUNCTION (gsl_matrix, calloc) (M, N);
    
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            FUNCTION (gsl_matrix, set) (a, i, j, (BASE)(3 + i +  5 * j));
            FUNCTION (gsl_matrix, set) (b, i, j, (BASE)(3 + 2 * i + 4 * j));
          }
      }
    
    FUNCTION(gsl_matrix, memcpy) (m, a);
    FUNCTION(gsl_matrix, add) (m, b);
    
    {
      int status = 0;
      
      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = FUNCTION(gsl_matrix,get) (b,i,j);
              BASE z = x + y;
              if (r != z)
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_add matrix addition");
    }


    FUNCTION(gsl_matrix, memcpy) (m, a);
    FUNCTION(gsl_matrix, sub) (m, b);
    
    {
      int status = 0;
      
      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = FUNCTION(gsl_matrix,get) (b,i,j);
              BASE z = x - y;
              if (r != z)
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_sub matrix subtraction");
    }

    FUNCTION(gsl_matrix, memcpy) (m, a);
    FUNCTION(gsl_matrix, mul_elements) (m, b);
    
    {
      int status = 0;
      
      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = FUNCTION(gsl_matrix,get) (b,i,j);
              BASE z = x * y;
              if (r != z)
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_mul_elements multiplication");
    }

    FUNCTION(gsl_matrix, memcpy) (m, a);
    FUNCTION(gsl_matrix, div_elements) (m, b);
    
    {
      int status = 0;
      
      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = FUNCTION(gsl_matrix,get) (b,i,j);
              BASE z = x / y;
              if (fabs(r - z) > 2 * GSL_FLT_EPSILON * fabs(z))
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_div_elements division");
    }


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

    {
      int status = 0;

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              if (r !=  (ATOMIC)(2*x))
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_scale");
    }

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

    {
      int status = 0;

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = x + 3.0;
              if (fabs(r - y) > 2 * GSL_FLT_EPSILON * fabs(y))
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_add_constant");
    }

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

    {
      int status = 0;

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE r = FUNCTION(gsl_matrix,get) (m,i,j);
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = (i == j) ? (x + 5.0) : x;
              if (fabs(r - y) > 2 * GSL_FLT_EPSILON * fabs(y))
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_add_diagonal");
    }


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

    {
      int status = 0;

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < N; j++)
            {
              BASE x = FUNCTION(gsl_matrix,get) (a,i,j);
              BASE y = FUNCTION(gsl_matrix,get) (b,i,j);
              if (y != (BASE)(3 + i +  5 * j) || x != (BASE)(3 + 2 * i + 4 * j))
                status = 1;
            }
        }
      gsl_test (status, NAME (gsl_matrix) "_swap");
    }
      

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

 FUNCTION (gsl_matrix, free) (m);
 FUNCTION (gsl_vector, free) (v);
}

#if !(USES_LONGDOUBLE && !HAVE_PRINTF_LONGDOUBLE)
void
FUNCTION (test, text) (void)
{
  TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (M, N);

  size_t i, j;
  int k = 0;

  {
    FILE *f = fopen ("test.txt", "w");

    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            FUNCTION (gsl_matrix, set) (m, i, j, (BASE) k);
          }
      }

    FUNCTION (gsl_matrix, fprintf) (f, m, OUT_FORMAT);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.txt", "r");
    TYPE (gsl_matrix) * mm = FUNCTION (gsl_matrix, alloc) (M, N);
    status = 0;

    FUNCTION (gsl_matrix, fscanf) (f, mm);
    k = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            if (mm->data[i * N + j] != (BASE) k)
              status = 1;
          }
      }

    gsl_test (status, NAME (gsl_matrix) "_fprintf and fscanf");

    fclose (f);
    FUNCTION (gsl_matrix, free) (mm);
  }

  FUNCTION (gsl_matrix, free) (m);
}
#endif

void
FUNCTION (test, binary) (void)
{
  TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, calloc) (M, N);

  size_t i, j;
  size_t k = 0;

  {
    FILE *f = fopen ("test.dat", "wb");
    k = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            FUNCTION (gsl_matrix, set) (m, i, j, (BASE) k);
          }
      }

    FUNCTION (gsl_matrix, fwrite) (f, m);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "rb");
    TYPE (gsl_matrix) * mm = FUNCTION (gsl_matrix, alloc) (M, N);
    status = 0;

    FUNCTION (gsl_matrix, fread) (f, mm);
    k = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            if (mm->data[i * N + j] != (BASE) k)
              status = 1;
          }
      }

    gsl_test (status, NAME (gsl_matrix) "_write and read");

    fclose (f);
    FUNCTION (gsl_matrix, free) (mm);
  }

  FUNCTION (gsl_matrix, free) (m);
}

void
FUNCTION (test, trap) (void)
{
  TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (M, N);

  size_t i = 0, j = 0;
  double x;

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

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

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

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

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, i - 1, 0);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 1st index below lower bound");
  gsl_test (x != 0,
     NAME (gsl_matrix) "_get returns zero for 1st index below lower bound");

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, 0, j - 1);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 2nd index below lower bound");
  gsl_test (x != 0,
     NAME (gsl_matrix) "_get returns zero for 2nd index below lower bound");

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, M + 1, 0);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 1st index above upper bound");
  gsl_test (x != 0,
     NAME (gsl_matrix) "_get returns zero for 1st index above upper bound");

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, 0, N + 1);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 2nd index above upper bound");
  gsl_test (x != 0,
     NAME (gsl_matrix) "_get returns zero for 2nd index above upper bound");

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, M, 0);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 1st index at upper bound");
  gsl_test (x != 0,
        NAME (gsl_matrix) "_get returns zero for 1st index at upper bound");

  status = 0;
  x = FUNCTION (gsl_matrix, get) (m, 0, N);
  gsl_test (!status,
            NAME (gsl_matrix) "_get traps 2nd index at upper bound");
  gsl_test (x != 0,
        NAME (gsl_matrix) "_get returns zero for 2nd index at upper bound");

  FUNCTION (gsl_matrix, free) (m);
}