diff --git a/2025.04.04/14Ex/main.c b/2025.04.04/14Ex/main.c index 53a4bdd..ffbd2b5 100644 --- a/2025.04.04/14Ex/main.c +++ b/2025.04.04/14Ex/main.c @@ -23,7 +23,7 @@ int main(int argc, char *argv[]) printf("Usage: %s n p k [filename]\n", argv[0]); return 1; } - if (argc == 5) name = argv[4]; + if (argc == 5 && k == 0) name = argv[4]; a = (double *)malloc(n * n * sizeof(double)); if (!a) diff --git a/2025.04.04/dist/Kochubei/Makefile b/2025.04.04/dist/Kochubei/Makefile index cc56208..2c5a187 100644 --- a/2025.04.04/dist/Kochubei/Makefile +++ b/2025.04.04/dist/Kochubei/Makefile @@ -1,12 +1,12 @@ OPTS = -mfpmath=sse -fstack-protector-all -W -Wall -Wextra -Wunused -Wcast-align -Werror -pedantic -pedantic-errors -Wfloat-equal -Wpointer-arith -Wformat-security -Wmissing-format-attribute -Wformat=1 -Wwrite-strings -Wcast-align -Wno-long-long -std=gnu99 -Wstrict-prototypes -Wmissing-prototypes -Wmissing-declarations -Wold-style-definition -Wdeclaration-after-statement -Wbad-function-cast -Wnested-externs -Wmaybe-uninitialized -Wempty-body -Wlogical-op -Wold-style-declaration -Wmissing-parameter-type -Wignored-qualifiers -Winit-self -Wshadow -Wtype-limits -Wstrict-prototypes -Wmissing-field-initializers -Wno-pointer-sign -Wswitch-default -fopenmp -O3 -all: a.exe +all: a.out %.o: %.c gcc -c $(OPTS) $< -%.exe: %.o +%.out: %.o gcc $(OPTS) $^ -o $@ solve.o: solve.c solve.h @@ -14,4 +14,7 @@ add.o: add.c add.h a.o: a.c solve.h add.h -a.exe: a.o solve.o add.o +a.out: a.o solve.o add.o + +clean: + rm *.o *.out diff --git a/2025.04.04/dist/Krivoruchenko_SK/Makefile b/2025.04.04/dist/Krivoruchenko_SK/Makefile new file mode 100644 index 0000000..f3b65cb --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/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 = a.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.04.04/dist/Krivoruchenko_SK/array_io.c b/2025.04.04/dist/Krivoruchenko_SK/array_io.c new file mode 100644 index 0000000..188d677 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/array_io.c @@ -0,0 +1,50 @@ +#include +#include +#include +#include "array_io.h" + +io_status read_matrix(double *a, int n, const char *name) +{ + int i, j; + FILE *fp; + if (!(fp = fopen(name, "r"))) return ERROR_OPEN; + for (i = 0; i < n; 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 n, int p) +{ + int np = (n > p ? p : n); + int i, j; + + for (i = 0; i < np; i++) + { + for (j = 0; j < np; j++) + printf(" %10.3e", a[i * n + j]); + printf("\n"); + } +} + +void init_matrix(double *a, 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 < n; i++) + for (j = 0; j < n; j++) + a[i * n + j] = q(n, n, i+1, j+1); +} + +void init_identity_matrix(double *a, int n) +{ + a = memset(a, 0, n * n * sizeof(double)); + + #pragma omp simd + for (int i = 0; i < n; ++i) + a[i*n + i] = 1.0; +} diff --git a/2025.04.04/dist/Krivoruchenko_SK/array_io.h b/2025.04.04/dist/Krivoruchenko_SK/array_io.h new file mode 100644 index 0000000..77b4508 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/array_io.h @@ -0,0 +1,12 @@ +#ifndef ARRAY_IO_H +#define ARRAY_IO_H + +#include "io_status.h" +#include "init_f.h" + +io_status read_matrix(double *a, int n, const char *name); +void print_matrix(const double *a, int n, int p); +void init_matrix(double *a, int n, int k); +void init_identity_matrix(double *a, int n); + +#endif diff --git a/2025.04.04/dist/Krivoruchenko_SK/init_f.c b/2025.04.04/dist/Krivoruchenko_SK/init_f.c new file mode 100644 index 0000000..925afb5 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/init_f.c @@ -0,0 +1,30 @@ +#include "init_f.h" +#include + +#define MAX(n, m) (n < m ? m : n) + +double f1(int n, int m, int i, int j) +{ + return MAX(n, m) - MAX(i, j) + 1; +} + +double f2(int n, int m, int i, int j) +{ + (void)n; + (void)m; + return MAX(i, j); +} + +double f3(int n, int m, int i, int j) +{ + (void)n; + (void)m; + return abs(i - j); +} + +double f4(int n, int m, int i, int j) +{ + (void)n; + (void)m; + return 1./(i+j-1); +} diff --git a/2025.04.04/dist/Krivoruchenko_SK/init_f.h b/2025.04.04/dist/Krivoruchenko_SK/init_f.h new file mode 100644 index 0000000..691aace --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/init_f.h @@ -0,0 +1,9 @@ +#ifndef INIT_F_H +#define INIT_F_H + +double f1(int n, int m, int i, int j); +double f2(int n, int m, int i, int j); +double f3(int n, int m, int i, int j); +double f4(int n, int m, int i, int j); + +#endif diff --git a/2025.04.04/dist/Krivoruchenko_SK/io_status.h b/2025.04.04/dist/Krivoruchenko_SK/io_status.h new file mode 100644 index 0000000..4ec83a5 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/io_status.h @@ -0,0 +1,13 @@ +#ifndef IO_STATUS_H +#define IO_STATUS_H + +#define SINGULAR 2 + +typedef enum _io_status +{ + SUCCESS, + ERROR_OPEN, + ERROR_READ +} io_status; + +#endif diff --git a/2025.04.04/dist/Krivoruchenko_SK/main.c b/2025.04.04/dist/Krivoruchenko_SK/main.c new file mode 100644 index 0000000..ffbd2b5 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/main.c @@ -0,0 +1,130 @@ +#include +#include +#include +#include +#include "array_io.h" +#include "io_status.h" +#include "matrix.h" +#include "solve.h" + +/* ./a.out n p k [filename] */ +int main(int argc, char *argv[]) +{ + double t, *a, *x, r1 = 0, r2 = 0; + int n, p, k, res = 0, *c, task = 14; + char *name = 0; + + if (!((argc == 4 || argc == 5) && + sscanf(argv[1], "%d", &n) == 1 && + sscanf(argv[2], "%d", &p) == 1 && + sscanf(argv[3], "%d", &k) == 1 && + k >= 0 && k <= 4 && (!(k == 0 && argc != 5)))) + { + printf("Usage: %s n p k [filename]\n", argv[0]); + return 1; + } + if (argc == 5 && k == 0) name = argv[4]; + + a = (double *)malloc(n * n * sizeof(double)); + if (!a) + { + printf("Not enough memory\n"); + return 2; + } + + x = (double *)malloc(n * n * sizeof(double)); + if (!x) + { + free(a); + printf("Not enough memory\n"); + return 2; + } + c = (int *)malloc(n * sizeof(int)); + if (!c) + { + free(a); + free(x); + printf("Not enough memory\n"); + return 2; + } + + if (name) + { /* из файла */ + io_status ret; + ret = read_matrix(a, n, name); + do { + switch (ret) + { + case SUCCESS: + continue; + case ERROR_OPEN: + printf("Cannot open %s\n", name); + break; + case ERROR_READ: + printf("Cannot read %s\n", name); + } + free(a); + free(x); + free(c); + return 3; + } while (0); + } else init_matrix(a, n, k); + + init_identity_matrix(x, n); + + #pragma omp simd + for (int i = 0; i < n; ++i) + c[i] = i; + + printf("Initial matrix:\n"); + print_matrix(a, n, p); + + t = omp_get_wtime(); + res = t14_solve(n, a, x, c); + t = omp_get_wtime() - t; + + if (res == SINGULAR) + { + free(a); + free(x); + free(c); + + printf("The matrix is degenerate\n"); + return 4; + } + + if (name) + { /* из файла */ + io_status ret; + ret = read_matrix(a, n, name); + do { + switch (ret) + { + case SUCCESS: + continue; + case ERROR_OPEN: + printf("Cannot open %s\n", name); + break; + case ERROR_READ: + printf("Cannot read %s\n", name); + } + free(a); + free(x); + free(c); + return 3; + } while (0); + } else init_matrix(a, n, k); + + r1 = get_r1(n, a, x); + r2 = get_r2(n, a, x); + + printf("Inverse matrix:\n"); + print_matrix(x, n, p); + printf("%s : Task = %d Res1 = %e Res2 = %e Elapsed = %.2f K = %d N = %d\n", argv[0], task, r1, r2, t, k, n); + + free(a); + free(x); + free(c); + + return 0; +} diff --git a/2025.04.04/dist/Krivoruchenko_SK/matrix.c b/2025.04.04/dist/Krivoruchenko_SK/matrix.c new file mode 100644 index 0000000..149ac2c --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/matrix.c @@ -0,0 +1,93 @@ +#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 = 0; 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; + } +} + +double get_r1(const int n, const double * restrict A, const double * restrict X) +{ + double maximum = 0; + + if (n > 4000) return 0; + + #pragma omp parallel for reduction(max:maximum) + for (int j = 0; j < n; ++j) + { + double sum = 0; + + #pragma omp simd reduction(+:sum) + for (int i = 0; i < n; ++i) + { + int in = i*n; + double sum_ij = - (i == j); + + #pragma omp simd reduction(+:sum_ij) + for (int k = 0; k < n; ++k) + sum_ij += A[in + k] * X[k*n + j]; + + sum += fabs(sum_ij); + } + + if (maximum < sum) + maximum = sum; + } + + return maximum; +} + +double get_r2(const int n, const double * restrict A, const double * restrict X) +{ + double maximum = 0; + + if (n > 4000) return 0; + + #pragma omp parallel for reduction(max:maximum) + for (int j = 0; j < n; ++j) + { + double sum = 0; + + #pragma omp simd reduction(+:sum) + for (int i = 0; i < n; ++i) + { + int in = i*n; + double sum_ij = - (i == j); + + #pragma omp simd reduction(+:sum_ij) + for (int k = 0; k < n; ++k) + sum_ij += X[in + k] * A[k*n + j]; + + sum += fabs(sum_ij); + } + + if (maximum < sum) + maximum = sum; + } + + return maximum; +} + diff --git a/2025.04.04/dist/Krivoruchenko_SK/matrix.h b/2025.04.04/dist/Krivoruchenko_SK/matrix.h new file mode 100644 index 0000000..72e27c0 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/matrix.h @@ -0,0 +1,9 @@ +#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); +double get_r1(const int n, const double * restrict A, const double * restrict X); +double get_r2(const int n, const double * restrict A, const double * restrict X); + +#endif diff --git a/2025.04.04/dist/Krivoruchenko_SK/solve.c b/2025.04.04/dist/Krivoruchenko_SK/solve.c new file mode 100644 index 0000000..ec62a74 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/solve.c @@ -0,0 +1,216 @@ +#include "solve.h" +#include "io_status.h" +#include +#include +#include "array_io.h" +#include + +// c - changes in rows +int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) +{ + // Проходимся по главным минорам + for (int k = 0; k < n; ++k) { + double maximum = -1.; + int max_i = 0, max_j = 0; + +// printf("\n--------- K = %d ---------\n", k); + + // Ищем максимальный элемент минора + #pragma omp parallel + { + double local_max = -1.; + int loc_i = 0, loc_j = 0; + + #pragma omp for collapse(2) nowait + for (int i = k; i < n; ++i) + for (int j = k; j < n; ++j) + { + double aij = fabs(A[i * n + j]); + if (aij > local_max) { + local_max = aij; + loc_i = i; + loc_j = j; + } + } + + #pragma omp critical + { + if (local_max > maximum) { + maximum = local_max; + max_i = loc_i; + max_j = loc_j; + } + } + } + + + // Если максимальный по модулю элемент равен нулю, значит матрица вырождена + if (fabs(maximum) < DBL_EPSILON) + return SINGULAR; + +// printf("Maximum = %lf for i = %d, j = %d\n", maximum, max_i, max_j); + + // Меняем строки местами, если максимум находится не в k строке + if (max_i != k) + { + int kn = k*n; + int in = max_i*n; + + #pragma omp simd + for (int i = 0; i < k; ++i) + { + int kni = kn+i, ini = in+i; + double swap = X[kni]; + X[kni] = X[ini]; + X[ini] = swap; + } + + #pragma omp parallel for simd + for (int i = k; i < n; ++i) + { + int kni = kn+i, ini = in+i; + double swap = X[kni]; + X[kni] = X[ini]; + X[ini] = swap; + + swap = A[kni]; + A[kni] = A[ini]; + A[ini] = swap; + } + } + +// print_matrix(A, n, n); +// printf("\n"); + + // Меняем столбцы местами + if (max_j != k) + { + int swap_temp = c[max_j]; + c[max_j] = c[k]; + c[k] = swap_temp; + + #pragma omp simd + for (int in = 0; in < n*n; in+=n) + { + double swap = A[in + k]; + A[in + k] = A[in + max_j]; + A[in + max_j] = swap; + } + } + +// print_matrix(A, n, n); +// printf("\n"); + + gauss_inverse(n, k, A, X); + +// print_matrix(A, n, n); +// printf("Inverse matrix:\n"); +// print_matrix(X, n, n); + + } + + gauss_back_substitution(n, A, X); + + // Возвращаем строки назад + for (int k = 0; k < n; ++k) + { + int pnt_cur = c[k]; + + if (pnt_cur != k) + { + int pnt_nxt = 0; + + #pragma omp parallel for + for (int j = 0; j < n; ++j) + { + int loc_cur = pnt_cur; + double temp_cur = X[k*n + j]; + double temp_nxt = 0; + + do { + temp_nxt = X[loc_cur*n + j]; + X[loc_cur*n + j] = temp_cur; + temp_cur = temp_nxt; + + loc_cur = c[loc_cur]; + } while (loc_cur != k); + + X[k*n + j] = temp_cur; + } + + do { + pnt_nxt = c[pnt_cur]; + c[pnt_cur] = pnt_cur; + pnt_cur = pnt_nxt; + } while (pnt_nxt != k); + + c[k] = k; + } + } + + return 0; +} + +void gauss_inverse(const int n, const int k, double * restrict A, double * restrict X) +{ + const int kn = k*n; + const int kk = kn + k; + const double inv_akk = 1./A[kn + k]; + A[kn + k] = 1.; + + for (int ij = kn; ij <= kn+k; ++ij) + { + double xij = X[ij]; + if (fabs(xij) > DBL_EPSILON) X[ij] = xij*inv_akk; + } + + for (int ij = kn + k+1; ij < kn+n; ++ij) + { + double aij = A[ij], xij = X[ij]; + if (fabs(aij) > DBL_EPSILON) A[ij] = aij*inv_akk; + if (fabs(xij) > DBL_EPSILON) X[ij] = xij*inv_akk; + } + + #pragma omp parallel for + for (int i = k+1; i < n; ++i) + { + const int in = i*n; + const double aik = A[in + k]; + A[in + k] = 0; + X[in + k] -= X[kk] * aik; + + #pragma omp simd + for (int j = 0; j < k; ++j) + X[in + j] -= X[kn + j] * aik; + + #pragma omp simd + for (int j = k+1; j < n; ++j) + { + A[in + j] -= A[kn + j] * aik; + X[in + j] -= X[kn + j] * aik; + } + } +} + +// Обратный ход метода Гаусса +void gauss_back_substitution(const int n, double * restrict A, double * restrict X) +{ + // Идём с последней строки и вычитаем её из последующих + for (int k = n-1; k > 0; --k) + { + const int kn = k * n; + + #pragma omp parallel for + for (int i = 0; i < k; ++i) + { + const int in = i*n; + const double aik = A[in + k]; + A[in + k] = 0; + + #pragma omp simd + for (int j = 0; j < n; ++j) + X[in + j] -= X[kn + j] * aik; + } + } +} + diff --git a/2025.04.04/dist/Krivoruchenko_SK/solve.h b/2025.04.04/dist/Krivoruchenko_SK/solve.h new file mode 100644 index 0000000..2e54ac9 --- /dev/null +++ b/2025.04.04/dist/Krivoruchenko_SK/solve.h @@ -0,0 +1,8 @@ +#ifndef SOLVE_H +#define SOLVE_H + +int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c); +void gauss_inverse(const int n, const int k, double * restrict A, double * restrict X); +void gauss_back_substitution(const int n, double * restrict A, double * restrict X); + +#endif