gauss seidel descending

Makefile

@@ -3,6 +3,7 @@
 # --------------------------------------------------
 CC       := gcc
 CFLAGS   := -Wall -fopenmp -O3
+LDFLAGS  := -lm
 SRCDIR   := src
 OBJDIR   := out
 DATAPATH := data/web-Google/web-Google.mtx
@@ -27,9 +28,9 @@ sparse: $(OBJDIR)/sparse | $(OBJDIR)
 # --------------------------------------------------
 $(OBJDIR)/sparse: $(OBJS) $(DATAPATH) | $(OBJDIR)
 	@echo "→ Copying input data"
-	cp $(DATAPATH) $(OBJDIR)/input.rb
+	cp $(DATAPATH) $(OBJDIR)/input.mtx
 	@echo "→ Linking $@"
-	$(CC) $(CFLAGS) -o $@ $(OBJS)
+	$(CC) $(CFLAGS) $(LDFLAGS) -o $@ $(OBJS)
 
 # --------------------------------------------------
 # Compile

src/main.c

@@ -2,33 +2,58 @@
 #include <stdlib.h>
 #include "matrix_operation.h"
 #include "sparse_matrix.h"
-#include "read_from_rb.h"
+#include "read_from_mtx.h"
+#include "vector.h"
 #include "power_algorithm.h"
+#include "gauss_seidel.h"
 #include "time_helper.h"
 #include <sys/time.h>
 
-void test_pagerank(const char *path) {
-    struct timeval tvstart, tv;
-    gettimeofday(&tvstart, NULL);
+void test_stationary_distribution(const char *path) {
+    struct timeval tvstart, tvend;
 
+    // Read and prepare the matrix
+    gettimeofday(&tvstart, NULL);
     SparseMatrix *matrix = read_sparse_matrix_from_mtx(path);
     convert_to_stochastic(matrix);
-    // Time 2
-    gettimeofday(&tv, NULL);
-    print_time_diff("read matrix", &tvstart, &tv);
-    double *result = malloc(matrix->num_nodes * sizeof(double));
-    if (result == NULL) {
-        fprintf(stderr, "Memory allocation failed\n");
-        exit(EXIT_FAILURE);
+    gettimeofday(&tvend, NULL);
+    print_time_diff("Read and convert matrix", &tvstart, &tvend);
+
+    double epsilon = 1e-10;
+    double alpha = 0.8;
+
+    // pagerank
+    gettimeofday(&tvstart, NULL);
+    double *pi_power = pagerank(matrix, epsilon, alpha);
+    gettimeofday(&tvend, NULL);
+    print_time_diff("Pagerank", &tvstart, &tvend);
+
+    // nable delta
+    gettimeofday(&tvstart, NULL);
+    double *pi_new = gauss_seidel_pagerank(matrix, epsilon, alpha);
+
+    gettimeofday(&tvend, NULL);
+    print_time_diff("Gauss Seidel algorithm", &tvstart, &tvend);
+
+    // Compare results
+    double diff = diff_norm_vector(pi_power, pi_new, matrix->num_nodes);
+    printf("Difference between power method and new algorithm: %.16f\n", diff);
+
+
+    for (int i = 0; i < 10; ++i) {
+        printf("pi_old: %.16f\n", pi_power[i]);
+        printf("pi_new: %.16f\n", pi_new[i]);
+        printf("\n");
     }
-    result = pagerank(matrix, 1e-10, 0.35);
-    // Time 3
-    gettimeofday(&tv, NULL);
-    print_time_diff("finish", &tvstart, &tv);
+
+
+    // Clean up
+    free(pi_power);
+    free(pi_new);
+    free_sparse_matrix(matrix);
 }
 
 int main() {
-    test_pagerank("./out/input.rb");
-
+    test_stationary_distribution("./out/input.mtx");
     return 0;
 }

src/matrix_operation.c

@@ -66,3 +66,9 @@ void convert_to_stochastic(SparseMatrix *matrix) {
     free(non_zero);
 }
 
+int compare_arcs(const void *a, const void *b) {
+    Arc *arc_a = (Arc *)a;
+    Arc *arc_b = (Arc *)b;
+    return arc_a->origin - arc_b->origin;
+}
+

src/matrix_operation.h

@@ -5,5 +5,6 @@
 void multiply_vector_matrix(const double *vector, const SparseMatrix *matrix, double *result);
 void multiply_vector_matrix_parallel(const double *vector, const SparseMatrix *matrix, double *result);
 void convert_to_stochastic(SparseMatrix *matrix);
+int compare_arcs(const void *a, const void *b);
 
 #endif

src/read_from_mtx.c

@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include "sparse_matrix.h"
+#include "matrix_operation.h"
 
 int read_dims_ignore_comment(SparseMatrix *matrix, FILE *file) {
     char buffer[1024];
@@ -66,6 +67,22 @@ SparseMatrix* read_sparse_matrix_from_mtx(const char *filename) {
 
     parse_arcs(matrix, file);
 
+    // Sort arcs by origin
+    qsort(matrix->arcs, matrix->num_arcs, sizeof(Arc), compare_arcs);
+    // Build row_ptr
+    matrix->row_ptr = (int *)malloc((matrix->num_nodes + 1) * sizeof(int));
+    int current_origin = 0;
+    matrix->row_ptr[0] = 0;
+    for (int k = 0; k < matrix->num_arcs; k++) {
+        while (matrix->arcs[k].origin > current_origin) {
+            current_origin++;
+            matrix->row_ptr[current_origin] = k;
+        }
+    }
+    for (int i = current_origin + 1; i <= matrix->num_nodes; i++) {
+        matrix->row_ptr[i] = matrix->num_arcs;
+    }
+
     fclose(file);
     return matrix;
 }

src/read_from_rb.c

@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include "read_from_rb.h"
 #include "sparse_matrix.h"
+#include "matrix_operation.h"
 
 SparseMatrix* read_sparse_matrix_from_rb(const char *filename) {
     FILE *file = fopen(filename, "r");
@@ -97,6 +98,24 @@ SparseMatrix* read_sparse_matrix_from_rb(const char *filename) {
         }
     }
 
+
+    // Sort arcs by origin
+    qsort(matrix->arcs, matrix->num_arcs, sizeof(Arc), compare_arcs);
+
+    // Build row_ptr
+    matrix->row_ptr = (int *)malloc((matrix->num_nodes + 1) * sizeof(int));
+    int current_origin = 0;
+    matrix->row_ptr[0] = 0;
+    for (int k = 0; k < matrix->num_arcs; k++) {
+        while (matrix->arcs[k].origin > current_origin) {
+            current_origin++;
+            matrix->row_ptr[current_origin] = k;
+        }
+    }
+    for (int i = current_origin + 1; i <= matrix->num_nodes; i++) {
+        matrix->row_ptr[i] = matrix->num_arcs;
+    }
+
     free(col_ptr);
     free(row_ind);
     free(values);

src/sparse_matrix.c

@@ -4,6 +4,7 @@
 void free_sparse_matrix(SparseMatrix *matrix) {
     if (matrix) {
         free(matrix->arcs);
+        free(matrix->row_ptr);
         free(matrix);
     }
 }

src/sparse_matrix.h

@@ -11,6 +11,7 @@ typedef struct {
     int num_nodes;
     int num_arcs;
     Arc *arcs;
+    int *row_ptr;
 } SparseMatrix;
 
 void free_sparse_matrix(SparseMatrix *matrix);