Fshl0

Gaussian Elimination

Mar 3rd, 2021
575
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. const double EPS = 1e-9;
  2. const int INF = 2; // it does not have to be infinity or a big number
  3.  
  4. /*  - this function return the number of solutions of the system (0, 1, or infinity).
  5. if at least one solution exists, then it is returned in the vector Ans.
  6.     - The function uses two pointers - the current column "col" and the current
  7.     row "row"
  8.     - For each variable xi, the value where(i) is the line where this column is not
  9.     zero. This vector is needed because some variables can be independent.
  10.     - In this implementation, the current ith line is not divided by aii as
  11.     described above, so in the end the matrix is not identity matrix (though
  12.     apparently dividing the ith line can help reducing errors)
  13.     - After finding solution, it is inserted back into the matrix - to check
  14.     whether the system has at least one solution or not. if the test solution is
  15.     successful then the function returns 1 or inf, depending on whether there is
  16.     at least one independent variable.
  17. */
  18.  
  19. int GaussianElimination (vector <vector <double> > A,
  20.                          vector <double>& Ans) {
  21.     int n = (int) A.size();
  22.     int m = (int) A[0].size() - 1; //ignore the last column which is the vector of B;
  23.     vector <int> where (m, -1);
  24.     for (int col = 0, row = 0; col < m && row < n; col++) {
  25.         int sel = row;
  26.         for (int i = row; i < n; i++)
  27.             if (abs (A[i][col]) > abs (A[sel][col]))
  28.                 sel = i;
  29.         if (abs (A[sel][col]) < EPS) continue;
  30.         for (int i = col; i <= m; i++)
  31.             swap (A[sel][i], A[row][i]);
  32.         where[col] = row;
  33.         for (int i = 0; i < n; i++) {
  34.             if (i != row) {
  35.                 double c = A[i][col] / A[row][col];
  36.                 for (int j = col; j <= m; j++)
  37.                     A[i][j] -= A[row][j] * c;
  38.             }
  39.         }
  40.         row++;
  41.     }
  42.     Ans.assign (m, 0);
  43.     for (int i = 0; i < m; i++) {
  44.         if (where[i] != -1)
  45.             Ans[i] = A[where[i]][m] / A[where[i]][i];
  46.     }
  47.     for (int i = 0; i < n; i++) {
  48.         double sum = 0;
  49.         for (int j = 0; j < m; j++)
  50.             sum += Ans[j] * A[i][j];
  51.         if (abs (sum - A[i][m]) > EPS)
  52.             return 0;
  53.     }
  54.     return 1;
  55. }
RAW Paste Data