labs-3semester/Lab4/lab.c

#define _USE_MATH_DEFINES
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

// Необязательный флаг для дебага (можно удалить)
#define DEBUG 0

// Фиксированная размерность матрицы (обрабатывается препроцессором, просто alias)
#define N 10

//
struct N_eps {
	double eps;
	int I;	
};

// Тип вектора
typedef double vec[N];

// Тип матрицы
typedef double mat[N][N];

// Скалярное произведение
double dot(vec v, vec u) {
    double res = 0;
    for(int i = 0; i < N; i++) res += v[i]*u[i];

    return res;
}

// Прибавление вектора
void add(vec v, vec u, vec result) {
    for(int i = 0; i < N; i++) {
        result[i] = v[i] + u[i];
    }
}

// Вычитание вектора
void sub(vec v, vec u, vec result) {
    for(int i = 0; i < N; i++) {
        result[i] = v[i] - u[i];
    }
}

// Умножение на скаляр
void mul(vec v, double a) {
    for(int i = 0; i < N; i++) {
        v[i] *= a;
    }
}

// Умножение матрицы на вектор (res = m*v)
void apply(mat m, vec v, vec res) {   
    for(int i = 0; i < N; i++) {
        res[i] = 0;
        for(int j = 0; j < N; j++) {
            res[i] += m[i][j] * v[j]; 
        }
    }
}


// Копирование матрицы (copy = m)
void copy_vec(vec copy, vec x) {
    for(int i = 0; i < N; i++) {
        copy[i] = x[i];
    }
}

//Ошибка между двумя векторами
double error(vec a, vec b){
    double error = 0;
    for (int i = 0; i < N; i++) {
		if (error < fabs(a[i] - b[i])) {
			error = fabs(a[i] - b[i]);
		}
	}
    return error;
}

// Чтение матрицы из файла
void mat_from_file(FILE* f, mat m) {
	double d;
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			fscanf(f, "%lf,", &d);
			m[i][j] = d;
		}
	}
}
// Чтение матрицы из файла
void min_max_from_file(FILE* f, vec v) {
	double d;
	for (int i = 0; i < 2; i++) {
		fscanf(f, "%lf,", &d);
		v[i] = d;
	}
}

// Красивая печать матрицы
void print_mat(mat m) {
    for(int j = 0; j < N; j++) {
        printf("-----------|");
    }
    printf("\n");
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            printf("%10.4lf |", m[i][j]);
        }
        printf("\n");
        for(int j = 0; j < N; j++) {
            printf("-----------|");
        }
        printf("\n");
    }
}

// Красивая печать вектора
void print_vec(vec v) {
    for(int j = 0; j < N; j++) {
        printf("%10.4lf |", v[j]);
    }
    printf("\n");
}

void mx_b(mat m, vec x, vec b, vec result){
    apply(m,x,result);
    sub(b,result,result);
}

// Решение СЛАУ методом Ричардсона
struct N_eps solve(mat m, double eps, int M, double min, double max) {
    struct N_eps result;
    result.I = 0;
    vec x = { 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9};
    vec x1;
    vec b;
    vec var;
    vec x_k = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
    apply(m,x_k,b);
    print_vec(b);

    for (int k = 0; k < M; k++) {
        


        double t = cos(M_PI * (2 * (k + 1) - 1) / (2 * M));
        printf("\n t = %lf\n", t);
        double alpha  = 2.0 / ((min + max) + ((max - min)*t));
    
        mx_b(m, x, b, var);
        //print_vec(var);

        mul(var,alpha);
        //print_vec(var);

        add(x,var,x1);
        print_vec(x1);

        if (error(x,x1) < eps){
            result.eps = error(x,x1);
            result.I += 1;
            break;
        }
        result.eps = error(x,x1);
        result.I += 1;
        copy_vec(x,x1);
        print_vec(x);
    }
    return result;
}

//взять из мат лаба
int main(){
    FILE* f;
    FILE* f2;
	if ((f = fopen("lab4_1.txt", "r")) == NULL) return 1;
    if ((f2 = fopen("min_max.txt", "r")) == NULL) return 1;
    mat m;
    vec min_max;
    FILE* file;
    FILE* file2;
    if ((file = fopen("lab_4_cond.csv", "w+")) == NULL) return 1;
    if ((file2 = fopen("lab_4_time.csv", "w+")) == NULL) return 1;
    for (int k = 0; k < 8; k++){
        mat_from_file(f, m);
        min_max_from_file(f2, min_max);
        clock_t begin = clock();
        struct N_eps result = solve(m, pow(10, -5), 1000, min_max[0], min_max[1]);
        clock_t end = clock();
        printf("\n\n\n");
        fprintf(file, "%.5lf, %e\n", pow(10, k), result.eps);
        fprintf(file2, "%.5lf, %e\n", pow(10, k), ((float)(end - begin) / CLOCKS_PER_SEC));
    }
    fclose(f);
    fclose(f2); 
    fclose(file);
    fclose(file2);

    FILE* fbad = fopen("lab_4_error_bad.txt", "r");
    FILE* fminbad = fopen("min_max2.txt", "r");
    FILE* fbade = fopen("lab_4_error_bad.csv", "w+");
    FILE* fbadm = fopen("lab_4_M_bad.csv", "w+");
    mat_from_file(fbad,m);
    min_max_from_file(fminbad, min_max);
    for (int i = 1; i < 6; i++) {
        struct N_eps result = solve(m, pow(10, -i), 1000, min_max[0], min_max[1]);
        fprintf(fbade, "%e, %e\n", pow(10, -i), result.eps);
        fprintf(fbadm, "%e, %d\n", pow(10, -i), result.I);
    }
    fclose(fbad);
    fclose(fminbad);
    fclose(fbade);
    fclose(fbadm);


    FILE* fgood = fopen("lab_4_error_good.txt", "r");\
    FILE* fmingood = fopen("min_max3.txt", "r");
    FILE* fgoode = fopen("lab_4_error_good.csv", "w+");
    FILE* fgoodm = fopen("lab_4_M_good.csv", "w+");
    mat_from_file(fgood,m);
    min_max_from_file(fmingood, min_max);
    for (int i = 1; i < 16; i++) {
        struct N_eps result = solve(m, pow(10, -i), 1000, min_max[0], min_max[1]);
        fprintf(fgoodm, "%e, %d\n", pow(10, -i), result.I);
        fprintf(fgoode, "%e, %e\n", pow(10, -i), result.eps);
    }
    fclose(fgood);
    fclose(fmingood);
    fclose(fgoode);
    fclose(fgoodm);


    FILE* fbad2 = fopen("lab_4_error_bad.txt", "r");
    FILE* fminbad2 = fopen("min_max2.txt", "r");
    FILE* fbadmt = fopen("lab_4_M_bad_tol.csv", "w+");
    mat_from_file(fbad2,m);
    min_max_from_file(fminbad2, min_max);
    for (int i = 0; i < 41; i++) {
        struct N_eps result = solve(m, pow(10, -15), i, min_max[0], min_max[1]);
        fprintf(fbadmt, "%e, %d\n", result.eps, result.I);
    }
    fclose(fbad2);
    fclose(fminbad2);
    fclose(fbadmt);

    FILE* fgood2 = fopen("lab_4_error_good.txt", "r");
    FILE* fmindoof2 = fopen("min_max3.txt", "r");
    FILE* fgoodmt2 = fopen("lab_4_M_good_tol.csv", "w+");
    mat_from_file(fgood2,m);
    min_max_from_file(fmindoof2, min_max);
    for (int i = 0; i < 41; i++) {
        struct N_eps result = solve(m, pow(10, -15), i, min_max[0], min_max[1]);
        fprintf(fgoodmt2, "%e, %d\n", result.eps, i);
    }
    fclose(fgood2);
    fclose(fmindoof2);
    fclose(fgoodmt2);
}