gaussj.cpp File Reference

00004 {
00005   int i, j, k, d = 0; m += n;
00006 
00007   /* Forward elimination. */
00008   for (j = 0; j < n; j++)
00009     {
00010       int im = j, j1 = j + 1; double t, pivot = fabs(a[j][j]);
00011 
00012       /* Search for maximum coefficient in column. */
00013       for (i = j1; i < n; i++) if ((t = fabs(a[i][j])) > pivot) pivot = t, im = i;
00014 
00015       /* Test for pivot less than tolerance (singular matrix). */
00016       if (pivot <= tol) return 0;
00017 
00018       /* Interchange rows if necessary. */
00019       if (im != j)
00020    {
00021      double *t = a[j]; a[j] = a[im]; a[im] = t; d ^= 1;
00022 
00023      for (k = n; k < m; k++) { double t = a[k][im]; a[k][im] = a[k][j]; a[k][j] = t; }
00024    }
00025 
00026       /* Divide equation by leading coefficient. */
00027       pivot = 1.0 / a[j][j]; if (a[j][j] < 0.0) d ^= 1;
00028 
00029       for (k = j; k < n; k++) if (a[j][k]) a[j][k] *= pivot;
00030       for (k = n; k < m; k++) if (a[k][j]) a[k][j] *= pivot;
00031 
00032       /* Eliminate next variable. */
00033       for (i = j1; i < n; i++)
00034    if (t = a[i][j])
00035      {
00036        for (k = j1; k < n; k++) if (a[j][k]) a[i][k] -= t * a[j][k];
00037        for (k =  n; k < m; k++) if (a[k][j]) a[k][i] -= t * a[k][j];
00038      }
00039     }
00040 
00041   /* Backward substitution. */
00042   for (i = n - 2; i >= 0; i--)
00043     for (j = i + 1; j < n; j++)
00044       if (a[i][j]) for (k = n; k < m; k++) if (a[k][j]) a[k][i] -= a[i][j] * a[k][j];
00045 
00046   return d ? -1 : 1;
00047 }
C:/temp/src/j2k/Beta/Math/Vector/NotDone/Gauss/gaussj.cpp File Reference

Functions

Function Documentation
