labs-3semester/Lab6/lab.c

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

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

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

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

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

// Структура для удобного сбора данных
struct result {
	double lambda;
    int M;
};

// Скалярное произведение
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) {
    for(int i = 0; i < N; i++) {
        v[i] += u[i];
    }
}

// Функция для вычисления сдвига
void shift(mat m, double s) {
    for(int i = 0; i < N; i++) {
        m[i][i] =  m[i][i] - s; 
    }
}

// Умножение матрицы на вектор (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]; 
        }
    }
}

// Функция нормализации
void normalize(vec v){
    double norm = sqrt(dot(v, v));
    for(int i = 0; i < N; i++) {
        v[i] /= norm;
    }
}

// Функция для вычисления ошибки между векторами
double error_vector(vec x, vec x1) {
    double error = 0;
    for (int i = 0; i < N; i++) {
        error = fmax(error, fabs(fabs(x[i]) - fabs(x1[i])));
    }
    return error;
}

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

// Чтение матрицы из файла
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 vec_from_file(FILE* f, vec v) {
	double d;
	for (int i = 0; i < N; 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");
}

struct result solve(mat m, vec x, double shift_v, double eps, vec eigen_vector) {
    struct result res;
    shift(m, shift_v);
    res.M = 0;
    vec y;
    double lambda2;
    normalize(x);
    for (int i = 0; i < 100; i++) {
        apply(m, x, y);
        res.lambda = shift_v + dot(x, y)/dot(x,x);
        res.M++;
        normalize(y);
        copy_vec(y, eigen_vector);
        copy_vec(y, x);
        if (fabs((res.lambda-lambda2)/lambda2) < eps){
            break;
        }
        lambda2 = res.lambda;
    }
    return res;
}


int main(){
    vec x;
    FILE* f;
	if ((f = fopen("matrix.txt", "r")) == NULL) return 1;
    FILE* fs;
	if ((fs = fopen("value_shift.txt", "r")) == NULL) return 1;
    FILE* vf;
    if ((vf = fopen("vector_shift.txt", "r")) == NULL) return 1;

    FILE* file;
    if ((file = fopen("lab_6_error.csv", "w+")) == NULL) return 1;
    mat m;
    vec eigen_vector;
    for (int i = 0; i < 9; i++) {
        for(int j = 0; j < N; j++) {
            x[j] = 1;
        }
        mat_from_file(f, m);
        vec_from_file(vf, eigen_vector);

        double shift_v;
        fscanf(fs, "%lf\n", &shift_v);

        vec eigen_approximate;
        struct result solu = solve(m, x, shift_v, pow(10,-10), eigen_approximate);
        print_vec(eigen_vector);
        print_vec(eigen_approximate);
        printf("\n");
        fprintf(file, "%e, %e, %e\n", pow(10,(i-2)), fabs(pow(10,i-1) - solu.lambda)/pow(10,i-1), error_vector(eigen_vector,eigen_approximate));
    }
    fclose(f);
    fclose(fs);
    fclose(file);
}