#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double *cholesky(double *A, int n) {
double *L = (double*)calloc(n * n, sizeof(double));
if (L == NULL)
exit(EXIT_FAILURE);
for (int i = 0; i < n; i++)
for (int j = 0; j < (i + 1); j++) {
double s = 0;
for (int k = 0; k < j; k++)
s += L[i * n + k] * L[j * n + k];
L[i * n + j] = (i == j) ?
sqrt(A[i * n + i] - s) :
(1.0 / L[j * n + j] * (A[i * n + j] - s));
}
return L;
}
int res;
void show_matrix(double *A, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
printf("%lf ", A[i*n + j]);
printf("\n");
}
}
void changeMatr(double *A, int n, char* file_name) {
FILE *input = fopen(file_name, "r");
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
fscanf(input, "%f", A+i*n+j);
}
fclose(input);
}
int main() {
int n = 10;
double* m = new double[n*n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
m[i*n + j] = 1.0 / (i + j + 1.0);
}
}
double* b = new double[n];
for (int i = 0; i < n; i++) {
b[i] = 1;
}
show_matrix(m, n);
double *res = cholesky(m, n);
printf("---------\n");
show_matrix(res, n);
double* mT = new double[n*n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
mT[i*n + j] = res[j*n + i];
}
}
printf("---------\n");
show_matrix(mT, n);
printf("=========\n");
changeMatr(m, 10, "input1.txt");
show_matrix(m, n);
printf("=========\n");
double* y = new double[n];
double* x = new double[n];
for (int i = 0; i < n; i++) {
y[i] = b[i];
for (int j = 0; j < i; j++) {
y[i] -= res[i*n + j] * y[j];
}
y[i] /= res[i*n + i];
}
for (int i = n - 1; i >= 0; i--) {
x[i] = y[i];
for (int j = n - 1; j > i; j--) {
x[i] -= mT[i*n + j] * x[j];
}
x[i] /= mT[i*n + i];
}
printf("---------\n");
for (int i = 0; i < n; i++) {
printf("[ ");
for (int j = 0; j < n; j++)
printf("%lf ", res[i*n + j]);
printf("] [%lf] = [%lf]\n", y[i], b[i]);
}
printf("---------\n");
for (int i = 0; i < n; i++) {
printf("[ ");
for (int j = 0; j < n; j++)
printf("%lf ", mT[i*n + j]);
printf("] [%lf] = [%lf]\n", x[i], y[i]);
}
printf("---------\n");
for (int i = 0; i < n; i++) {
printf("[ ");
for (int j = 0; j < n; j++)
printf("%lf ", m[i*n + j]);
printf("] [%lf] = [%lf]\n", x[i], b[i]);
}
free(res);
free(m);
free(mT);
free(b);
free(x);
free(y);
//getchar();
return 0;
}