fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3.  
  4. // 行列出力用関数，nは行列サイズ，Mは先頭ポインタ
  5. void printM(int n, double* M) {
  6. for (int i = 0; i < n; i++) {
  7. for (int j = 0; j < n; j++) {
  8. printf("%lf ", M[i * n + j]);
  9. }
  10. printf("\n");
  11. }
  12. }
  13. // ベクトル出力用関数，nはベクトルサイズ，Vは先頭ポインタ
  14. void printV(int n, double* V) {
  15. for (int j = 0; j < n; j++) {
  16. printf("%5.3f ", V[j]);
  17. }
  18. printf("\n");
  19. }
  20.  
  21. int main() {
  22. // 乱数セット
  23. srand(29);
  24. // 行列サイズ設定
  25. int n = 3;
  26. // A,x,bを格納するメモリの確保
  27. double* A = (double*)malloc(n * n * sizeof(double));
  28. double* x = (double*)malloc(n * sizeof(double));
  29. double* xAns = (double*)malloc(n * sizeof(double));
  30. double* b = (double*)malloc(n * sizeof(double));
  31.  
  32. // ベクトルの初期値を入れる
  33. for (int i = 0; i < n; i++) {
  34. x[i] = 0;
  35. b[i] = 0;
  36. xAns[i] = 5 + i; // 最終的に答えはx=5,6,7...となる
  37. }
  38.  
  39. // Aの行列を生成
  40. for (int i = 0; i < n; i++) {
  41. for (int j = 0; j < n; j++) {
  42. A[i * n + j] = 1 + (double)(rand() % 100) / 100; // 1~2の乱数で埋める
  43. }
  44. // 対角成分はその行の要素の総和にする(対角成分が大きい方が安定するため)
  45. for (int j = 0; j < n; j++) {
  46. double tmp = 0;
  47. tmp += A[i * n + j];
  48. A[i * n + i] = tmp;
  49. }
  50. }
  51.  
  52. // A,xAnsからbベクトルを逆算
  53. for (int i = 0; i < n; i++) {
  54. for (int j = 0; j < n; j++) {
  55. b[i] += A[i * n + j] * xAns[j];
  56. }
  57. }
  58.  
  59. /*------ここより下でAx=bを解く------*/
  60. 	double* L = (double*)malloc(n * n * sizeof(double));
  61. 	double* U = (double*)malloc(n * n * sizeof(double));
  62. 	double* L2 = (double*)malloc(n * n * sizeof(double));
  63. 	double* U2 = (double*)malloc(n * n * sizeof(double));
  64. 	double* Z = (double*)malloc(n * sizeof(double));
  65. 	for(int i = 0; i < n; i++){
  66. 		for(int j = 0; j < n; j++){
  67. 			L[i * n + j] = 0;
  68. 			U[i * n + j] = A[i * n + j];
  69. 			L2[i * n + j] = 0;
  70. 			U2[i * n + j] = 0;
  71. 		}
  72. 		Z[i] = b[i];
  73. 	}
  74.  
  75. 	for(int i = 0; i < n; i++){
  76. 		L[i * n + i] = U[i * n + i];
  77. 		for(int j = i; j < n; j++){
  78. 			U[i * n + j] /= L[i * n + i];
  79. 		}
  80. 		for(int j = i+1; j < n; j++){
  81. 			L[j * n + i] = U[j * n + i];
  82. 			for(int k = i; k < n; k++){
  83. 			U[j * n + k] -= U[i * n + k] * L[j * n + i];
  84. 			}
  85. 		}
  86. 	}
  87. 	for(int i = 0; i < n; i++){
  88. 		for(int j = 0; j < n; j++){
  89. 			L2[i * n + j] = L[i * n + j];
  90. 			U2[i * n + j] = U[i * n + j];
  91. 		}
  92. 	}
  93. 	for(int i = 0; i < n; i++){
  94. 		double ii = L2[i * n + i];
  95. 		L2[i * n + i] /= ii;
  96. 		Z[i] /= ii;
  97. 		for(int j = i+1; j < n; j++){
  98. 			double ji = L2[j * n + i];
  99. 			L2[j * n + i] -= L2[i * n + i] * ji;
  100. 			Z[j] -= Z[i] * ji;
  101. 		}
  102. 	}
  103. 	for(int i = 0; i < n; i++){
  104. 		x[i] = Z[i];
  105. 	}
  106. 	for(int i = n-1; i > 0; i--){
  107. 		for(int j = i-1; j >= 0; j--){
  108. 			double ji = U2[j * n + i];
  109. 			U2[j * n + i] -= U2[i * n + i] * ji;
  110. 			x[j] -= x[i] * ji;
  111. 		}
  112. 	}
  113.  
  114. /*------ここより上でAx=bを解く------*/
  115.  
  116. // 最後に残ったA,b,x,xの答えを出力
  117. printf("A:\n");
  118. printM(n, A);
  119. printf("b:\n");
  120. printV(n, b);
  121. printf("L:\n");
  122. printM(n, L);
  123. printf("U:\n");
  124. printM(n, U);
  125. printf("Z:\n");
  126. printV(n, Z);
  127. printf("x:\n");
  128. printV(n, x);
  129. printf("xAns:\n");
  130. printV(n, xAns);
  131. }
  132.  
Success #stdin #stdout 0s 5276KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
A:
1.040000 1.870000 1.040000 
1.450000 1.670000 1.670000 
1.950000 1.720000 1.720000 
b:
23.700 28.960 32.110 
L:
1.040000 0.000000 0.000000 
1.450000 -0.937212 0.000000 
1.950000 -1.786250 -0.649302 
U:
1.000000 1.798077 1.000000 
0.000000 1.000000 -0.234739 
0.000000 0.000000 1.000000 
Z:
22.788 4.357 7.000 
x:
5.000 6.000 7.000 
xAns:
5.000 6.000 7.000
https://ideone.com/iPhzlR
language:
C++ (gcc 8.3)
created:
1 year 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!