fork download
copy 
  1. #include <stdio.h>
  2.  
  3. #define N 4
  4.  
  5. int main() {
  6.  
  7. double A[N][N] = {
  8.     {9, 2, 1, 1},
  9.     {18, 4, 2, 2},
  10.     {1, -1, 3, 2},
  11.     {2, 1, 1, 5}
  12. };
  13.  
  14. double b[N] = {10, 25, 7, 8};
  15.  
  16. double L[N][N] = {0};
  17. double U[N][N];
  18.  
  19. double y[N] = {0};
  20. double x[N] = {0};
  21.  
  22. double LU[N][N] = {0};
  23. double Ax[N] = {0};
  24.  
  25. int i, j, k;
  26.  
  27. /* Copy A into U */
  28. for(i = 0; i < N; i++) {
  29. for(j = 0; j < N; j++) {
  30. U[i][j] = A[i][j];
  31. }
  32. }
  33.  
  34. /* Initialize L as identity matrix */
  35. for(i = 0; i < N; i++) {
  36. L[i][i] = 1.0;
  37. }
  38.  
  39. /* LU decomposition */
  40. for(i = 0; i < N - 1; i++) {
  41.  
  42. for(j = i + 1; j < N; j++) {
  43.  
  44. double factor = U[j][i] / U[i][i];
  45.  
  46. L[j][i] = factor;
  47.  
  48. for(k = 0; k < N; k++) {
  49. U[j][k] = U[j][k] - factor * U[i][k];
  50. }
  51. }
  52. }
  53.  
  54. /* Forward substitution: Ly = b */
  55. for(i = 0; i < N; i++) {
  56.  
  57. y[i] = b[i];
  58.  
  59. for(j = 0; j < i; j++) {
  60. y[i] -= L[i][j] * y[j];
  61. }
  62. }
  63.  
  64. /* Back substitution: Ux = y */
  65. for(i = N - 1; i >= 0; i--) {
  66.  
  67. x[i] = y[i];
  68.  
  69. for(j = i + 1; j < N; j++) {
  70. x[i] -= U[i][j] * x[j];
  71. }
  72.  
  73. x[i] /= U[i][i];
  74. }
  75.  
  76. /* Display L matrix */
  77. printf("L行列:\n");
  78.  
  79. for(i = 0; i < N; i++) {
  80. for(j = 0; j < N; j++) {
  81. printf("%8.3f ", L[i][j]);
  82. }
  83. printf("\n");
  84. }
  85.  
  86. /* Display U matrix */
  87. printf("\nU行列:\n");
  88.  
  89. for(i = 0; i < N; i++) {
  90. for(j = 0; j < N; j++) {
  91. printf("%8.3f ", U[i][j]);
  92. }
  93. printf("\n");
  94. }
  95.  
  96. /* Display y vector */
  97. printf("\nyベクトル:\n");
  98.  
  99. for(i = 0; i < N; i++) {
  100. printf("%8.3f\n", y[i]);
  101. }
  102.  
  103. /* Display solution x */
  104. printf("\n連立方程式の解 x:\n");
  105.  
  106. for(i = 0; i < N; i++) {
  107. printf("x%d = %8.3f\n", i + 1, x[i]);
  108. }
  109.  
  110. /* Check LU = A */
  111. for(i = 0; i < N; i++) {
  112. for(j = 0; j < N; j++) {
  113.  
  114. LU[i][j] = 0;
  115.  
  116. for(k = 0; k < N; k++) {
  117. LU[i][j] += L[i][k] * U[k][j];
  118. }
  119. }
  120. }
  121.  
  122. printf("\nLU行列（LU = A の確認）:\n");
  123.  
  124. for(i = 0; i < N; i++) {
  125. for(j = 0; j < N; j++) {
  126. printf("%8.3f ", LU[i][j]);
  127. }
  128. printf("\n");
  129. }
  130.  
  131. /* Check Ax = b */
  132. for(i = 0; i < N; i++) {
  133.  
  134. Ax[i] = 0;
  135.  
  136. for(j = 0; j < N; j++) {
  137. Ax[i] += A[i][j] * x[j];
  138. }
  139. }
  140.  
  141. printf("\nAxベクトル（Ax = b の確認）:\n");
  142.  
  143. for(i = 0; i < N; i++) {
  144. printf("%8.3f\n", Ax[i]);
  145. }
  146.  
  147. printf("\n元のbベクトル:\n");
  148.  
  149. for(i = 0; i < N; i++) {
  150. printf("%8.3f\n", b[i]);
  151. }
  152.  
  153. return 0;
  154. }
  155.  
Success #stdin #stdout 0s 5328KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
L行列:
   1.000    0.000    0.000    0.000 
   2.000    1.000    0.000    0.000 
   0.111     -inf    1.000    0.000 
   0.222      inf     -nan    1.000 

U行列:
   9.000    2.000    1.000    1.000 
   0.000    0.000    0.000    0.000 
    -nan     -nan     -nan     -nan 
    -nan     -nan     -nan     -nan 

yベクトル:
  10.000
   5.000
     inf
    -nan

連立方程式の解 x:
x1 =     -nan
x2 =     -nan
x3 =     -nan
x4 =     -nan

LU行列（LU = A の確認）:
    -nan     -nan     -nan     -nan 
    -nan     -nan     -nan     -nan 
    -nan     -nan     -nan     -nan 
    -nan     -nan     -nan     -nan 

Axベクトル（Ax = b の確認）:
    -nan
    -nan
    -nan
    -nan

元のbベクトル:
  10.000
  25.000
   7.000
   8.000
https://ideone.com/fLp5fI
language:
C (gcc 8.3)
created:
1 day 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!