gaussJordan.h

Go to the documentation of this file.

 
#pragma once
#include "../../Matrix.h"
#include "../../matrix_utils.h"
#include "LU.h"
#include "backwardSub.h"
#include "forwardSub.h"
#include <vector>
 
Matrix<double> gaussJordan(const Matrix<double>& A) {
  size_t i, j, k;
  size_t n = A.columns();
  double d;
  Matrix<double> mat = HorizontalConcat(A, eye(A.rows(), n));
 
  // Partial pivoting
  for(i = n - 1; i > 1; i--) {
    if(mat(i - 1, 1) < mat(i, 1)) {
      for(j = 0; j < (2 * n); ++j) {
        d             = mat(i, j);
        mat(i, j)     = mat(i - 1, j);
        mat(i - 1, j) = d;
      }
    }
  }
 
  // Reducing To Diagonal Matrix
  for(i = 0; i < n; i++) {
    for(j = 0; j < n; j++) {
      if(j != i) {
        d = mat(j, i) / mat(i, i);
        for(k = 0; k < n * 2; ++k) { mat(j, k) -= mat(i, k) * d; }
      }
    }
    d = mat(i, i);
    for(j = 0; j < (2 * n); j++) { mat(i, j) = mat(i, j) / d; }
  }
 
  return mat.GetSlice(0, mat.rows() - 1, n, mat.columns() - 1);
}