C言語で逆行列を求める

Python なら numpy で簡単に求められる逆行列ですが、C ではそうもいきません。

そこで、先人たちの知恵を借りるべくいろいろとインターネットを漁ったのですが、自分にとってちょうどよいサンプルプログラムがなかったので、近しいサンプルを参考にして作り直してみました。

バグなどがあれば教えてください🙏

プログラム

inverse.c

#include <stdio.h>
#include <stdlib.h>

#define N 3

void inverse_matrix(long double **, int, long double **);
void get_minor(long double **, long double **, int, int, int);
long double calculate_determinant(long double **, int);

void inverse_matrix(long double **matrix, int order, long double **inv_matrix) {

  int i, j;
  long double inv_det;
  long double **minor;

  inv_det = 1.0 / calculate_determinant(matrix, order);
  
  minor = (long double **)calloc(order - 1, sizeof(long double *));
  for (i = 0; i < order - 1; i++)
    minor[i] = (long double *)calloc(order - 1, sizeof(long double));
  
  for (j = 0; j < order; j++) {
    for (i = 0; i < order; i++) {
      get_minor(matrix, minor, j, i, order);
      inv_matrix[i][j] = inv_det * calculate_determinant(minor, order - 1);
      if ( (i + j) % 2 == 1)
        inv_matrix[i][j] = -inv_matrix[i][j];
    }
  }
  
  for (i = 0; i < order - 1; i++)
    free(minor[i]);
  free(minor);
  
}


void get_minor(long double **src, long double **dest, int row, int col, int order) {
  
  int i, j;
  int col_count = 0;
  int row_count = 0;

  for (j = 0; j < order; j++) {
    if (j != row) {
      col_count = 0;
      for (i = 0; i < order; i++) {
        if (i != col) {
          dest[row_count][col_count] = src[j][i];
          col_count++;
        }
      }
      row_count++;
    }
  }

}


long double calculate_determinant(long double **matrix, int order) {

  int i;
  long double det = 0.0;
  long double **minor;
  
  if (order == 1)
    return matrix[0][0];
  
  minor = (long double **)calloc(order - 1, sizeof(long double *));
  for (i = 0; i < order - 1; i++)
    minor[i] = (long double *)calloc(order - 1, sizeof(long double));

  for (i = 0; i < order; i++) {
    get_minor(matrix, minor, 0, i, order);
    det += (i % 2 == 1 ? -1.0 : 1.0) * matrix[0][i] * calculate_determinant(minor, order - 1);
  }

  for (i = 0; i < order - 1; i++)
    free(minor[i]);
  free(minor);

  return det;
  
}


int main(void) {
  
  int i, j;
  long double **matrix, **inv_matrix;

  matrix = (long double **)calloc(N, sizeof(long double *));
  for (i = 0; i < N; i++)
    matrix[i] = (long double *)calloc(N, sizeof(long double));

  inv_matrix = (long double **)calloc(N, sizeof(long double *));
  for (i = 0; i < N; i++)
    inv_matrix[i] = (long double *)calloc(N, sizeof(long double));
  
  matrix[0][0] = 1;  matrix[0][1] = 1;  matrix[0][2] = 3;
  matrix[1][0] = 1;  matrix[1][1] = 3;  matrix[1][2] = -3;
  matrix[2][0] = -2; matrix[2][1] = -4; matrix[2][2] = -4;

  for (j = 0; j < N; j++) {
    for (i = 0; i < N; i++) {
      printf("%Lf, ", matrix[j][i]);
    }
    printf("\n");
  }

  /* 逆行列を求める */
  inverse_matrix(matrix, N, inv_matrix);
  
  for (j = 0; j < N; j++) {
    for (i = 0; i < N; i++) {
      printf("%Lf, ", inv_matrix[j][i]);
    }
    printf("\n");
  }

  for (i = 0; i < N; i++) {
    free(matrix[i]);
    free(inv_matrix[i]);
  }
  free(matrix);
  free(inv_matrix);
  
  return 0;
  
}

動作結果

1.000000, 1.000000, 3.000000,
1.000000, 3.000000, -3.000000,
-2.000000, -4.000000, -4.000000,

3.000000, 1.000000, 1.500000,
-1.250000, -0.250000, -0.750000,
-0.250000, -0.250000, -0.250000,