Merge pull request #33 from nikitavlaev/Task01

Task01 Vlaev

CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.15)
+project(Task01 C)
+
+set(CMAKE_C_STANDARD 11)
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
+
+add_executable(Task01 main.c gauss_algos.h gauss_algos.c)
+
+find_package(GSL REQUIRED)    # See below (2)
+target_link_libraries(Task01 GSL::gsl GSL::gslcblas)

gauss_algos.c

@@ -0,0 +1,75 @@
+#include "gauss_algos.h"
+
+double *GSL_gauss(int n, double *a, double *b) {
+    gsl_matrix_view gsl_a = gsl_matrix_view_array(a, n, n);
+    gsl_vector_view gsl_b = gsl_vector_view_array(b, n);
+    gsl_vector *gsl_x = gsl_vector_alloc(n);
+    int s;
+    gsl_permutation *p = gsl_permutation_alloc(n);
+    gsl_linalg_LU_decomp(&gsl_a.matrix, p, &s);
+    gsl_linalg_LU_solve(&gsl_a.matrix, p, &gsl_b.vector, gsl_x);
+    gsl_permutation_free(p);
+    double *x = malloc(n * sizeof(double));
+    for (int i = 0; i < n; i++)
+        x[i] = gsl_vector_get(gsl_x, i);
+    gsl_vector_free(gsl_x);
+    return x;
+}
+
+double *seq_gauss(int n, double *a, double *b) {
+    double *x = malloc(sizeof(*x) * n);
+    // Прямой ход -- O(n^3)
+    for (int k = 0; k < n - 1; k++) {
+        // Исключение x_i из строк k+1...n-1
+        double pivot = a[k * n + k];
+        for (int i = k + 1; i < n; i++) {
+            // Из уравнения (строки) i вычитается уравнение k
+            double heads_ratio = a[i * n + k] / pivot;
+            for (int j = k; j < n; j++) {
+                a[i * n + j] -= heads_ratio * a[k * n + j];
+            }
+            b[i] -= heads_ratio * b[k];
+        }
+    }
+    // Обратный ход -- O(n^2)
+    for (int k = n - 1; k >= 0; k--) {
+        x[k] = b[k];
+        for (int i = k + 1; i < n; i++)
+            x[k] -= a[k * n + i] * x[i];
+        x[k] /= a[k * n + k];
+    }
+    return x;
+}
+
+double *parallel_gauss(int n, double *a, double *b) {
+    double *x = malloc(sizeof(*x) * n);
+    // Прямой ход -- O(n^3)
+    int i, j, k;
+    omp_set_num_threads(omp_get_max_threads());
+    for (k = 0; k < n - 1; k++) {
+        // Исключение x_i из строк k+1...n-1
+        double pivot = a[k * n + k];
+
+//#pragma omp parallel for shared(n, a, b, x, pivot, k) private(i, j) default(none) schedule(dynamic)
+        for (i = k + 1; i < n; i++) {
+            // Из уравнения (строки) i вычитается уравнение k
+            double heads_ratio = a[i * n + k] / pivot;
+#pragma omp simd
+            for (j = k; j < n; j++) {
+                a[i * n + j] -= a[k * n + j] * heads_ratio;
+            }
+            b[i] -= heads_ratio * b[k];
+        }
+    }
+
+    // Обратный ход -- O(n^2)
+    for (k = n - 1; k >= 0; k--) {
+        x[k] = b[k];
+#pragma omp simd
+        for (j = k + 1; j < n; j++)
+            x[k] -= a[k * n + j] * x[j];
+        x[k] /= a[k * n + k];
+    }
+    return x;
+}
+

gauss_algos.h

@@ -0,0 +1,13 @@
+#ifndef TASK01_GAUSS_ALGOS_H
+#define TASK01_GAUSS_ALGOS_H
+    #pragma once
+    #include <omp.h>
+    #include <stdlib.h>
+    #include <gsl/gsl_linalg.h>
+    double *GSL_gauss(int n, double *a, double *b);
+
+    double *seq_gauss(int n, double *a, double *b);
+
+    double *parallel_gauss(int n, double *a, double *b);
+
+#endif //TASK01_GAUSS_ALGOS_H

gauss_results.txt

@@ -0,0 +1,28 @@
+Gaussian Elimination (sequential): n 10, mean-time (sec) 0.000002
+Gaussian Elimination (gsl): n 10, mean-time (sec) 0.000003
+Gaussian Elimination (parallel): n 10, mean-time (sec) 0.000002
+
+Gaussian Elimination (sequential): n 100, mean-time (sec) 0.001078
+Gaussian Elimination (gsl): n 100, mean-time (sec) 0.000243
+Gaussian Elimination (parallel): n 100, mean-time (sec) 0.001081
+
+Gaussian Elimination (sequential): n 200, mean-time (sec) 0.008343
+Gaussian Elimination (gsl): n 200, mean-time (sec) 0.001807
+Gaussian Elimination (parallel): n 200, mean-time (sec) 0.008964
+
+Gaussian Elimination (sequential): n 500, mean-time (sec) 0.129559
+Gaussian Elimination (gsl): n 500, mean-time (sec) 0.030214
+Gaussian Elimination (parallel): n 500, mean-time (sec) 0.135428
+
+Gaussian Elimination (sequential): n 1000, mean-time (sec) 1.239578
+Gaussian Elimination (gsl): n 1000, mean-time (sec) 0.286485
+Gaussian Elimination (parallel): n 1000, mean-time (sec) 1.174270
+
+Gaussian Elimination (sequential): n 2000, mean-time (sec) 10.298086
+Gaussian Elimination (gsl): n 2000, mean-time (sec) 2.707311
+Gaussian Elimination (parallel): n 2000, mean-time (sec) 11.302454
+
+Gaussian Elimination (sequential): n 4000, mean-time (sec) 84.347705
+Gaussian Elimination (gsl): n 4000, mean-time (sec) 23.611058
+Gaussian Elimination (parallel): n 4000, mean-time (sec) 80.878901
+

main.c

@@ -0,0 +1,78 @@
+#include <stdio.h>
+#include "gauss_algos.h"
+#include <time.h>
+
+double *generate_matrix(int n) {
+    double *a = malloc(sizeof(*a) * n * n);
+    for (int i = 0; i < n; i++) { // Инициализация
+        // srand(i * (n + 1));
+        for (int j = 0; j < n; j++)
+            a[i * n + j] = rand() % 100 + 1;
+    }
+    return a;
+}
+
+double *generate_column(int n) {
+    double *b = malloc(sizeof(*b) * n);
+    for (int i = 0; i < n; i++) {
+        b[i] = rand() % 100 + 1;
+    }
+    return b;
+}
+
+int main(int argc, char *argv[]) {
+
+    char need_cmp = 0;
+    if (argc < 4) {
+        exit(1);
+    }
+    int n = atoi(argv[1]);
+    int exp_n = atoi(argv[2]);
+    char *algo_name = argv[3];
+
+    clock_t start, end;
+    double clock_sum;
+    double *(*algo)(int, double *, double *);
+    switch (algo_name[0]) {
+        case 's' : {
+            algo = &seq_gauss;
+            need_cmp = 1;
+            break;
+        }
+        case 'g' : {
+            algo = &GSL_gauss;
+            break;
+        }
+        case 'p' : {
+            algo = &parallel_gauss;
+            need_cmp = 1;
+            break;
+        }
+    }
+    clock_sum = 0;
+    for (int j = 0; j < exp_n; j++) {
+        double *a = generate_matrix(n);
+        double *b = generate_column(n);
+        start = clock();
+        double *x = algo(n, a, b);
+        end = clock();
+        clock_sum += ((double)(end - start)) / CLOCKS_PER_SEC;
+
+        if (need_cmp) {
+            double *gsl_x = GSL_gauss(n, a, b);
+            // Сравнение векторов
+            for (int i = 0; i < n; i++) {
+                if (fabs(x[i] - gsl_x[i]) > 0.0001) {
+                    fprintf(stderr, "Invalid result: elem %d: %f %f\n", i, x[i], gsl_x[i]);
+                    break;
+                }
+            }
+            free(gsl_x);
+        }
+        free(b);
+        free(a);
+        free(x);
+    }
+    printf("Gaussian Elimination (%s): n %d, mean-time (sec) %.6f\n", algo_name, n, ((double) clock_sum) / (double)exp_n);
+    return 0;
+}

start.sh

@@ -0,0 +1,16 @@
+#!/bin/bash
+
+touch gauss_results.txt
+> gauss_results.txt
+rm gauss
+gcc main.c gauss_algos.c -o gauss -lgsl -lgslcblas -fopenmp -O0
+exp_N=10
+for size_N in {10,100,200,500,1000,2000,4000}
+do
+  {
+    ./gauss "$size_N" "$exp_N" "sequential"
+    ./gauss "$size_N" "$exp_N" "gsl"
+    ./gauss "$size_N" "$exp_N" "parallel"
+    echo ' '
+  } >> gauss_results.txt
+done