C++实现n阶矩阵求逆：伴随矩阵法的高效优化方案问询

我通过伴随矩阵和行列式的方法，用C++实现了n阶矩阵的求逆操作，其中行列式采用递归方式计算。当前代码测试可正常运行，但受作业限制无法使用第三方库，想确认是否存在更优的实现方式。

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;
}

示例输入 #

示例输出 #

内容的提问来源于stack exchange，提问作者陳逸軒