#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#define M 20
#define Q 1.001 - 2 * M * 0.001
void scan_mat(FILE *file, double **mat, int size);
void print_mat(FILE *output, double **mat, int size, int flag);
void cpy_mat(double **dst_mat, double **src_mat, int size);
void print_ans(FILE *output, double **mat, int size, int *ord);
double **create_mat(int size);
void free_mat(double **mat, int size);
void swap(void *ptr_1, void *ptr_2, int size_of_cells);
double norm_diff_vec(double *vec_1, double *vec_2, int size);
void scan_mat(FILE *input, double **a, int n) {
int i, j;
for (i = 0; i < n; i++)
for (j = 0; j <= n; j++)
fscanf(input, "%lf", &a[i][j]);
}
void print_mat(FILE *output, double **a, int n, int flag) {
int i, j;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
fprintf(output, "%8.3lf ", a[i][j]);
if (flag)
fprintf(output, "| %8.3lf", a[i][n]);
fprintf(output, "\n");
}
fprintf(output, "\n");
}
void cpy_mat(double **m, double **a, int n) {
int i;
for (i = 0; i < n; i++)
memcpy(m[i], a[i], (n + 1) * sizeof(**m));
}
void print_ans(FILE *output, double **a, int n, int *ord) {
int i, j;
for (i = 0; i < n; i++) {
for (j = 0; ord[j] != i; j++)
;
fprintf(output, "x%d = %lf\n", i + 1, a[j][n]);
}
fprintf(output, "\n");
}
double **create_mat(int n) {
double **a = malloc(n * sizeof(*a));
int i;
for (i = 0; i < n; i++)
a[i] = calloc(n + 1, sizeof(**a));
return a;
}
void free_mat(double **a, int n) {
int i;
for (i = 0; i < n; i++)
free(a[i]);
free(a);
}
void swap(void *a, void *b, int size) {
void *tmp = malloc(size);
memcpy(tmp, a, size);
memcpy(a, b, size);
memcpy(b, tmp, size);
free(tmp);
}
double norm_diff_vec(double *a, double *b, int n) {
int i;
double tmp = 0;
for (i = 0; i < n; i++) {
tmp += (a[i] - b[i]) * (a[i] - b[i]);
}
return sqrt(tmp);
}
int main(void) {
FILE *output = fopen("out.txt", "w");
int n, i, j, k;
double **bcp, **m, **rev, det = 1, w, e, tmp;
printf("Choose way entering matrix:\n1. From file\n2. According to formulae\n\n");
char c;
do {
c = getch();
} while (c != '1' && c != '2');
if (c == '1') {
FILE *input = fopen("in.txt", "r");
if (!input) {
printf("File 'in.txt' not found!\n");
exit(1);
}
fscanf(input, "%d", &n);
bcp = create_mat(n);
scan_mat(input, bcp, n);
printf("Matrix readed from file.\n\n");
close(input);
} else {
n = 40;
bcp = create_mat(n);
printf("Enter x:\n");
scanf("%lf", &tmp);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (i == j) {
bcp[i][j] = pow(Q - 1, i + j);
} else {
bcp[i][j] = pow(Q, i + j) + 0.1 * (j - i);
}
}
}
for (i = 0; i < n; i++) {
bcp[i][n] = fabs(tmp - n / 10) * i * sin(tmp);
}
printf("Matrix filled.\n\n");
}
fprintf(output, "Source matrix:\n\n");
print_mat(output, bcp, n, 1);
printf("Enter w:\n");
scanf("%lf", &w);
printf("\nEnter e:\n");
scanf("%lf", &e);
m = create_mat(n);
rev = create_mat(n);
int ord[n + 1];
for (i = 0; i < n; i++) {
ord[i] = i;
rev[i][i] = 1;
}
ord[n] = n;
/// ---------- First way ----------
printf("\nStd. Gauss'es method...\n\n");
fprintf(output, "\nStd. Gauss'es method:\n\n");
cpy_mat(m, bcp, n);
/// Stage 1
for (k = 0; k < n; k++) {
for(i = k; i < n && !m[i][k]; i++)
;
if (i == n) {
printf("Det = 0\n");
fprintf(output, "Det = 0\n");
exit(1);
}
swap(m + i, m + k, sizeof(*m));
if ((i - k) & 1)
det *= -1;
swap(rev + i, rev + k, sizeof(*m));
for (j = 0; j < n; j++)
rev[k][j] /= m[k][k];
det *= m[k][k];
for (j = n; j >= k; j--)
m[k][j] /= m[k][k];
for (i = k + 1; i < n; i++) {
for (j = 0; j < n; j++)
rev[i][j] -= rev[k][j] * m[i][k];
for (j = n; j >= k; j--)
m[i][j] -= m[k][j] * m[i][k];
}
}
/// Stage 2
for (i = n - 1; i >= 0; i--)
for (j = n - 1; j > i; j--) {
for (k = 0; k < n; k++)
rev[i][k] -= m[i][j] * rev[j][k];
m[i][n] -= m[i][j] * m[j][n];
}
print_ans(output, m, n, ord);
/// ---------- Second way ----------
printf("Gauss'es method with choice main elem...\n\n");
fprintf(output, "\nGauss'es method with choice main elem:\n\n");
cpy_mat(m, bcp, n);
/// Stage 1
for (k = 0; k < n; k++) {
i = k;
for(j = k; j < n; j++)
if (fabs(m[k][ord[j]]) > fabs(m[k][ord[i]]))
i = j;
if (!m[k][ord[i]]) {
printf("Det = 0\n");
return -1;
}
swap(ord + i, ord + k, sizeof(*ord));
for (j = n; j >= k; j--)
m[k][ord[j]] /= m[k][ord[k]];
for (i = k + 1; i < n; i++)
for (j = n; j >= k; j--)
m[i][ord[j]] -= m[k][ord[j]] * m[i][ord[k]];
}
/// Stage 2
for (i = n - 1; i >= 0; i--)
for (j = n - 1; j > i; j--)
m[i][ord[n]] -= m[i][ord[j]] * m[j][ord[n]];
print_ans(output, m, n, ord);
/// ---------- Third way ----------
printf("Zeydel's method...\n\n");
fprintf(output, "\nZeydel's method\nw = 1.000\ne = %lf\n", e);
double x[n], old[n], norm = 0;
for (i = 0; i < n; i++)
x[i] = 0;
for (i = 0; i < n; i++) {
for (j = 0; j <= n; j++) {
m[i][j] = 0;
for (k = 0; k < n; k++) {
m[i][j] += bcp[k][i] * bcp[k][j];
}
}
}
int cnt = 0;
do {
memcpy(old, x, n * sizeof(*old));
cnt++;
for (i = 0; i < n; i++) {
tmp = 0;
for (j = 0; j < n; j++) {
tmp += m[i][j] * x[j];
}
x[i] = x[i] + (m[i][n] - tmp) / m[i][i];
}
printf("%lf\n", norm_diff_vec(x, old, n));
} while (norm_diff_vec(x, old, n) >= e);
fprintf(output, "%d iter.:\n\n", cnt);
for (i = 0; i < n; i++)
fprintf(output, "x%d = %lf\n", i + 1, x[i]);
fprintf(output, "\n");
/// ---------- Fourth way ----------
printf("Relaxation method...\n");
fprintf(output, "\nRelaxation method.\nw = %.3lf\ne = %lf\n", w, e);
for (i = 0; i < n; i++)
x[i] = 0;
cnt = 0;
do {
memcpy(old, x, n * sizeof(*old));
cnt++;
for (i = 0; i < n; i++) {
tmp = 0;
for (j = 0; j < n; j++) {
tmp += m[i][j] * x[j];
}
x[i] = x[i] + w / m[i][i] * (m[i][n] - tmp);
}
} while (norm_diff_vec(x, old, n) >= e);
fprintf(output, "%d iter.:\n\n", cnt);
for (i = 0; i < n; i++)
fprintf(output, "x%d = %lf\n", i + 1, x[i]);
/// -----------------------------
fprintf(output, "\n\nReverse mat:\n\n");
print_mat(output, rev, n, 0);
fprintf(output, "\nDet = %lf\n\n", det);
free_mat(rev, n);
free_mat(bcp, n);
free_mat(m, n);
return 0;
}