diff --git a/2025.04.04/14Ex/array_io.c b/2025.04.04/14Ex/array_io.c index c40f8e1..8489429 100644 --- a/2025.04.04/14Ex/array_io.c +++ b/2025.04.04/14Ex/array_io.c @@ -1,29 +1,21 @@ #include #include #include "array_io.h" -#include "io_status.h" -io_status read_matrix(double *a, int n, const char *name, bool transpose) +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 (transpose) - { - if (fscanf(fp, "%lf", a + j * n + i) != 1) - {fclose(fp); return ERROR_READ;} - } else - if (fscanf(fp, "%lf", a + i * n + j) != 1) - {fclose(fp); return ERROR_READ;} - } + 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, bool transpose) +void print_matrix(const double *a, int n, int p) { int np = (n > p ? p : n); int i, j; @@ -31,17 +23,12 @@ void print_matrix(const double *a, int n, int p, bool transpose) for (i = 0; i < np; i++) { for (j = 0; j < np; j++) - { - if (transpose) - printf(" %10.3e", a[j * n + i]); - else - printf(" %10.3e", a[i * n + j]); - } + printf(" %10.3e", a[i * n + j]); printf("\n"); } } -void init_matrix(double *a, int n, int k, bool transpose) +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}; @@ -49,15 +36,9 @@ void init_matrix(double *a, int n, int k, bool transpose) q = f[k-1]; for (i = 0; i < n; i++) for (j = 0; j < n; j++) - { - if (transpose) - a[i * n + j] = q(n, n, j+1, i+1); - else - a[i * n + j] = q(n, n, j+1, i+1); - } + 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)); @@ -66,12 +47,12 @@ void init_identity_matrix(double *a, int n) a[i*n + i] = 1.0; } -int read_or_init_matrix(double *a, char *name, int k, int n, bool transpose) +int read_or_init_matrix(double *a, char *name, int k, int n) { if (name) { /* из файла */ io_status ret; - ret = read_matrix(a, n, name, transpose); + ret = read_matrix(a, n, name); do { switch (ret) { @@ -85,7 +66,7 @@ int read_or_init_matrix(double *a, char *name, int k, int n, bool transpose) } return 3; } while (0); - } else init_matrix(a, n, k, transpose); + } else init_matrix(a, n, k); return 0; } diff --git a/2025.04.04/14Ex/array_io.h b/2025.04.04/14Ex/array_io.h index 75dc475..6d75304 100644 --- a/2025.04.04/14Ex/array_io.h +++ b/2025.04.04/14Ex/array_io.h @@ -1,14 +1,13 @@ #ifndef ARRAY_IO_H #define ARRAY_IO_H -#include #include "io_status.h" #include "init_f.h" -io_status read_matrix(double *a, int n, const char *name, bool transpose); -void print_matrix(const double *a, int n, int p, bool transpose); -void init_matrix(double *a, int n, int k, bool transpose); +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); -int read_or_init_matrix(double *a, char *name, int k, int n, bool transpose); +int read_or_init_matrix(double *a, char *name, int k, int n); #endif diff --git a/2025.04.04/14Ex/main.c b/2025.04.04/14Ex/main.c index fe2ccfc..53ba17d 100644 --- a/2025.04.04/14Ex/main.c +++ b/2025.04.04/14Ex/main.c @@ -1,6 +1,5 @@ #include #include -#include #include #include "array_io.h" #include "io_status.h" @@ -48,7 +47,7 @@ int main(int argc, char *argv[]) return 2; } - res = read_or_init_matrix(a, name, k, n, true); + res = read_or_init_matrix(a, name, k, n); if (res) { free(a); @@ -64,12 +63,12 @@ int main(int argc, char *argv[]) c[i] = i; printf("Initial matrix:\n"); - print_matrix(a, n, p, true); + print_matrix(a, n, p); t = clock(); res = t14_solve(n, a, x, c); t = (clock() - t) / CLOCKS_PER_SEC; - + if (res == SINGULAR) { free(a); @@ -80,7 +79,7 @@ int main(int argc, char *argv[]) return 4; } - res = read_or_init_matrix(a, name, k, n, false); + res = read_or_init_matrix(a, name, k, n); if (res) { free(a); @@ -89,12 +88,12 @@ int main(int argc, char *argv[]) return res; } - + r1 = get_r1(n, a, x); r2 = get_r2(n, a, x); printf("Inverse matrix:\n"); - print_matrix(x, n, p, false); + 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); diff --git a/2025.04.04/14Ex/matrix.c b/2025.04.04/14Ex/matrix.c index ce8f0ba..ddea7e4 100644 --- a/2025.04.04/14Ex/matrix.c +++ b/2025.04.04/14Ex/matrix.c @@ -19,7 +19,6 @@ void matvec_mul(int n, const double * restrict A, const double * restrict x, dou for (int i = 0; i < n; i++) { double sum = 0; - for (int j = 0; j < n; j++) sum += A[i * n + j] * x[j]; x_k[i] = sum; @@ -94,6 +93,6 @@ double get_matrix_norm(const int n, const double * restrict A) if (sum > norm) norm = sum; } - + return norm; } diff --git a/2025.04.04/14Ex/solve.c b/2025.04.04/14Ex/solve.c index 768d835..dd583ef 100644 --- a/2025.04.04/14Ex/solve.c +++ b/2025.04.04/14Ex/solve.c @@ -1,15 +1,19 @@ #include "solve.h" #include "io_status.h" -#include -#include #include "array_io.h" #include "matrix.h" +#include +#include #include // c - changes in rows int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) { double norm = get_matrix_norm(n, A); + double eps = DBL_EPSILON*norm; + + if (norm < DBL_EPSILON) + return SINGULAR; // Проходимся по главным минорам for (int k = 0; k < n; ++k) { @@ -26,11 +30,13 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) max_i = i; max_j = j; } - } - + } + +// printf("\n------- K = %d -------\n", k); +// printf("Maximum = %lf i = %d j = %d\n", maximum, max_i, max_j); // Если максимальный по модулю элемент равен нулю, значит матрица вырождена - if (fabs(maximum) < DBL_EPSILON * norm) + if (fabs(maximum) < eps) return SINGULAR; // Меняем строки местами, если максимум находится не в k строке @@ -67,93 +73,79 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) c[max_j] = c[k]; c[k] = swap_temp; - for (int in = 0; in < n*n; in+=n) + for (int i = 0; i < n; i++) { + const int in = i*n; double swap = A[in + k]; A[in + k] = A[in + max_j]; A[in + max_j] = swap; } } + +// printf("BEFORE GAUSS\n"); +// printf("Original matrix:\n"); +// print_matrix(A, n, n); +// printf("Inverse matrix:\n"); +// print_matrix(X, n, n); gauss_inverse(n, k, A, X); + +// printf("AFTER GAUSS\n"); +// printf("Original matrix:\n"); +// 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; + const int kn = k*n; + int i = c[k]; + while (i != k) + { + const int in = i*n; + const int swap_int = c[i]; + c[i] = i; + i = swap_int; + 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; + double swap_temp = X[in+j]; + X[in+j] = X[kn+j]; + X[kn+j] = swap_temp; } - - 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.; + const double inv_akk = 1./A[kk]; - 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; - } + for (int ij = kk+1; ij < kn+n; ij++) + A[ij] *= inv_akk; + + for (int ij = kn; ij < kn+n; ij++) + X[ij] *= inv_akk; 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; - for (int j = 0; j < k; ++j) - X[in + j] -= X[kn + j] * aik; + for (int ij = in+k+1, kj = kk+1; ij < in+n; ij++, kj++) + A[ij] -= A[kj] * aik; - for (int j = k+1; j < n; ++j) - { - A[in + j] -= A[kn + j] * aik; - X[in + j] -= X[kn + j] * aik; - } + for (int ij = in, kj = kn; kj < kn + n; ij++, kj++) + X[ij] -= X[kj] * aik; } } @@ -169,11 +161,16 @@ void gauss_back_substitution(const int n, double * restrict A, double * restrict { const int in = i*n; const double aik = A[in + k]; - A[in + k] = 0; for (int j = 0; j < n; ++j) X[in + j] -= X[kn + j] * aik; } + +// printf("\n------- K = %d -------\n", k); +// printf("Original matrix:\n"); +// print_matrix(A, n, n); +// printf("Inverse matrix:\n"); +// print_matrix(X, n, n); } } diff --git a/2025.04.04/dist/Krivoruchenko_SK/solve.c b/2025.04.04/dist/Krivoruchenko_SK/solve.c index b5c1a76..7180473 100644 --- a/2025.04.04/dist/Krivoruchenko_SK/solve.c +++ b/2025.04.04/dist/Krivoruchenko_SK/solve.c @@ -36,7 +36,7 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) // printf("Maximum = %lf i = %d j = %d\n", maximum, max_i, max_j); // Если максимальный по модулю элемент равен нулю, значит матрица вырождена - if (maximum < eps) + if (fabs(maximum) < eps) return SINGULAR; // Меняем строки местами, если максимум находится не в k строке @@ -45,17 +45,22 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) int kn = k*n; int in = max_i*n; - for (int mj = in, kj = kn; mj < in+n; ++mj, ++kj) + for (int i = 0; i < k; ++i) { - double swap = X[mj]; - X[mj] = X[kj]; - X[kj] = swap; + int kni = kn+i, ini = in+i; + double swap = X[kni]; + X[kni] = X[ini]; + X[ini] = swap; } for (int i = k; i < n; ++i) { int kni = kn+i, ini = in+i; - double swap = A[kni]; + double swap = X[kni]; + X[kni] = X[ini]; + X[ini] = swap; + + swap = A[kni]; A[kni] = A[ini]; A[ini] = swap; } @@ -82,7 +87,7 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) // printf("Inverse matrix:\n"); // print_matrix(X, n, n); - gauss_inverse(n, k, A, X, eps); + gauss_inverse(n, k, A, X); // printf("AFTER GAUSS\n"); // printf("Original matrix:\n"); @@ -92,9 +97,30 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) } gauss_back_substitution(n, A, X); - - // Возвращаем строки назад + for (int k = 0; k < n; ++k) + { + const int kn = k*n; + int i = c[k]; + + while (i != k) + { + const int in = i*n; + const int swap_int = c[i]; + c[i] = i; + i = swap_int; + + for (int ij = in, kj = kn; ij < in+n; ++ij, ++kj) + { + double swap_temp = X[ij]; + X[ij] = X[kj]; + X[kj] = swap_temp; + } + } + } + + // Возвращаем строки назад +/* for (int k = 0; k < n; ++k) { int pnt_cur = c[k]; @@ -127,41 +153,24 @@ int t14_solve(int n, double * restrict A, double * restrict X, int * restrict c) c[k] = k; } - } + }*/ return 0; } // Прямой ход Го ----- йда -void gauss_inverse(const int n, const int k, double * restrict A, double * restrict X, double eps) +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[kk]; - - if (eps > DBL_EPSILON) - eps = DBL_EPSILON; - // Делим на A[k][k] - for (int kj = kk+1; kj < kn+n; kj++) - A[kj] *= inv_akk; - - // Меняем обратную матрицу - for (int j = 0; j < n; ++j) - { - const int BS = 32; - double xkj = X[kn + j]; - if (fabs(xkj) <= eps) - continue; - - xkj *= inv_akk; - X[kn + j] = xkj; - - for (int i_block = k + 1; i_block < n; i_block += BS) - for (int i = i_block; i < i_block + BS && i < n; ++i) - X[i * n + j] -= xkj * A[i * n + k]; - } - + for (int ij = kk+1; ij < kn+n; ij++) + A[ij] *= inv_akk; + + for (int ij = kn; ij < kn+n; ij++) + X[ij] *= inv_akk; + for (int i = k+1; i < n; ++i) { const int in = i*n; @@ -169,6 +178,9 @@ void gauss_inverse(const int n, const int k, double * restrict A, double * restr for (int ij = in+k+1, kj = kk+1; ij < in+n; ij++, kj++) A[ij] -= A[kj] * aik; + + for (int ij = in, kj = kn; kj < kn + n; ij++, kj++) + X[ij] -= X[kj] * aik; } } diff --git a/2025.04.04/dist/Krivoruchenko_SK/solve.h b/2025.04.04/dist/Krivoruchenko_SK/solve.h index a0f0e1c..2e54ac9 100644 --- a/2025.04.04/dist/Krivoruchenko_SK/solve.h +++ b/2025.04.04/dist/Krivoruchenko_SK/solve.h @@ -2,7 +2,7 @@ #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, double eps); +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