From 2058f7ed6f737f7c34e4fb948f0e4dd1d7efe066 Mon Sep 17 00:00:00 2001 From: AZEN-SGG Date: Wed, 26 Mar 2025 18:30:04 +0300 Subject: [PATCH] Task 5 is done! --- 2025.03.28/05Ex/Makefile | 21 +++++ 2025.03.28/05Ex/a.txt | 5 ++ 2025.03.28/05Ex/array_io.c | 40 ++++++++++ 2025.03.28/05Ex/array_io.h | 11 +++ 2025.03.28/05Ex/init_f.c | 30 +++++++ 2025.03.28/05Ex/init_f.h | 9 +++ 2025.03.28/05Ex/io_status.h | 13 +++ 2025.03.28/05Ex/main.c | 154 ++++++++++++++++++++++++++++++++++++ 2025.03.28/05Ex/matrix.c | 31 ++++++++ 2025.03.28/05Ex/matrix.h | 7 ++ 2025.03.28/05Ex/solve.c | 112 ++++++++++++++++++++++++++ 2025.03.28/05Ex/solve.h | 8 ++ 2025.03.28/05Ex/t.txt | 5 ++ 2025.03.28/05Ex/v.txt | 2 + 14 files changed, 448 insertions(+) create mode 100644 2025.03.28/05Ex/Makefile create mode 100644 2025.03.28/05Ex/a.txt create mode 100644 2025.03.28/05Ex/array_io.c create mode 100644 2025.03.28/05Ex/array_io.h create mode 100644 2025.03.28/05Ex/init_f.c create mode 100644 2025.03.28/05Ex/init_f.h create mode 100644 2025.03.28/05Ex/io_status.h create mode 100644 2025.03.28/05Ex/main.c create mode 100644 2025.03.28/05Ex/matrix.c create mode 100644 2025.03.28/05Ex/matrix.h create mode 100644 2025.03.28/05Ex/solve.c create mode 100644 2025.03.28/05Ex/solve.h create mode 100644 2025.03.28/05Ex/t.txt create mode 100644 2025.03.28/05Ex/v.txt diff --git a/2025.03.28/05Ex/Makefile b/2025.03.28/05Ex/Makefile new file mode 100644 index 0000000..9fe59ac --- /dev/null +++ b/2025.03.28/05Ex/Makefile @@ -0,0 +1,21 @@ +CFLAGS = -mfpmath=sse -fstack-protector-all -W -Wall -Wextra -Wunused \ +-Wempty-body -Wlogical-op -Wold-style-declaration -Wmissing-parameter-type \ +-Wignored-qualifiers -Winit-self -Wshadow -Wtype-limits \ +-Wpointer-arith -Wformat-security -Wmissing-format-attribute -Wformat=1 \ +-Wdeclaration-after-statement -Wbad-function-cast -Wnested-externs \ +-Wmissing-prototypes -Wmissing-declarations -Wold-style-definition \ +-Wcast-align -Werror -pedantic -pedantic-errors -Wfloat-equal \ +-Wwrite-strings -Wno-long-long -std=gnu99 -Wstrict-prototypes \ +-Wmissing-field-initializers -Wpointer-sign -fopenmp -O3 + +TARGET = a05.out +OBJ = main.o solve.o array_io.o init_f.o matrix.o + +%.o: %.c + gcc $(CFLAGS) -c $< -o $@ + +$(TARGET): $(OBJ) + gcc $^ -o $@ -lm -fopenmp + +clean: + rm *.o *.out diff --git a/2025.03.28/05Ex/a.txt b/2025.03.28/05Ex/a.txt new file mode 100644 index 0000000..171f0cf --- /dev/null +++ b/2025.03.28/05Ex/a.txt @@ -0,0 +1,5 @@ +4 1 0 0 0 +1 4 1 0 0 +0 1 4 1 0 +0 0 1 4 1 +0 0 0 1 4 diff --git a/2025.03.28/05Ex/array_io.c b/2025.03.28/05Ex/array_io.c new file mode 100644 index 0000000..f91fbac --- /dev/null +++ b/2025.03.28/05Ex/array_io.c @@ -0,0 +1,40 @@ +#include +#include "array_io.h" + +io_status read_matrix(double *a, int m, int n, const char *name) +{ + int i, j; + FILE *fp; + if (!(fp = fopen(name, "r"))) return ERROR_OPEN; + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + if (fscanf(fp, "%lf", a + i * n + j) != 1) + {fclose(fp); return ERROR_READ;} + fclose(fp); + return SUCCESS; +} + +void print_matrix(const double *a, int m, int n, int p) +{ + int mp = (m > p ? p : m); + int np = (n > p ? p : n); + int i, j; + + for (i = 0; i < mp; i++) + { + for (j = 0; j < np; j++) + printf(" %10.3e", a[i * n + j]); + printf("\n"); + } +} + +void init_matrix(double *a, int m, int n, int k) +{ + double (*q)(int, int, int, int); + double (*f[])(int, int, int, int) = {f1, f2, f3, f4}; + int i, j; + q = f[k-1]; + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + a[i * n + j] = q(m, n, i+1, j+1); +} diff --git a/2025.03.28/05Ex/array_io.h b/2025.03.28/05Ex/array_io.h new file mode 100644 index 0000000..c74c9a0 --- /dev/null +++ b/2025.03.28/05Ex/array_io.h @@ -0,0 +1,11 @@ +#ifndef ARRAY_IO_H +#define ARRAY_IO_H + +#include "io_status.h" +#include "init_f.h" + +io_status read_matrix(double *a, int m, int n, const char *name); +void print_matrix(const double *a, int m, int n, int p); +void init_matrix(double *a, int m, int n, int k); + +#endif diff --git a/2025.03.28/05Ex/init_f.c b/2025.03.28/05Ex/init_f.c new file mode 100644 index 0000000..def39b7 --- /dev/null +++ b/2025.03.28/05Ex/init_f.c @@ -0,0 +1,30 @@ +#include "init_f.h" +#include + +#define MAX(n, m) (n < m ? m : n) + +double f1(int m, int n, int i, int j) +{ + return MAX(n, m) - MAX(i, j) + 1; +} + +double f2(int m, int n, int i, int j) +{ + (void)n; + (void)m; + return MAX(i, j); +} + +double f3(int m, int n, int i, int j) +{ + (void)n; + (void)m; + return abs(i - j); +} + +double f4(int m, int n, int i, int j) +{ + (void)n; + (void)m; + return 1./(i+j-1); +} diff --git a/2025.03.28/05Ex/init_f.h b/2025.03.28/05Ex/init_f.h new file mode 100644 index 0000000..40e985b --- /dev/null +++ b/2025.03.28/05Ex/init_f.h @@ -0,0 +1,9 @@ +#ifndef INIT_F_H +#define INIT_F_H + +double f1(int m, int n, int i, int j); +double f2(int m, int n, int i, int j); +double f3(int m, int n, int i, int j); +double f4(int m, int n, int i, int j); + +#endif diff --git a/2025.03.28/05Ex/io_status.h b/2025.03.28/05Ex/io_status.h new file mode 100644 index 0000000..420a981 --- /dev/null +++ b/2025.03.28/05Ex/io_status.h @@ -0,0 +1,13 @@ +#ifndef IO_STATUS_H +#define IO_STATUS_H + +#define LEN 1234 + +typedef enum _io_status +{ + SUCCESS, + ERROR_OPEN, + ERROR_READ +} io_status; + +#endif diff --git a/2025.03.28/05Ex/main.c b/2025.03.28/05Ex/main.c new file mode 100644 index 0000000..c5984f2 --- /dev/null +++ b/2025.03.28/05Ex/main.c @@ -0,0 +1,154 @@ +#include +#include +#include +#include "array_io.h" +#include "io_status.h" +#include "matrix.h" +#include "solve.h" + +/* ./a.out m n p k_a [filename_a] k_x [filename_x] */ +int main(int argc, char *argv[]) +{ + double t, r1, r2, *a, *x_0, *b, *x, *r; + int n, m, p, k_a, k_x, task = 5; + char *name_a = 0, *name_x = 0; + if (!((argc == 6 || argc == 7 || argc == 8) && + sscanf(argv[1], "%d", &m) == 1 && + sscanf(argv[2], "%d", &n) == 1 && + sscanf(argv[3], "%d", &p) == 1 && + sscanf(argv[4], "%d", &k_a) == 1 && + (k_a >= 0 && k_a <= 4) && + (!(k_a == 0 && argc == 6)) && + ((k_a == 0 && sscanf(argv[6], "%d", &k_x) == 1) || + (k_a != 0 && sscanf(argv[5], "%d", &k_x) == 1)) && + (k_x >= 0 && k_x <= 4) && + (!(k_x == 0 && argc == 6)) && + (!((k_a == 0 && k_x == 0) && argc != 8)))) + { + printf("Usage: %s m n p k_a [filename_a] k_x [filename_x]\n", argv[0]); + return 1; + } + if (argc != 6) + { + int i_x = 6; + if (k_a == 0) {name_a = argv[5];i_x++;} + if (k_x == 0) name_x = argv[i_x]; + } + + a = (double *)malloc((size_t)n * (size_t)n * sizeof(double)); + if (!a) + { + printf("Not enough memory\n"); + return 2; + } + x_0 = (double *)malloc((size_t)n * sizeof(double)); + if (!x_0) + { + free(a); + printf("Not enough memory\n"); + return 2; + } + b = (double *)malloc((size_t)n * sizeof(double)); + if (!b) + { + free(a); + free(x_0); + printf("Not enough memory\n"); + return 2; + } + x = (double *)malloc((size_t)n * sizeof(double)); + if (!x) + { + free(a); + free(x_0); + free(b); + printf("Not enough memory\n"); + return 2; + } + r = (double *)malloc((size_t)n * sizeof(double)); + if (!r) + { + free(a); + free(x_0); + free(b); + free(x); + printf("Not enough memory\n"); + return 2; + } + + if (name_a) + { /* из файла */ + io_status ret; + ret = read_matrix(a, n, n, name_a); + do { + switch (ret) + { + case SUCCESS: + continue; + case ERROR_OPEN: + printf("Cannot open %s\n", name_a); + break; + case ERROR_READ: + printf("Cannot read %s\n", name_a); + break; + } + free(a); + free(x_0); + free(b); + free(x); + return 3; + } while (0); + } else init_matrix(a, n, n, k_a); + + if (name_x) + { + io_status ret; + ret = read_matrix(x_0, n, 1, name_x); + do { + switch (ret) + { + case SUCCESS: + continue; + case ERROR_OPEN: + printf("Cannot open %s\n", name_x); + break; + case ERROR_READ: + printf("Cannot read %s\n", name_x); + break; + } + free(a); + free(x_0); + free(b); + free(x); + return 3; + } while (0); + } else init_matrix(x_0, n, 1, k_x); + + init_vec_b(a, b, n); + + printf("Matrix A:\n"); + print_matrix(a, n, n, p); + printf("Vector x_0:\n"); + print_matrix(x_0, 1, n, p); + printf("Vector b:\n"); + print_matrix(b, 1, n, p); + + t = clock(); + t5_solve(a, x_0, b, x, r, n, m); + t = (clock() - t) / CLOCKS_PER_SEC; + + r1 = get_r1(a, x, b, n); + r2 = get_r2_value(x, n); + + printf("Vector x_m:\n"); + print_matrix(x, 1, n, p); + printf("%s : Task = %d Res1 = %e Res2 = %e Elapsed = %.2f\n", argv[0], task, r1, r2, t); + + free(a); + free(x_0); + free(b); + free(x); + free(r); + + return 0; +} diff --git a/2025.03.28/05Ex/matrix.c b/2025.03.28/05Ex/matrix.c new file mode 100644 index 0000000..f52cc14 --- /dev/null +++ b/2025.03.28/05Ex/matrix.c @@ -0,0 +1,31 @@ +#include "matrix.h" +#include + +void init_vec_b(const double * restrict a, double * restrict b, int n) +{ + #pragma omp parallel for + for (int i = 0; i < n; ++i) + { + double sum = 0; + + #pragma omp simd reduction(+:sum) + for (int k = 1; k < n; k+=2) + sum += a[i * n + k]; + + b[i] = sum; + } +} + +void matvec_mul(int n, const double * restrict A, const double * restrict x, double * restrict x_k) +{ + #pragma omp parallel for + for (int i = 0; i < n; i++) + { + double sum = 0; + #pragma omp simd reduction(+:sum) + for (int j = 0; j < n; j++) + sum += A[i * n + j] * x[j]; + x_k[i] = sum; + } +} + diff --git a/2025.03.28/05Ex/matrix.h b/2025.03.28/05Ex/matrix.h new file mode 100644 index 0000000..66351d7 --- /dev/null +++ b/2025.03.28/05Ex/matrix.h @@ -0,0 +1,7 @@ +#ifndef MATRIX_H +#define MATRIX_H + +void init_vec_b(const double * restrict a, double * restrict b, int n); +void matvec_mul(int n, const double * restrict A, const double * restrict x, double * restrict x_k); + +#endif diff --git a/2025.03.28/05Ex/solve.c b/2025.03.28/05Ex/solve.c new file mode 100644 index 0000000..814e3db --- /dev/null +++ b/2025.03.28/05Ex/solve.c @@ -0,0 +1,112 @@ +#include "solve.h" +#include "matrix.h" +#include + + +void t5_solve(const double * restrict A, double * restrict x_0, const double * restrict b, double * restrict x, double * restrict r, int n, int m) +{ + for (int k = 0; k < m; ++k) + { + double *swap_temp = 0; + double t = 0; + double dot_Dr_r = 0, dot_ADr_Dr = 0; + + #pragma omp parallel for + for (int i = 0; i < n; ++i) + { + double sum = 0; + #pragma omp simd reduction(+:sum) + for (int j = 0; j < n; ++j) + sum += A[i*n + j] * x_0[j]; + r[i] = sum - b[i]; + } + + #pragma omp parallel for reduction(+:dot_Dr_r) + for (int i = 0; i < n; ++i) + { + double d_ri = r[i] / A[i*n +i]; + dot_Dr_r += d_ri * r[i]; + r[i] = d_ri; + } + + #pragma omp parallel for reduction(+:dot_ADr_Dr) + for (int i = 0; i < n; ++i) + { + double sum = 0; + + #pragma omp simd reduction(+:sum) + for (int j = 0; j < n; ++j) + sum += A[i*n + j] * r[j]; + + dot_ADr_Dr += sum * r[i]; + } + + t = dot_Dr_r / dot_ADr_Dr; + + #pragma omp simd + for (int i = 0; i < n; ++i) + x[i] = x_0[i] - r[i]*t; + + swap_temp = x; + x = x_0; + x_0 = swap_temp; + } + + if (m % 2 == 0) // Проверил 100 раз + #pragma omp simd + for (int i = 0; i < n; i++) + { + double temp = x[i]; + x[i] = x_0[i]; + x_0[i] = temp; + } + else + { + double * swap_temp = x; + x = x_0; + x_0 = swap_temp; + } +} + +double get_r1(const double * restrict A, const double * restrict x_k, const double * restrict b, int n) +{ + double norm_r1x_1 = 0; + double residual_norm_1 = 0; + double r1 = 0; + + #pragma omp parallel for reduction(+:residual_norm_1, norm_r1x_1) + for (int i = 0; i < n; ++i) + { + double bi = b[i]; + double sum = 0; + + #pragma omp simd reduction(+:sum) + for (int j = 0; j < n; ++j) + sum += A[i*n + j] * x_k[j]; + + residual_norm_1 += fabs(sum - bi); + norm_r1x_1 += fabs(bi); + } + + r1 = residual_norm_1 / norm_r1x_1; + + return r1; +} + +double get_r2_value(const double * restrict x_k, int n) +{ + double relative_error = 0; + double total_diff = 0; + double template_sum = 0; + + #pragma omp parallel for reduction(+:total_diff, template_sum) + for (int i = 0; i < n; ++i) + { + short int modi = i ^ 1; + total_diff += fabs(x_k[i] - modi); + template_sum += modi; + } + + relative_error = total_diff / template_sum; + return relative_error; +} diff --git a/2025.03.28/05Ex/solve.h b/2025.03.28/05Ex/solve.h new file mode 100644 index 0000000..c4152b8 --- /dev/null +++ b/2025.03.28/05Ex/solve.h @@ -0,0 +1,8 @@ +#ifndef SOLVE_H +#define SOLVE_H + +void t5_solve(const double * restrict A, double * restrict x_0, const double * restrict b, double * restrict x, double * restrict r, int n, int m); +double get_r1(const double * restrict A, const double * restrict x_k, const double * restrict b, int n); +double get_r2_value(const double * restrict x_k, int n); + +#endif diff --git a/2025.03.28/05Ex/t.txt b/2025.03.28/05Ex/t.txt new file mode 100644 index 0000000..229972f --- /dev/null +++ b/2025.03.28/05Ex/t.txt @@ -0,0 +1,5 @@ +0 +0 +0 +0 +0 diff --git a/2025.03.28/05Ex/v.txt b/2025.03.28/05Ex/v.txt new file mode 100644 index 0000000..6ed281c --- /dev/null +++ b/2025.03.28/05Ex/v.txt @@ -0,0 +1,2 @@ +1 +1