#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);
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KIAovLyDooYzliJflh7rlipvnlKjplqLmlbDvvIxu44Gv6KGM5YiX44K144Kk44K677yMTeOBr+WFiOmgreODneOCpOODs+OCvwp2b2lkIHByaW50TShpbnQgbiwgZG91YmxlKiBNKSB7CmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CnByaW50ZigiJWxmICIsIE1baSAqIG4gKyBqXSk7Cn0KcHJpbnRmKCJcbiIpOwp9Cn0KLy8g44OZ44Kv44OI44Or5Ye65Yqb55So6Zai5pWw77yMbuOBr+ODmeOCr+ODiOODq+OCteOCpOOCuu+8jFbjga/lhYjpoK3jg53jgqTjg7Pjgr8Kdm9pZCBwcmludFYoaW50IG4sIGRvdWJsZSogVikgewpmb3IgKGludCBqID0gMDsgaiA8IG47IGorKykgewpwcmludGYoIiU1LjNmICIsIFZbal0pOwp9CnByaW50ZigiXG4iKTsKfQogCmludCBtYWluKCkgewovLyDkubHmlbDjgrvjg4Pjg4gKc3JhbmQoMjkpOwovLyDooYzliJfjgrXjgqTjgrroqK3lrpoKaW50IG4gPSAzOwovLyBBLHgsYuOCkuagvOe0jeOBmeOCi+ODoeODouODquOBrueiuuS/nQpkb3VibGUqIEEgPSAoZG91YmxlKiltYWxsb2MobiAqIG4gKiBzaXplb2YoZG91YmxlKSk7CmRvdWJsZSogeCA9IChkb3VibGUqKW1hbGxvYyhuICogc2l6ZW9mKGRvdWJsZSkpOwpkb3VibGUqIHhBbnMgPSAoZG91YmxlKiltYWxsb2MobiAqIHNpemVvZihkb3VibGUpKTsKZG91YmxlKiBiID0gKGRvdWJsZSopbWFsbG9jKG4gKiBzaXplb2YoZG91YmxlKSk7CiAKLy8g44OZ44Kv44OI44Or44Gu5Yid5pyf5YCk44KS5YWl44KM44KLCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CnhbaV0gPSAwOwpiW2ldID0gMDsKeEFuc1tpXSA9IDUgKyBpOyAvLyDmnIDntYLnmoTjgavnrZTjgYjjga94PTUsNiw3Li4u44Go44Gq44KLCn0KIAovLyBB44Gu6KGM5YiX44KS55Sf5oiQCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CkFbaSAqIG4gKyBqXSA9IDEgKyAoZG91YmxlKShyYW5kKCkgJSAxMDApIC8gMTAwOyAvLyAxfjLjga7kubHmlbDjgafln4vjgoHjgosKfQovLyDlr77op5LmiJDliIbjga/jgZ3jga7ooYzjga7opoHntKDjga7nt4/lkozjgavjgZnjgoso5a++6KeS5oiQ5YiG44GM5aSn44GN44GE5pa544GM5a6J5a6a44GZ44KL44Gf44KBKQpmb3IgKGludCBqID0gMDsgaiA8IG47IGorKykgewpkb3VibGUgdG1wID0gMDsKdG1wICs9IEFbaSAqIG4gKyBqXTsKQVtpICogbiArIGldID0gdG1wOwp9Cn0KIAovLyBBLHhBbnPjgYvjgoli44OZ44Kv44OI44Or44KS6YCG566XCmZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKSB7CmZvciAoaW50IGogPSAwOyBqIDwgbjsgaisrKSB7CmJbaV0gKz0gQVtpICogbiArIGpdICogeEFuc1tqXTsKfQp9CiAKLyotLS0tLS3jgZPjgZPjgojjgorkuIvjgadBeD1i44KS6Kej44GPLS0tLS0tKi8KCWRvdWJsZSogTCA9IChkb3VibGUqKW1hbGxvYyhuICogbiAqIHNpemVvZihkb3VibGUpKTsKCWRvdWJsZSogVSA9IChkb3VibGUqKW1hbGxvYyhuICogbiAqIHNpemVvZihkb3VibGUpKTsKCWRvdWJsZSogTDIgPSAoZG91YmxlKiltYWxsb2MobiAqIG4gKiBzaXplb2YoZG91YmxlKSk7Cglkb3VibGUqIFUyID0gKGRvdWJsZSopbWFsbG9jKG4gKiBuICogc2l6ZW9mKGRvdWJsZSkpOwoJZG91YmxlKiBaID0gKGRvdWJsZSopbWFsbG9jKG4gKiBzaXplb2YoZG91YmxlKSk7Cglmb3IoaW50IGkgPSAwOyBpIDwgbjsgaSsrKXsKCQlmb3IoaW50IGogPSAwOyBqIDwgbjsgaisrKXsKCQkJTFtpICogbiArIGpdID0gMDsKCQkJVVtpICogbiArIGpdID0gQVtpICogbiArIGpdOwoJCQlMMltpICogbiArIGpdID0gMDsKCQkJVTJbaSAqIG4gKyBqXSA9IDA7CgkJfQoJCVpbaV0gPSBiW2ldOwoJfQoKCWZvcihpbnQgaSA9IDA7IGkgPCBuOyBpKyspewoJCUxbaSAqIG4gKyBpXSA9IFVbaSAqIG4gKyBpXTsKCQlmb3IoaW50IGogPSBpOyBqIDwgbjsgaisrKXsKCQkJVVtpICogbiArIGpdIC89IExbaSAqIG4gKyBpXTsKCQl9CgkJZm9yKGludCBqID0gaSsxOyBqIDwgbjsgaisrKXsKCQkJTFtqICogbiArIGldID0gVVtqICogbiArIGldOwoJCQlmb3IoaW50IGsgPSBpOyBrIDwgbjsgaysrKXsKCQkJVVtqICogbiArIGtdIC09IFVbaSAqIG4gKyBrXSAqIExbaiAqIG4gKyBpXTsKCQkJfQoJCX0KCX0KCWZvcihpbnQgaSA9IDA7IGkgPCBuOyBpKyspewoJCWZvcihpbnQgaiA9IDA7IGogPCBuOyBqKyspewoJCQlMMltpICogbiArIGpdID0gTFtpICogbiArIGpdOwoJCQlVMltpICogbiArIGpdID0gVVtpICogbiArIGpdOwoJCX0KCX0KCWZvcihpbnQgaSA9IDA7IGkgPCBuOyBpKyspewoJCWRvdWJsZSBpaSA9IEwyW2kgKiBuICsgaV07CgkJTDJbaSAqIG4gKyBpXSAvPSBpaTsKCQlaW2ldIC89IGlpOwoJCWZvcihpbnQgaiA9IGkrMTsgaiA8IG47IGorKyl7CgkJCWRvdWJsZSBqaSA9IEwyW2ogKiBuICsgaV07CgkJCUwyW2ogKiBuICsgaV0gLT0gTDJbaSAqIG4gKyBpXSAqIGppOwoJCQlaW2pdIC09IFpbaV0gKiBqaTsKCQl9Cgl9Cglmb3IoaW50IGkgPSAwOyBpIDwgbjsgaSsrKXsKCQl4W2ldID0gWltpXTsKCX0KCWZvcihpbnQgaSA9IG4tMTsgaSA+IDA7IGktLSl7CgkJZm9yKGludCBqID0gaS0xOyBqID49IDA7IGotLSl7CgkJCWRvdWJsZSBqaSA9IFUyW2ogKiBuICsgaV07CgkJCVUyW2ogKiBuICsgaV0gLT0gVTJbaSAqIG4gKyBpXSAqIGppOwoJCQl4W2pdIC09IHhbaV0gKiBqaTsKCQl9Cgl9CgovKi0tLS0tLeOBk+OBk+OCiOOCiuS4iuOBp0F4PWLjgpLop6PjgY8tLS0tLS0qLwogCi8vIOacgOW+jOOBq+aui+OBo+OBn0EsYix4LHjjga7nrZTjgYjjgpLlh7rlipsKcHJpbnRmKCJBOlxuIik7CnByaW50TShuLCBBKTsKcHJpbnRmKCJiOlxuIik7CnByaW50VihuLCBiKTsKcHJpbnRmKCJMOlxuIik7CnByaW50TShuLCBMKTsKcHJpbnRmKCJVOlxuIik7CnByaW50TShuLCBVKTsKcHJpbnRmKCJaOlxuIik7CnByaW50VihuLCBaKTsKcHJpbnRmKCJ4OlxuIik7CnByaW50VihuLCB4KTsKcHJpbnRmKCJ4QW5zOlxuIik7CnByaW50VihuLCB4QW5zKTsKfQo=