fork(2) download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3.  
  4. #define N1 3
  5. #define N2 4
  6.  
  7. void make_square_matrix(double ***a, int n)
  8. {
  9.   int i;
  10.  
  11.   if ((*a = (double **)malloc(n * sizeof(double *))) == NULL)
  12.     exit(1);
  13.   for (i = 0; i < n; i++)
  14.     if (((*a)[i] = (double *)malloc(n * sizeof(double))) == NULL)
  15.       exit(1);
  16. }
  17.  
  18. void free_square_matrix(double **a, int n)
  19. {
  20.   int i;
  21.  
  22.   for (i = 0; i < n; i++)
  23.     free(a[i]);
  24.   free(a);
  25. }
  26.  
  27. /* LU分解(ドゥーリトル法) */
  28. void LU_decomp(double **a, double **l, double **u, int n)
  29. {
  30.   int i, j, k;
  31.  
  32.   for (i = 0; i < n; i++)
  33.     l[i][i] = 1.0;
  34.   for (i = 0; i < n; i++)
  35.     u[0][i] = a[0][i];
  36.  
  37.   for (i = 1; i < n; i++) {
  38.     if (u[i - 1][i - 1] == 0.0)
  39.       exit(1);
  40.     for (j = i; j < n; j++) {
  41.       // L成分の算出
  42.       l[j][i - 1] = a[j][i - 1];
  43.       for (k = 0; k < i - 1; k++)
  44.         l[j][i - 1] -= l[j][k] * u[k][i - 1];
  45.       l[j][i - 1] /= u[i - 1][i - 1];
  46.       // U成分の算出
  47.       u[i][j] = a[i][j];
  48.       for (k = 0; k < i; k++)
  49.         u[i][j] -= l[i][k] * u[k][j];
  50.     }
  51.   }
  52. }
  53.  
  54. void print_matrix(double **a, int n, char *s)
  55. {
  56.   int i, j;
  57.  
  58.   printf("*** %s ***\n", s);
  59.  
  60.   for (i = 0; i < n; i++) {
  61.     for (j = 0; j < n; j++)
  62.       printf("% 11.7f  ", a[i][j]);
  63.     putchar('\n');
  64.   }
  65. }
  66.  
  67. /* 三角行列の行列式は対角成分の積 */
  68. double det_triangular_matrix(double **a, int n)
  69. {
  70.   int i;
  71.   double det = 1.0;
  72.  
  73.   for (i = 0; i < n; i++)
  74.     det *= a[i][i];
  75.  
  76.   return det;
  77. }
  78.  
  79. int main(void)
  80. {
  81.   double **a, **l, **u, det;
  82.   int i, j;
  83.  
  84.   double a1[N1][N1] = {{2.0, 3.0, -1.0}, {4.0, 4.0, -3.0}, {2.0, -3.0, 1.0}};
  85.   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}};
  86.  
  87.   make_square_matrix(&a, N1);
  88.   make_square_matrix(&l, N1);
  89.   make_square_matrix(&u, N1);
  90.  
  91.   for (i = 0; i < N1; i++)
  92.     for (j = 0; j < N1; j++)
  93.       a[i][j] = a1[i][j];
  94.  
  95.   LU_decomp(a, l, u, N1);
  96.  
  97.   print_matrix(a, N1, "行列 A");
  98.   print_matrix(l, N1, "行列 L");
  99.   print_matrix(u, N1, "行列 U");
  100.  
  101.   /* |A| = |L||U| */
  102.   det = det_triangular_matrix(l, N1) * det_triangular_matrix(u, N1);
  103.   printf("行列式 |A| = %f\n", det);
  104.  
  105.   free_square_matrix(a, N1);
  106.   free_square_matrix(l, N1);
  107.   free_square_matrix(u, N1);
  108.  
  109.   putchar('\n');
  110.  
  111.   make_square_matrix(&a, N2);
  112.   make_square_matrix(&l, N2);
  113.   make_square_matrix(&u, N2);
  114.  
  115.   for (i = 0; i < N2; i++)
  116.     for (j = 0; j < N2; j++)
  117.       a[i][j] = a2[i][j];
  118.  
  119.   LU_decomp(a, l, u, N2);
  120.  
  121.   print_matrix(a, N2, "行列 A");
  122.   print_matrix(l, N2, "行列 L");
  123.   print_matrix(u, N2, "行列 U");
  124.  
  125.   det = det_triangular_matrix(l, N2) * det_triangular_matrix(u, N2);
  126.   printf("行列式 |A| = %f\n", det);
  127.  
  128.   free_square_matrix(a, N2);
  129.   free_square_matrix(l, N2);
  130.   free_square_matrix(u, N2);
  131.  
  132.   return 0;
  133. }
Success #stdin #stdout 0.01s 1852KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
*** 行列 A ***
  2.0000000    3.0000000   -1.0000000  
  4.0000000    4.0000000   -3.0000000  
  2.0000000   -3.0000000    1.0000000  
*** 行列 L ***
  1.0000000    0.0000000    0.0000000  
  2.0000000    1.0000000    0.0000000  
  1.0000000    3.0000000    1.0000000  
*** 行列 U ***
  2.0000000    3.0000000   -1.0000000  
  0.0000000   -2.0000000   -1.0000000  
  0.0000000    0.0000000    5.0000000  
行列式 |A| = -20.000000

*** 行列 A ***
  1.0000000    1.0000000    0.0000000    3.0000000  
  2.0000000    1.0000000   -1.0000000    1.0000000  
  3.0000000   -1.0000000   -1.0000000    2.0000000  
 -1.0000000    2.0000000    3.0000000   -2.0000000  
*** 行列 L ***
  1.0000000    0.0000000    0.0000000    0.0000000  
  2.0000000    1.0000000    0.0000000    0.0000000  
  3.0000000    4.0000000    1.0000000    0.0000000  
 -1.0000000   -3.0000000    0.0000000    1.0000000  
*** 行列 U ***
  1.0000000    1.0000000    0.0000000    3.0000000  
  0.0000000   -1.0000000   -1.0000000   -5.0000000  
  0.0000000    0.0000000    3.0000000   13.0000000  
  0.0000000    0.0000000    0.0000000  -14.0000000  
行列式 |A| = 42.000000
https://ideone.com/Duslx
language:
C (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!