0

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
陳逸軒
  • 11
  • 2
  • 1
    No you still do too much dynamic memory allocation. Make your matrix a template class with fixed sizes. Having said that, if you really want to write your own matrix class, the only thing you really should be doing is to profile your code by running a lot of matrix operations and finding your bottlenecks from there (an they can be platform specific : depend on cache, branch prediction, vectorization by compiler etc.). All in all if you really need fast matrix operations, use a library – Pepijn Kramer Jun 08 '23 at 05:30
  • 1
    Note recursion is usually not the best at producing fastest code possible either, change to an iterative algorithm (and yes every recursive algorithm can be transformed into an iterative on, or the other way around). – Pepijn Kramer Jun 08 '23 at 05:31
  • Due to the constraints of the assignment, I cannot use any other libraries. I will explore alternative methods.Thank you for your suggestions! – 陳逸軒 Jun 08 '23 at 05:42
  • 1
    The assignment, and I cannot use an external library" part is a big constraint that you should mention in your question. Another optimization point, avoid if statements within your loops if you can (doing a bit more calculations can even speed up things), missed branch predictions will slow down your execution too. – Pepijn Kramer Jun 08 '23 at 05:47
  • You could ask if there is a worse way, and the answer would be "probably not". – n. m. could be an AI Jun 08 '23 at 05:59
  • 1
    Your `Matrix` class violates the rule of 3 due to it not having a user-defined assignment operator. Thus anything to do with the copy-semantics of `Matrix`, such as passing and returning `Matrix` by value, is suspect. `{Matrix m(6, 6); Matrix m2(5,5); m = m2;}` -- If your destructor attempts to `delete[]` the memory for `Matrix`, you now have a double-delete error. – PaulMcKenzie Jun 08 '23 at 06:30
  • Then you have this: `Matrix *result = new Matrix(n, n);`. There is no need to create a `Matrix` dynamically here. All that should be done is `Matrix result(n,n);` and simply `return result;`. *Due to the constraints of the assignment* -- Whoever assigned this to you either didn't explain in detail what is required to write such a class so that it has correct copy semantics, or you were taught this but are leaving out this important part of your `Matrix` class. – PaulMcKenzie Jun 08 '23 at 06:35
  • 1
    determinant based inverse is fast only for small matrices for bigger ones use GEM (Gauss elimination method) its much faster however its a bit tricky to make it robust due to precision issues ...however it can be done inplace without any recursions see [Gaussian elimination without result for acceleration](https://stackoverflow.com/a/55439139/2521214) for some inspiration both GEM and subdeterminant versions are there even template version – Spektre Jun 08 '23 at 07:44
  • 1
    Agree with @Spektre - the adjugate / determinant method is good for proving theorems, but terrible for computation. The Gauss-elimination-type method (applying the same row operations to the identity matrix as you do to your original matrix whilst putting it in row-echelon form) is faster and can (within precision constraints) also give you the rank of the matrix if it turns out to be singular. – lastchance Jun 08 '23 at 08:09
  • LU decomposition is the way to go. No large scale linear algebra package solves matrix equations this way. – duffymo Jun 08 '23 at 09:47

0 Answers0