fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4
  5.  
  6. void pivot(double x[][N+1], int i) {
  7.     double c[N+1] = {0};
  8.     for (int j = i + 1; j < N; j++) {
  9.         if (fabs(x[i][i]) < fabs(x[j][i])) {
  10.             for (int n = 0; n < N + 1; n++) {
  11.                 c[n] = x[j][n];
  12.                 x[j][n] = x[i][n];
  13.                 x[i][n] = c[n];
  14.             }
  15.         }
  16.     }
  17. }
  18.  
  19. void lu(double x[][N+1], double l[][N], double u[][N], double y[N]) {
  20.     for (int i=0;i<N;i++) {
  21.         pivot(x, i);
  22.         l[i][i] = x[i][i];
  23.         for (int j=i;j<N;j++) {
  24.             u[i][j] = x[i][j];
  25.         }
  26.         for (int a=i+1;a<N;a++) {
  27.             l[a][i] = x[a][i];
  28.             double det = x[a][i] / x[i][i];
  29.             for (int b=i;b<N+1;b++) {
  30.                 x[a][b] -= det * x[i][b];
  31.             }
  32.         }
  33.     }
  34.     for (int i=0;i<N;i++) {
  35.         y[i] = x[i][N];
  36.     }
  37. }
  38.  
  39. void sol(double x[N], double u[][N], double y[N]) {
  40.     for (int i=N-1;i>=0;i--) {
  41.         x[i] = y[i];
  42.         for (int j=i+1;j<N;j++) {
  43.             x[i] -= u[i][j] * x[j];
  44.         }
  45.         x[i] /= u[i][i];
  46.     }
  47.     printf("x:\n");
  48.     for (int i=0;i<N;i++) {
  49.         printf("x%d = %f\n", i+1, x[i]);
  50.     }
  51.     printf("\n");
  52. }
  53.  
  54. void con(double a[][N], double x[N]){
  55. double result[N] = {0};
  56.     for (int i = 0; i < N; i++) {
  57.         for (int j = 0; j < N; j++) {
  58.             result[i] += a[i][j] * x[j];
  59.         }
  60.     }
  61.     for (int i = 0; i < N; i++) { 
  62.     	printf("%f\n", result[i]); 
  63.     } 
  64. 	printf("\n");
  65. }
  66.  
  67. int main(void) {
  68.     double x[][N+1] = {
  69.         {3, 1.5, -6, 4.8, 1.2},
  70.         {1, 1.5, -2, -2.4, 0.6},
  71.         {0, -1.5, -2, -1, -2.4},
  72.         {2, 4, -1.8, -0.6, 0}
  73.     };
  74.  
  75.     double l[N][N] = {0};
  76.     double u[N][N] = {0};
  77.     double y[N] = {0};
  78.     double s[N] = {0};
  79.  
  80.     double a[4][4];
  81. 	for(int i=0;i<4;i++){
  82. 		for(int j=0;j<4;j++){
  83. 			a[i][j]=x[i][j];
  84. 		}
  85. 	}
  86.  
  87.     lu(x, l, u, y);
  88.     printf("U行列:\n");
  89.     for (int i=0;i<N;i++) {
  90.         for (int j=0;j<N;j++) {
  91.             printf("%f ", u[i][j]);
  92.         }
  93.         printf("\n");
  94.     }
  95.     printf("\n");
  96.  
  97.     printf("L行列:\n");
  98.     for (int i=0;i<N;i++) {
  99.         for (int j=0;j<N;j++) {
  100.             printf("%f ", l[i][j]);
  101.         }
  102.         printf("\n");
  103.     }
  104.     printf("\n");
  105.  
  106.     printf("y:\n");
  107.     for (int i=0;i<N;i++){
  108.         printf("y%d = %f\n", i+1, y[i]);
  109.     }
  110.     printf("\n");
  111.  
  112.     sol(s, u, y);
  113.     con(a, s);
  114.     return 0;
  115. }
Success #stdin #stdout 0.01s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
U行列:
3.000000 1.500000 -6.000000 4.800000 
0.000000 3.000000 2.200000 -3.800000 
0.000000 0.000000 -0.900000 -2.900000 
0.000000 0.000000 0.000000 -0.370370 

L行列:
3.000000 0.000000 0.000000 0.000000 
1.000000 3.000000 0.000000 0.000000 
0.000000 1.000000 -0.900000 0.000000 
2.000000 -1.500000 -0.733333 -0.370370 

y:
y1 = 1.200000
y2 = -0.800000
y3 = -2.800000
y4 = 2.748148

x:
x1 = 81.052000
x2 = -29.480000
x3 = 27.020000
x4 = -7.420000

1.200000
0.600000
-2.400000
0.000000
https://ideone.com/zhmu55
language:
C (gcc 8.3)
created:
2 days 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!