Spectra, Eigen, EigenSolver C++. Not giving smallest Eigenvalue

Question

It gives largest eigenvalue

Does not give smallest one

I'll write the post again, copying all the code this time.You can run it in less than 3 seconds. about 10 to compile.

As you can see, it does return the largest value, but not the smallest.

I printed a small matrix just to show it to you. It seems there are no issue with a 2^3 x 2^3 ones though. The example in the documentation works just fine.

I'll write the code:

#include <iostream>
#include <random>
#include <Eigen/Sparse>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
using namespace std;
using uint = unsigned int;

#define ALLOCATE_ARRAY(type, count) ((type*) malloc(sizeof(type) * (count)))
#define REALLOC_ARRAY(array, type, count) ((type*) realloc(array, sizeof(type) * count))
#define MATRIX(array, i,j) array[j+i*total_sites]


uint number_sites = 4;
uint total_sites = 0;
uint total_Hilbert_elements = 0;
uint number_spin_states = 2;
int* computational_base = 0;

default_random_engine ran_generator;

void inizialize_Hamiltonian_matrix (float trasverse_magnetic_field, float longitudinal_magnetic_field, uint boundary_conditions_flag, Eigen::MatrixXf& Hamiltonian_config_matrix);

int main(){

    float beta = 0.3;
    float trasverse_magnetic_field = 1;
    float longitudinal_magnetic_field = 1;

    uint dimensions = 1;
    uint initial_state = 1;
    uint boundary_conditions_flag = 1;
    uint number_eigenvalues = 1;
    total_sites = pow(number_sites, dimensions);

    total_Hilbert_elements = pow(number_spin_states, total_sites);
    computational_base =  ALLOCATE_ARRAY(int, total_Hilbert_elements*total_sites);

    for (int i=0; i<total_Hilbert_elements; ++i){
        uint i_holder = i;
        for (int j=0; j<total_sites; ++j){
            MATRIX(computational_base, i,j) = i_holder%2;
            i_holder = i_holder/2;
        }
    }

    // Build Matrix and print on shell
    Eigen::MatrixXf Hamiltonian_config_matrix(total_Hilbert_elements, total_Hilbert_elements);
    inizialize_Hamiltonian_matrix (trasverse_magnetic_field, longitudinal_magnetic_field, boundary_conditions_flag, Hamiltonian_config_matrix);

    for (int i=0; i<total_Hilbert_elements; ++i){
        for (int j=0; j<total_Hilbert_elements; ++j){
            if (Hamiltonian_config_matrix(i,j)>=0){
                cout << " " << Hamiltonian_config_matrix(i,j);
            }
            else{
                cout << Hamiltonian_config_matrix(i,j);
            }
        }
        cout << endl;
    }

    // Converting in Sparse Matrix and Spectra stuff
    Eigen::SparseMatrix<float> Sparse_Hamiltonian_config;
    Sparse_Hamiltonian_config = Hamiltonian_config_matrix.sparseView();
    Spectra::SparseGenMatProd<float> Operation_object_wrap(Sparse_Hamiltonian_config);

    // Initialize and Compute
    Spectra::GenEigsSolver<Spectra::SparseGenMatProd<float>> Ready_for_Lanczos(Operation_object_wrap, number_eigenvalues, 2*number_eigenvalues+1);
    Ready_for_Lanczos.init();

    uint which_eigenvalues = 0;
    cout << "Type 1 if Largest Eigenvalue, Type 0 if Smallest" << endl;
    cin >> which_eigenvalues;

    if (which_eigenvalues==0){
        int nconv = Ready_for_Lanczos.compute(Spectra::SortRule::SmallestMagn);
    }
    else{
        int nconv = Ready_for_Lanczos.compute(Spectra::SortRule::LargestMagn);
    }

    // Retrieve results
    Eigen::VectorXcf evalues;
    if(Ready_for_Lanczos.info() == Spectra::CompInfo::Successful)
        evalues = Ready_for_Lanczos.eigenvalues();
    cout << "Eigenvalues found:\n" << evalues << endl;




    return 0;
}



void inizialize_Hamiltonian_matrix (float trasverse_magnetic_field, float longitudinal_magnetic_field, uint boundary_conditions_flag, Eigen::MatrixXf& Hamiltonian_config_matrix){
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==MATRIX(computational_base, i,j+1)){
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -1.0;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +1.0;
            }

        }
    }
    if (boundary_conditions_flag==1){
        uint j = total_sites-1;
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==MATRIX(computational_base, i,0)){
                    Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -1.0;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +1.0;
            }
        }
    }
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==1){
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) -longitudinal_magnetic_field;
            }
            else{
                Hamiltonian_config_matrix(i,i) = Hamiltonian_config_matrix(i,i) +longitudinal_magnetic_field;
            }

        }
    }
    uint new_i_binary = 0;
    for (uint j=0; j<total_sites-1; ++j){
        for (uint i=0; i<total_Hilbert_elements; ++i){
            if (MATRIX(computational_base, i,j)==1){
                new_i_binary = i - pow(2, j);
            }
            else{
                new_i_binary = i + pow(2, j);
            }
            Hamiltonian_config_matrix(new_i_binary,i) = Hamiltonian_config_matrix(new_i_binary,i) -trasverse_magnetic_field;
        }
    }
}

If you can show me what i did wrong, it could be a huge help. Thanks