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

#define N1 3
#define N2 4

void make_square_matrix(double ***a, int n)
{
  int i;
  
  if ((*a = (double **)malloc(n * sizeof(double *))) == NULL)
    exit(1);
  for (i = 0; i < n; i++)
    if (((*a)[i] = (double *)malloc(n * sizeof(double))) == NULL)
      exit(1);
}

void free_square_matrix(double **a, int n)
{
  int i;
  
  for (i = 0; i < n; i++)
    free(a[i]);
  free(a);
}

/* LU分解(ドゥーリトル法) */
void LU_decomp(double **a, double **l, double **u, int n)
{
  int i, j, k;

  for (i = 0; i < n; i++)
    l[i][i] = 1.0;
  for (i = 0; i < n; i++)
    u[0][i] = a[0][i];

  for (i = 1; i < n; i++) {
    if (u[i - 1][i - 1] == 0.0)
      exit(1);
    for (j = i; j < n; j++) {
      // L成分の算出
      l[j][i - 1] = a[j][i - 1];
      for (k = 0; k < i - 1; k++)
        l[j][i - 1] -= l[j][k] * u[k][i - 1];
      l[j][i - 1] /= u[i - 1][i - 1];
      // U成分の算出
      u[i][j] = a[i][j];
      for (k = 0; k < i; k++)
        u[i][j] -= l[i][k] * u[k][j];
    }
  }
}

void print_matrix(double **a, int n, char *s)
{
  int i, j;
  
  printf("*** %s ***\n", s);
  
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
      printf("% 11.7f  ", a[i][j]);
    putchar('\n');
  }
}

/* 三角行列の行列式は対角成分の積 */
double det_triangular_matrix(double **a, int n)
{
  int i;
  double det = 1.0;
  
  for (i = 0; i < n; i++)
    det *= a[i][i];
  
  return det;
}

int main(void)
{
  double **a, **l, **u, det;
  int i, j;
  
  double a1[N1][N1] = {{2.0, 3.0, -1.0}, {4.0, 4.0, -3.0}, {2.0, -3.0, 1.0}};
  double a2[N2][N2] = {{1.0, 1.0, 0.0, 3.0}, {2.0, 1.0, -1.0, 1.0}, {3.0, -1.0, -1.0, 2.0}, {-1.0, 2.0, 3.0, -2.0}};
  
  make_square_matrix(&a, N1);
  make_square_matrix(&l, N1);
  make_square_matrix(&u, N1);
  
  for (i = 0; i < N1; i++)
    for (j = 0; j < N1; j++)
      a[i][j] = a1[i][j];
  
  LU_decomp(a, l, u, N1);
  
  print_matrix(a, N1, "行列 A");
  print_matrix(l, N1, "行列 L");
  print_matrix(u, N1, "行列 U");

  /* |A| = |L||U| */
  det = det_triangular_matrix(l, N1) * det_triangular_matrix(u, N1);
  printf("行列式 |A| = %f\n", det);
  
  free_square_matrix(a, N1);
  free_square_matrix(l, N1);
  free_square_matrix(u, N1);

  putchar('\n');
  
  make_square_matrix(&a, N2);
  make_square_matrix(&l, N2);
  make_square_matrix(&u, N2);
  
  for (i = 0; i < N2; i++)
    for (j = 0; j < N2; j++)
      a[i][j] = a2[i][j];
  
  LU_decomp(a, l, u, N2);
  
  print_matrix(a, N2, "行列 A");
  print_matrix(l, N2, "行列 L");
  print_matrix(u, N2, "行列 U");
  
  det = det_triangular_matrix(l, N2) * det_triangular_matrix(u, N2);
  printf("行列式 |A| = %f\n", det);

  free_square_matrix(a, N2);
  free_square_matrix(l, N2);
  free_square_matrix(u, N2);

  return 0;
}