diff --git a/src/gauss_seidel.c b/src/gauss_seidel.c
new file mode 100644
index 0000000..aa9a8de
--- /dev/null
+++ b/src/gauss_seidel.c
@@ -0,0 +1,88 @@
+/*
+ * Gauss-Seidel Descendant PageRank Implementation
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "sparse_matrix.h"
+#include "vector.h"
+
+
+double* gauss_seidel_pagerank(const SparseMatrix *matrix, double epsilon, double alpha) {
+    int N = matrix->num_nodes;
+    size_t vec_size = N * sizeof(double);
+
+    double *x = malloc(vec_size);
+    double *x_old = malloc(vec_size);
+    double *f = malloc(vec_size);
+    double right_const = (1.0 - alpha) / N;
+
+    if (!x || !x_old || !f) {
+        fprintf(stderr, "Memory allocation failed\n");
+        exit(EXIT_FAILURE);
+    }
+
+    init_vector(x, N, 1.0 / N);
+    generate_f(matrix, f);
+
+    double diff;
+    int iter = 0;
+
+    do {
+        // 1. store x from previous step
+        memcpy(x_old, x, vec_size);
+
+        // 2. alpha/N * (x * f)
+        double x_f = vec_product(x_old, f, N);
+        double right_var = (alpha / (double)N) * x_f;
+
+        // 3. alpha*(pi*M) + (right_const+alpha/N * (pi * f))*e
+        //    splitted in 2 parts (2 sums instead of one)
+        //    Can't use previously defined multiply_vector_matrix_parallel function here
+        //
+        //    sum from 1 to i-1: use old x
+        //    i=j: use value M(i,i)
+        //    sum from i+1 to n: use x
+        //
+        //    Gauss Seidel descending, so i from n to 0
+        //
+        for (int i = N - 1; i >= 0; i--) {
+            double result_mult = 0.0;
+            double M_ii = 0.0;
+            for (int k = matrix->row_ptr[i]; k < matrix->row_ptr[i + 1]; k++) {
+                int j = matrix->arcs[k].dest;
+                double value = matrix->arcs[k].value;
+
+                if (j == i) {
+                    M_ii = value;
+                } else if (j > i) {
+                    result_mult += value * x[j];
+                } else {
+                    result_mult += value * x_old[j];
+                }
+            }
+            double diag_dependant_term = 1.0 - alpha * M_ii;
+
+            x[i] = (alpha * result_mult + right_const + right_var) / diag_dependant_term;
+        }
+
+        // 4. renormalization step
+        normalize_vector(x, N);
+
+        // 5. convergence
+        diff = diff_norm_vector(x, x_old, N);
+
+        // if ((++iter) % 10 == 0) {
+        //     printf("Gauss-Seidel Iteration %d: diff = %.16f\n", iter, diff);
+        // }
+        ++iter;
+    } while (diff > epsilon);
+
+    printf("Gauss-Seidel: %d Iterations\n", iter);
+
+    free(x_old);
+    free(f);
+
+    // x ~ pi (convergence)
+    return x;
+}
diff --git a/src/main.c b/src/main.c
index 23f239f..a6d0999 100644
--- a/src/main.c
+++ b/src/main.c
@@ -8,8 +8,9 @@
 #include "gauss_seidel.h"
 #include "time_helper.h"
 #include <sys/time.h>
+#include <string.h>
 
-void test_stationary_distribution(const char *path) {
+void test_stationary_distribution(const char *path, double epsilon, double alpha) {
     struct timeval tvstart, tvend;
 
     gettimeofday(&tvstart, NULL);
@@ -18,9 +19,6 @@ void test_stationary_distribution(const char *path) {
     gettimeofday(&tvend, NULL);
     print_time_diff("Read and convert matrix", &tvstart, &tvend);
 
-    double epsilon = 1e-10;
-    double alpha = 0.8;
-
     gettimeofday(&tvstart, NULL);
     double *pi_power = pagerank(matrix, epsilon, alpha);
     gettimeofday(&tvend, NULL);
@@ -30,7 +28,7 @@ void test_stationary_distribution(const char *path) {
     double *pi_new = gauss_seidel_pagerank(matrix, epsilon, alpha);
 
     gettimeofday(&tvend, NULL);
-    print_time_diff("Gauss-Seidel:", &tvstart, &tvend);
+    print_time_diff("Gauss-Seidel", &tvstart, &tvend);
 
     double diff = diff_norm_vector(pi_power, pi_new, matrix->num_nodes);
     printf("Difference norm Pagerank and Gauss-Seidel: %.16f\n", diff);
@@ -40,7 +38,29 @@ void test_stationary_distribution(const char *path) {
     free_sparse_matrix(matrix);
 }
 
-int main() {
-    test_stationary_distribution("./out/input.mtx");
+int main(int argc, char *argv[]) {
+    const char *matrix_path = NULL;
+    double epsilon = -1.0;
+    double alpha = -1.0;
+
+    for (int i = 1; i < argc; i++) {
+        if ((strcmp(argv[i], "--matrix") == 0 || strcmp(argv[i], "-matrix") == 0) && i + 1 < argc) {
+            matrix_path = argv[++i];
+        } else if ((strcmp(argv[i], "--epsilon") == 0 || strcmp(argv[i], "-epsilon") == 0) && i + 1 < argc) {
+            epsilon = atof(argv[++i]);
+        } else if ((strcmp(argv[i], "--alpha") == 0 || strcmp(argv[i], "-alpha") == 0) && i + 1 < argc) {
+            alpha = atof(argv[++i]);
+        } else {
+            fprintf(stderr, "Unknown or incomplete argument: %s\n", argv[i]);
+            return EXIT_FAILURE;
+        }
+    }
+
+    if (matrix_path == NULL || epsilon <= 0 || alpha < 0 || alpha > 1) {
+        fprintf(stderr, "Usage: %s --matrix <path> --epsilon <value> --alpha <value>\n", argv[0]);
+        return EXIT_FAILURE;
+    }
+
+    test_stationary_distribution(matrix_path, epsilon, alpha);
     return 0;
 }