From 27a323466130e02c249a0a7666e2eb0e85d66f93 Mon Sep 17 00:00:00 2001
From: AZEN-SGG <krivoruchenko.sergei@gmail.com>
Date: Wed, 2 Apr 2025 12:47:26 +0300
Subject: [PATCH] =?UTF-8?q?=D0=94=D0=BE=D0=B1=D0=B0=D0=B2=D0=B8=D0=BB=20?=
 =?UTF-8?q?=D0=BA=D0=BE=D0=B4=20=D0=9A=D0=BE=D1=87=D1=83=D0=B1=D0=B5=D1=8F?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 2025.04.04/dist/Kochubei/Makefile |  17 ++
 2025.04.04/dist/Kochubei/a.c      | 130 +++++++++++++++
 2025.04.04/dist/Kochubei/add.c    | 267 ++++++++++++++++++++++++++++++
 2025.04.04/dist/Kochubei/add.h    |  50 ++++++
 2025.04.04/dist/Kochubei/solve.c  |  54 ++++++
 2025.04.04/dist/Kochubei/solve.h  |   8 +
 6 files changed, 526 insertions(+)
 create mode 100644 2025.04.04/dist/Kochubei/Makefile
 create mode 100644 2025.04.04/dist/Kochubei/a.c
 create mode 100644 2025.04.04/dist/Kochubei/add.c
 create mode 100644 2025.04.04/dist/Kochubei/add.h
 create mode 100644 2025.04.04/dist/Kochubei/solve.c
 create mode 100644 2025.04.04/dist/Kochubei/solve.h

diff --git a/2025.04.04/dist/Kochubei/Makefile b/2025.04.04/dist/Kochubei/Makefile
new file mode 100644
index 0000000..cc56208
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/Makefile
@@ -0,0 +1,17 @@
+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
+
+%.o: %.c
+	gcc -c $(OPTS) $<
+
+%.exe: %.o
+	gcc $(OPTS) $^ -o $@
+	
+solve.o: solve.c solve.h
+add.o: add.c add.h
+
+a.o: a.c solve.h add.h
+
+a.exe: a.o solve.o add.o
diff --git a/2025.04.04/dist/Kochubei/a.c b/2025.04.04/dist/Kochubei/a.c
new file mode 100644
index 0000000..990afa2
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/a.c
@@ -0,0 +1,130 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <omp.h>
+#include "solve.h"
+
+#define task 10
+
+int main(int argc, char** argv)
+{
+    int n=0;
+    int p=0;
+    int k=0;
+    char* name;
+    mat A;
+    mat B;
+    mat X;
+    int* r;
+    matrix_status result;
+    double r1=0, r2=0;
+    double t=0;
+    if (argc<4)
+    {
+        printf("Usage: %s n p k [name if k=0]\n", argv[0]);
+        return -1;
+    }
+    if (sscanf(argv[1], "%d", &n)!=1 || sscanf(argv[2], "%d", &p)!=1 || sscanf(argv[3], "%d", &k)!=1)
+    {
+        printf("Usage: %s n p k [name if k=0]\n", argv[0]);
+        return -1;
+    }
+    A=(mat)malloc(n*n*sizeof(double));
+    if (!A)
+    {
+        printf("Not enough memory\n");
+        return -1;
+    }
+    if (k==0)
+    {
+        io_status ret;
+        if (argc<5)
+        {
+            printf("Usage: %s n p k [name if k=0]\n", argv[0]);
+            free(A);
+            return -1;
+        }
+        name=argv[4];
+        ret=read_matrix(A,n,n,name);
+        if (ret!=SUCCESS)
+        {
+            printf("Read error\n");
+            free(A);
+            return -1;
+        }
+    }
+    else
+    {
+        init_matrix(A,n,n,k);
+    }
+    B=(mat)malloc(n*sizeof(double));
+    if (!B)
+    {
+        printf("Not enough memory\n");
+        free(A);
+        return -1;
+    }
+    X=(mat)malloc(n*sizeof(double));
+    if (!X)
+    {
+        printf("Not enough memory\n");
+        free(A);
+        free(B);
+        return -1;
+    }
+    r=(int*)malloc(n*sizeof(int));
+    if (!r)
+    {
+        printf("Not enough memory\n");
+        free(A);
+        free(B);
+        free(X);
+        return -1;
+    }
+    init_vector(A,B,n);
+    printf("Matrix:\n");
+    print_matrix(A,n,n,p);
+    printf("Vector:\n");
+    print_matrix(B,n,1,p);
+    t=omp_get_wtime();
+    result=gauss(A,B,X,n,r);
+    t=(omp_get_wtime()-t);
+    if (result==FAILED)
+    {
+        printf("Matrix is singular.\n");
+        free(A);
+        free(B);
+        free(X);
+        free(r);
+        return 0;
+    }
+    printf("Solution:\n");
+    print_matrix(X,n,1,p);
+    if (k==0)
+    {
+        io_status ret;
+        name=argv[4];
+        ret=read_matrix(A,n,n,name);
+        if (ret!=SUCCESS)
+        {
+            printf("Read error\n");
+            free(A);
+            free(B);
+            free(X);
+            free(r);
+            return -1;
+        }
+    }
+    else
+    {
+        init_matrix(A,n,n,k);
+    }
+    init_vector(A,B,n);
+    r1=R1(A,X,B,n);
+    r2=R2(X,n);
+    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(B);
+    free(X);
+    free(r);
+    return 0;
+}
diff --git a/2025.04.04/dist/Kochubei/add.c b/2025.04.04/dist/Kochubei/add.c
new file mode 100644
index 0000000..9b98f34
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/add.c
@@ -0,0 +1,267 @@
+#include <stdio.h>
+#include <math.h>
+#include "add.h"
+
+#define eps 1e-16
+
+// Вычисление минимума, максимума и модуля от ЦЕЛОГО числа
+// ВАЖНО: Не путать с fmax и т.д, а то снова всё ДЗ улетит в FPE...
+int max(int a, int b)
+{
+    if (a>=b) return a;
+    else return b;
+}
+
+int min(int a, int b)
+{
+    if (a<=b) return a;
+    else return b;
+}
+
+int abs(int a)
+{
+    if (a>=0) return 0;
+    else return -a;
+}
+
+// Вычисление элемента матрицы
+// k - номер формулы, n, m - размеры матрицы, i, j - индексы элемента
+// Значения k=5 и k=6 не описаны в условиях и нужны исключительно для большего количества тестов
+double f(int k, int n, int m, int i, int j)
+{
+    i=i+1;
+    j=j+1;
+    if (k==1)
+    {
+        return max(n,m)-max(i,j)+1;
+    }
+    else if (k==2)
+    {
+        return max(i,j);
+    }
+    else if (k==3)
+    {
+        return abs(i-j);
+    }
+    // Матрица Гильберта
+    else if (k==4)
+    {
+        return 1.0/(i+j-1);
+    }
+
+    // Эксперименты с другими матрицами
+
+    // Матрица Лемера
+    else if (k==5)
+    {
+        return (min(i,j))*1.0/max(i,j);
+    }
+    // Матрица Редхеффера
+    else if (k==6)
+    {
+        if (j==1 || j%i==0) return 1;
+        else return 0;
+    }
+    else return 0;
+}
+
+// Чтение матрицы
+// A - массив выделенный под матрицу, n, m - размеры матрицы, name - имя файла
+io_status read_matrix(restrict mat A, int n, int m, char* name)
+{
+    FILE* f=fopen(name, "rt");
+    if (!f) return ERROR_OPEN;
+    for (int i=0; i<n*m; i++)
+    {
+        if (fscanf(f, "%lf", &A[i])!=1)
+        {
+            fclose(f);
+            return ERROR_READ;
+        }
+    }
+    fclose(f);
+    return SUCCESS;
+}
+
+// Вывод матрицы
+// A - матрица, n, m - размеры матрицы, p - максимальное количество строк и столбцов для вывода
+void print_matrix(restrict mat A, int n, int m, int p)
+{
+    for (int i=0; i<n && i<p; i++)
+    {
+        for (int j=0; j<m && j<p; j++)
+        {
+            printf("%10.3e ", A[i*m+j]);
+        }
+        printf("\n");
+    }
+}
+
+// Инициализация матрицы
+// A - матрица, n, m - размеры матрицы, k - номер формулы
+// При нулевом k из main должна вызываться функция read_matrix()
+void init_matrix(restrict mat A, int n, int m, int k)
+{
+    for (int i=0; i<n; i++)
+    {
+        for (int j=0; j<m; j++)
+        {
+            A[i*m+j]=f(k, n, m, i, j);
+        }
+    }
+}
+
+// Инициализация вектора
+// A - заполненная матрица размера n*n, V - вектор, n - длина вектора
+void init_vector(restrict mat A, restrict mat V, int n)
+{
+    for (int i=0; i<n; i++)
+    {
+        double sum=0;
+        for (int k=0; k<=(n-1)/2; k++)
+        {
+            sum=sum+A[i*n+2*k];
+        }
+        V[i]=sum;
+    }
+}
+
+// Норма вектора
+// V - вектор, n - длина вектора
+// Вычисляет норму как сумму модулей компонент вектора
+double vector_norm(restrict mat V, int n)
+{
+    double sum=0;
+    for (int i=0; i<n; i++)
+    {
+        sum=sum+fabs(V[i]);
+    }
+    return sum;
+}
+
+// Вычисление невязки (r1)
+// A - матрица, X - вектор численно найденного решения, B - свободный столбец, n - размер матрицы и длина вектора
+// ВАЖНО: после вызова портит вектор B
+double R1(restrict mat A, restrict mat X, restrict mat B, int n)
+{
+    double norm1=0;
+    double norm2=0;
+    norm2=vector_norm(B, n);
+    if (0<=norm2 && norm2<=0) return 1e300;
+    for (int i=0; i<n; i++)
+    {
+        double sum=0;
+        for (int j=0; j<n; j++)
+        {
+            sum=sum+A[i*n+j]*X[j];
+        }
+        B[i]=sum-B[i];
+    }
+    norm1=vector_norm(B,n);
+    return norm1/norm2;
+}
+
+// Вычисление r2
+// X - вектор, n - его длина
+double R2(restrict mat X, int n)
+{
+    double sum1=0;
+    double sum2=0;
+    if (n<=0) return 0;
+    for (int i=1; i<=n; i++)
+    {
+        sum1=sum1+fabs(X[i-1]-(i&1));
+    }
+    for (int i=1; i<=n; i++)
+    {
+        sum2=sum2+(i&1);
+    }
+    return sum1/sum2;
+}
+
+// Проверка на равенство элемента матрицы нулю, относительно ее нормы
+// x - элемент матрицы, norm - норма матрицы
+int equals_zero(double x, double norm)
+{
+    if (fabs(x)<eps*norm) return 1;
+    else return 0;
+}
+
+// Норма матрицы A_inf
+// A - квадратная матрица, n - ее размер
+double matrix_norm(restrict mat A, int n)
+{
+    double norm=0;
+    for (int i=0; i<n; i++)
+    {
+        double sum=0;
+        for (int j=0; j<n; j++)
+        {
+            sum=sum+fabs(A[i*n+j]);
+        }
+        if (sum>norm) norm=sum;
+    }
+    return norm;
+}
+
+// Элементарное преобразование для метода Гаусса
+// A - квадратная матрица, n - ее размер, k - номер строки
+// ВАЖНО: предполагается, что элемент A[k*n+k] уже выбран ненулевым
+// ВАЖНО: предполагается, что элементы, стоящие левее A[k*n+k], уже являются нулевыми
+void divide_row(restrict mat A, int n, int k)
+{
+    double el=A[k*n+k];
+    for (int i=k+1; i<n; i++)
+    {
+        A[k*n+i]=A[k*n+i]/el;
+    }
+    A[k*n+k]=1;
+}
+
+int find_max(restrict mat A, int n, int k)
+{
+    int res=k;
+    double mx=fabs(A[k*n+k]);
+    for (int i=k+1; i<n; i++)
+    {
+        double tmp=fabs(A[k*n+i]);
+        if (tmp>mx)
+        {
+            res=i;
+            mx=tmp;
+        }
+    }
+    return res;
+}
+
+void subtract_row(restrict mat A, int n, int a, int b)
+{
+    double el=A[a*n+b];
+    for (int i=b+1; i<n; i++)
+    {
+        A[a*n+i]=A[a*n+i]-A[b*n+i]*el;
+    }
+    A[a*n+b]=0;
+}
+
+void swap(double* a, double* b)
+{
+    double tmp=*a;
+    *a=*b;
+    *b=tmp;
+}
+
+void swap_int(int* a, int* b)
+{
+    int tmp=*a;
+    *a=*b;
+    *b=tmp;
+}
+
+void swap_columns(restrict mat A, int n, int a, int b)
+{
+    for (int i=0; i<n; i++)
+    {
+        swap(A+i*n+a, A+i*n+b);
+    }
+}
diff --git a/2025.04.04/dist/Kochubei/add.h b/2025.04.04/dist/Kochubei/add.h
new file mode 100644
index 0000000..8a38c1d
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/add.h
@@ -0,0 +1,50 @@
+typedef enum io_status_ {
+SUCCESS,
+ERROR_OPEN,
+ERROR_READ,
+} io_status;
+
+typedef enum matrix_status_ {
+SOLVED,
+FAILED,
+} matrix_status;
+
+typedef double* mat;
+
+int max(int a, int b);
+
+int min(int a, int b);
+
+int abs(int a);
+
+double f(int k, int n, int m, int i, int j);
+
+io_status read_matrix(restrict mat A, int n, int m, char* name);
+
+void print_matrix(restrict mat A, int n, int m, int p);
+
+void init_matrix(restrict mat A, int n, int m, int k);
+
+void init_vector(restrict mat A, restrict mat V, int n);
+
+double vector_norm(restrict mat V, int n);
+
+double R1(restrict mat A, restrict mat X, restrict mat B, int n);
+
+double R2(restrict mat X, int n);
+
+int equals_zero(double x, double norm);
+
+double matrix_norm(restrict mat A, int n);
+
+void divide_row(restrict mat A, int n, int k);
+
+int find_max(restrict mat A, int n, int k);
+
+void subtract_row(restrict mat A, int n, int a, int b);
+
+void swap(double* a, double* b);
+
+void swap_columns(restrict mat A, int n, int a, int b);
+
+void swap_int(int* a, int* b);
diff --git a/2025.04.04/dist/Kochubei/solve.c b/2025.04.04/dist/Kochubei/solve.c
new file mode 100644
index 0000000..b1eefae
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/solve.c
@@ -0,0 +1,54 @@
+#include "solve.h"
+#include <stdio.h>
+
+void init_array(int* r, int n)
+{
+    for (int i=0; i<n; i++)
+    {
+        r[i]=i;
+    }
+}
+
+matrix_status gauss_step(restrict mat A, restrict mat B, int n, int k, int* r, double norm)
+{
+    int max_pos=find_max(A, n, k);
+    if (max_pos!=k)
+    {
+        swap_columns(A, n, k, max_pos);
+        swap_int(r+k, r+max_pos);
+    }
+    if (equals_zero(A[k*n+k], norm))
+    {
+        return FAILED;
+    }
+    B[k]=B[k]/A[k*n+k];
+    divide_row(A, n, k);
+    #pragma omp parallel for
+    for (int i=k+1; i<n; i++)
+    {
+        double el=A[i*n+k];
+        B[i]=B[i]-el*B[k];
+        subtract_row(A, n, i, k);
+    }
+    return SOLVED;
+}
+
+matrix_status gauss(restrict mat A, restrict mat B, restrict mat X, int n, int* r)
+{
+    double norm=matrix_norm(A,n);
+    init_array(r,n);
+    for (int i=0; i<n; i++)
+    {
+        if (gauss_step(A, B, n, i, r, norm)==FAILED) return FAILED;
+    }
+    for (int i=n-1; i>=0; i--)
+    {
+        double sum=0;
+        for (int j=n-1; j>i; j--)
+        {
+            sum=sum+A[i*n+j]*X[r[j]];
+        }
+        X[r[i]]=B[i]-sum;
+    }
+    return SOLVED;
+}
diff --git a/2025.04.04/dist/Kochubei/solve.h b/2025.04.04/dist/Kochubei/solve.h
new file mode 100644
index 0000000..a33c74b
--- /dev/null
+++ b/2025.04.04/dist/Kochubei/solve.h
@@ -0,0 +1,8 @@
+#include "add.h"
+
+void init_array(int* r, int n);
+
+matrix_status gauss_step(restrict mat A, restrict mat B, int n, int k, int* r, double norm);
+
+matrix_status gauss(restrict mat A, restrict mat B, restrict mat X, int n, int* r);
+