Search  

---

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

Go to the documentation of this file.

#include <time.h>
#include <iostream>
#include "Ludecomp.h"

size_t smatrix<double>::PrintWidth = 8;
size_t smatrix<double>::PrintPrecision = 3;

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

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

  double *a[4096], x[4096], e; 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;

  smatrix<double> sA(A, n, n), sx(x, n, 1), sb = sA * sx;

  t0 = clock();
  TLUDecomp<double> LU(sA, true); sb = LU.Solve(sb);
  t += clock() - t0;

  e = (sb - sx).maxCoef;

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

  free(A);
}

---