Search  

---

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

Go to the documentation of this file.

#include <math.h>

int gaussj(double **a, int n, int m, double tol)
{
  int i, j, k, d = 0; m += n;

  /* Forward elimination. */
  for (j = 0; j < n; j++)
    {
      int im = j, j1 = j + 1; double t, pivot = fabs(a[j][j]);

      /* Search for maximum coefficient in column. */
      for (i = j1; i < n; i++) if ((t = fabs(a[i][j])) > pivot) pivot = t, im = i;

      /* Test for pivot less than tolerance (singular matrix). */
      if (pivot <= tol) return 0;

      /* Interchange rows if necessary. */
      if (im != j)
   {
     double *t = a[j]; a[j] = a[im]; a[im] = t; d ^= 1;

     for (k = n; k < m; k++) { double t = a[k][im]; a[k][im] = a[k][j]; a[k][j] = t; }
   }

      /* Divide equation by leading coefficient. */
      pivot = 1.0 / a[j][j]; if (a[j][j] < 0.0) d ^= 1;

      for (k = j; k < n; k++) if (a[j][k]) a[j][k] *= pivot;
      for (k = n; k < m; k++) if (a[k][j]) a[k][j] *= pivot;

      /* Eliminate next variable. */
      for (i = j1; i < n; i++)
   if (t = a[i][j])
     {
       for (k = j1; k < n; k++) if (a[j][k]) a[i][k] -= t * a[j][k];
       for (k =  n; k < m; k++) if (a[k][j]) a[k][i] -= t * a[k][j];
     }
    }

  /* Backward substitution. */
  for (i = n - 2; i >= 0; i--)
    for (j = i + 1; j < n; j++)
      if (a[i][j]) for (k = n; k < m; k++) if (a[k][j]) a[k][i] -= a[i][j] * a[k][j];

  return d ? -1 : 1;
}

---