Search  

---

C:/temp/src/j2k/Beta/Math/Vector/NotDone/Gauss/test_gaussj.cpp

Go to the documentation of this file.

#include <time.h>
#include <stdio.h>
#include <stdlib.h>

void main()
{
  int i, j, k, n = 2000, n2 = n * n, t0, t = 0;

  double *A = (double *) malloc(n2 * sizeof(double));

  double *a[4096], b[4096], x[4096], e, ee; a[0] = A;

  for (i = 1; i < n; i++) a[i] = a[i-1] + n;

  srand(clock());

  for (i = 0; i < n2; i++) A[i] = 0.0;

  for (i = 0; i < n; i++) a[i][i] = 1.0;

  for (k = n; k < 6000; k++)
    {
      do i = rand() * n / 32768, j = rand() * n / 32768; while (i <= j || a[i][j]);

      a[i][j] = (1 + rand()) / 32768.0;
    }

  for (k = 0; k < n2; k++)
    {
      i = rand() * n / 32768; j = rand() * n / 32768;

      e = a[i][j]; a[i][j] = a[j][i]; a[j][i] = e;
    }

  for (k = i = 0; i < n; i++) for (j = 0; j < n; j++) if (a[i][j]) k++;

  printf("k = %d\n", k);

  for (i = 0; i < n; i++) x[i] = rand() / 32767.0;

  for (i = 0; i < n; i++) for (b[i] = j = 0; j < n; j++) b[i] += a[i][j]*x[j];

  t0 = clock();
  a[n] = b; if (!gaussj(a, n, 1, 1e-7)) return;
  t += clock() - t0;

  for (e = i = 0; i < n; i++) if ((ee = fabs(x[i] - b[i])) > e) e = ee;

  printf("Error = %e, Time = %d\n", e, t);

  free(A);
}

---