I have written a code which should be doing the same in both languages.
It generates a matrix of a given size, populates it with random numbers, then makes it hermitian, and last it calculates its eigenvalues and eigenvectors.

I have ran the code for different matrix sizes both in python and c++, and got following results:

I don’t understand why from around matrix size = 200 the c++ code starts to get slower and slower.

Here is my c++ code:

#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <random>
#include <ctime>

using namespace std;

int main(int argc, char* argv[]) {
    if (argc != 2) {
        cout << "Usage: " << argv[0] << " <matrix_size>" << endl;
        return 1;
    }

    const size_t matrix_size = atoi(argv[1]);

    gsl_matrix* matrix = gsl_matrix_calloc(matrix_size, matrix_size);
    gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(matrix_size);
    gsl_vector* eigenvalues = gsl_vector_alloc(matrix_size);
    gsl_matrix* eigenvectors = gsl_matrix_alloc(matrix_size, matrix_size);
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<> dis(0.0, 1.0);

    // Initialize matrix with random numbers
    for (size_t i = 0; i < matrix_size; ++i) {
        for (size_t j = 0; j < matrix_size; ++j) {
            gsl_matrix_set(matrix, i, j, dis(gen));
        }
    }

    // Make the matrix symmetric
    for (size_t i = 0; i < matrix_size; ++i) {
        for (size_t j = i + 1; j < matrix_size; ++j) {
            double value = gsl_matrix_get(matrix, i, j);
            gsl_matrix_set(matrix, j, i, value);
        }
    }

    // Calculate eigenvalues and eigenvectors
    clock_t start_time = clock();
    gsl_eigen_symmv(matrix, eigenvalues, eigenvectors, workspace);
    clock_t end_time = clock();
    double elapsed_time = double(end_time - start_time) / CLOCKS_PER_SEC;

    // Print eigenvalues and eigenvectors
    cout << elapsed_time << endl;

    // Free allocated memory
    gsl_matrix_free(matrix);
    gsl_eigen_symmv_free(workspace);
    gsl_vector_free(eigenvalues);
    gsl_matrix_free(eigenvectors);

    return 0;
}

compiled with g++ -O3 gsltests.cpp -o gsltests -lgsl -lgslcblas.

And my python code:

import numpy as np
from scipy.linalg import eigh
import time


def main(matrix_size):
    matrix = np.random.rand(matrix_size, matrix_size)
    matrix = (matrix + matrix.T) / 2  # Make the matrix symmetric

    start_time = time.time()
    eigenvalues, eigenvectors = eigh(matrix)
    end_time = time.time()
    elapsed_time = end_time - start_time

    print(elapsed_time)

if __name__ == "__main__":
    import sys
    if len(sys.argv) != 2:
        print("Usage:", sys.argv[0], "<matrix_size>")
        sys.exit(1)
    main(int(sys.argv[1]))