#include <stdio.h>
#include <stdlib.h>
// 行列出力用関数,nは行列サイズ,Mは先頭ポインタ
void printM(int n, double* M) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
printf("%lf ", M[i * n + j]);
}
printf("\n");
}
}
// ベクトル出力用関数,nはベクトルサイズ,Vは先頭ポインタ
void printV(int n, double* V) {
for (int j = 0; j < n; j++) {
printf("%5.3f ", V[j]);
}
printf("\n");
}
int main() {
// 乱数セット
srand(29);
// 行列サイズ設定
int n = 3;
// A,x,bを格納するメモリの確保
double* A = (double*)malloc(n * n * sizeof(double));
double* x = (double*)malloc(n * sizeof(double));
double* xAns = (double*)malloc(n * sizeof(double));
double* b = (double*)malloc(n * sizeof(double));
// ベクトルの初期値を入れる
for (int i = 0; i < n; i++) {
x[i] = 0;
b[i] = 0;
xAns[i] = 5 + i; // 最終的に答えはx=5,6,7...となる
}
// Aの行列を生成
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i * n + j] = 1 + (double)(rand() % 100) / 100; // 1~2の乱数で埋める
}
// 対角成分はその行の要素の総和にする(対角成分が大きい方が安定するため)
for (int j = 0; j < n; j++) {
double tmp = 0;
tmp += A[i * n + j];
A[i * n + i] = tmp;
}
}
// A,xAnsからbベクトルを逆算
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
b[i] += A[i * n + j] * xAns[j];
}
}
/*------ここより下でAx=bを解く------*/
double* L = (double*)malloc(n * n * sizeof(double));
double* U = (double*)malloc(n * n * sizeof(double));
double* L2 = (double*)malloc(n * n * sizeof(double));
double* U2 = (double*)malloc(n * n * sizeof(double));
double* Z = (double*)malloc(n * sizeof(double));
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
L[i * n + j] = 0;
U[i * n + j] = A[i * n + j];
L2[i * n + j] = 0;
U2[i * n + j] = 0;
}
Z[i] = b[i];
}
for(int i = 0; i < n; i++){
L[i * n + i] = U[i * n + i];
for(int j = i; j < n; j++){
U[i * n + j] /= L[i * n + i];
}
for(int j = i+1; j < n; j++){
L[j * n + i] = U[j * n + i];
for(int k = i; k < n; k++){
U[j * n + k] -= U[i * n + k] * L[j * n + i];
}
}
}
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
L2[i * n + j] = L[i * n + j];
U2[i * n + j] = U[i * n + j];
}
}
for(int i = 0; i < n; i++){
double ii = L2[i * n + i];
L2[i * n + i] /= ii;
Z[i] /= ii;
for(int j = i+1; j < n; j++){
double ji = L2[j * n + i];
L2[j * n + i] -= L2[i * n + i] * ji;
Z[j] -= Z[i] * ji;
}
}
for(int i = 0; i < n; i++){
x[i] = Z[i];
}
for(int i = n-1; i > 0; i--){
for(int j = i-1; j >= 0; j--){
double ji = U2[j * n + i];
U2[j * n + i] -= U2[i * n + i] * ji;
x[j] -= x[i] * ji;
}
}
/*------ここより上でAx=bを解く------*/
// 最後に残ったA,b,x,xの答えを出力
printf("A:\n");
printM(n, A);
printf("b:\n");
printV(n, b);
printf("L:\n");
printM(n, L);
printf("U:\n");
printM(n, U);
printf("Z:\n");
printV(n, Z);
printf("x:\n");
printV(n, x);
printf("xAns:\n");
printV(n, xAns);
}