Calculate the inverse of an n x n matrix

239 Views Asked by At

I want to find the inverse of an n-order matrix using the method of finding the adjugate matrix and determinant, and then deducing the inverse. I'm using C++ language, and for calculating the determinant, I'm utilizing the concept of recursion. I want to know if there is a better way.

p.s. Due to the constraints of the assignment, other libraries are not allowed.

Matrix.h

class Matrix
{
private:
  friend ostream &operator<<(ostream &, const Matrix &);
  friend istream &operator>>(istream &, Matrix &);

  friend Matrix getADJ(const Matrix);
  friend float getDet(const Matrix);

public:
  Matrix();
  Matrix(const int, const int, const float = 0.0f);
  Matrix(const Matrix &);
  ~Matrix();

  // operation
  Matrix operator/(const float) const;
  Matrix &operator=(const Matrix &);
  Matrix inverse() const;

private:
  float **data;
  int row;
  int column;
};

Matrix.cpp

#include "Matrix.h"

// self define func.
static float getRandFloat(const float lower = 0.0F, const float upper = 10.0F)
{
  float temp = (float)rand() / RAND_MAX;
  return lower + temp * (upper - lower);
}
static float **new_2DArray(const int _r, const int _c, const float init_val = 0.0f)
{
  if (_r < 0 || _c < 0)
  {
    cout << "[ERROR] Size should be a positive number!!" << endl;

    return nullptr;
  }
  const unsigned long r_ul = (unsigned long)_r;
  const unsigned long c_ul = (unsigned long)_c;

  float **temp = new float *[r_ul];
  for (unsigned long r = 0; r < r_ul; r++)
    temp[r] = new float[c_ul];

  for (unsigned long r = 0; r < r_ul; r++)
    for (unsigned long c = 0; c < c_ul; c++)
      temp[r][c] = init_val;

  return temp;
}

// matrix constructor
Matrix::Matrix()
{
  data = new_2DArray(row = 2, column = 2);
  // set row to 2, column to 2
  // all values in matrix are set to default value 0.0f
}
Matrix::Matrix(const int row, const int col, const float init_val)
{
  data = new_2DArray(this->row = row, this->column = col, init_val);
}
Matrix::Matrix(const Matrix &obj)
{
  row = obj.row;
  column = obj.column;
  data = new_2DArray(row, column);
  // the new matrix has the same value for row and column
  // create a new data space for the new matrix, assign all value equals to matrix obj

  for (int r = 0; r < row; r++)
    for (int c = 0; c < column; c++)
      data[r][c] = obj.data[r][c];
}
Matrix::~Matrix()
{
  for (int i = 0; i < row; i++)
    delete[] data[i];

  delete data;
}

//operators
Matrix &Matrix::operator=(const Matrix &obj)
{
  if (row != obj.row || column != obj.column)
  {
    cout << "[ERROR] Invalid operation: Size different!!" << endl;
    return *this;
  }
  Matrix *temp = new Matrix(obj);

  data = temp->data;

  return *this;
}

Matrix Matrix::inverse() const
{
  if (row != column || row <= 1)
  {
    cout << "[ERROR] Invalid operation: should be a nxn matrix(n >= 2)!!" << endl;
    return *this;
  }
  Matrix temp(*this);
  float det = getDet(temp);
  if (det == 0.0f)
  {
    cout << "[ERROR] DET should`t be zero!!" << endl;
    return *this;
  }

  return getADJ(temp) / det;
};

Matrix getADJ(const Matrix m)
{
  const int n = m.row;
  Matrix result(n, n);

  for (int c = 0; c < n; c++)
  {
    const int limited_column = c;
    bool if_positive = (c % 2 ? false : true);
    for (int r = 0; r < n; r++)
    {
      const int limited_row = r;
      Matrix temp(n - 1, n - 1);

      for (int sr = 0; sr < n; sr++)
      {
        if (sr == limited_row)
          continue;
        for (int sc = 0; sc < n; sc++)
        {
          if (sc == limited_column)
            continue;

          temp.data[sr > limited_row ? sr - 1 : sr][sc > limited_column ? sc - 1 : sc] = m.data[sr][sc];
        }
      }
      result.data[c][r] = getDet(temp) * (if_positive ? 1.0f : -1.0f);
      if_positive = !if_positive;
    }
  }

  return result;
}

float getDet(const Matrix m)
{
  const int n = m.row;
  if (n == 1)
    return m.data[0][0];

  float curLevel = 0.0f;
  for (int i = 0; i < n; i++)
  {
    const int limited_column = i;
    // limit row is always equal to zero, when calculating DET
    Matrix temp(n - 1, n - 1);
    for (int r = 1; r < n; r++)
    {
      for (int c = 0; c < n; c++)
      {
        if (c == limited_column)
          continue;
        temp.data[r - 1][c > limited_column ? c - 1 : c] = m.data[r][c];
      }

      curLevel += (m.data[0][i] * getDet(temp) * (i % 2 ? -1.0f : 1.0f));
    }
  }

  return curLevel;
}

// friend function: overloading << and >>
ostream &operator<<(ostream &out, const Matrix &m)
{
  for (int r = 0; r < m.row; r++)
  {
    for (int c = 0; c < m.column; c++)
      out << setw(10) << fixed << setprecision(5) << m.data[r][c];
    if (r != m.row - 1)
      out << endl;
  }

  return out;
}
istream &operator>>(istream &in, Matrix &m)
{
  for (int r = 0; r < m.row; r++)
    for (int c = 0; c < m.column; c++)
      m.data[r][c] = getRandFloat();
  return in;
}

main.cpp

#include <iostream>
#include "Matrix.h"

int main()
{
  srand((unsigned)time(0)); // set random seed

  Matrix m(6, 6); //m is a 6x6 matrix with all value set to 0
  std::cin >> m; // give each element in the matrix a random float value, when using >> operator

  cout << "inverse of m:" << endl;
  cout << m.inverse() << endl;
  //"<<" operator is use to output all value int the matrix
 
  return 0;
}

The method I'm currently testing seems to be working fine. However, I want to confirm if it is the most efficient solution possible.

example

input

8.0867 3.15236 1.66545 1.26714 6.88726 4.24564
6.41356 2.64673 3.65965 7.65898 4.52873 4.3264
3.72865 7.36083 3.49968 9.13076 0.72390 6.62454
8.59957 2.96699 6.25450 9.29979 1.59475 2.93691
0.57955 0.50900 4.75473 2.79461 8.99372 7.43932
2.59571 6.15638 0.20169 9.80159 5.25058 6.5434

output

0.01317 0.64093 0.12824 -0.34249 -0.13077 -0.25974
0.24609 -1.50393 -0.20641 0.81144 0.13351 0.52768
0.07594 -0.94817 -0.13510 0.61556 0.19399 0.21758
-0.10456 0.13971 -0.06629 0.002489 -0.02853 0.07390
0.13616 -0.92044 -0.32327 0.55861 0.14560 0.43126
-0.19173 1.71926 0.50620 -1.09852 -0.15382 -0.70406
0

There are 0 best solutions below