use dsyevd to diagonalize

diagonalize.cpp

/**
 * Helper for exact diagonalization: diagonalize the matrix given in mat with dimension dim, stores
 * the eigenvalues in eigs and optionally the eigenvectors in mat
 * @param dim the dimension of the matrix
 * @param mat the actual matrix (size: dim*dim)
 * @param eigs array of size dim to store eigenvalues
 * @param calc_eigenvectors whether or not the calculate the eigenvectors
 */
void Diagonalize(int dim, double *mat, double *eigs, bool calc_eigenvectors)
{
    assert(mat && "mat must be allocated");
    assert(eigs && "eigs must be allocated");

    char jobz;

    if(calc_eigenvectors)
        jobz = 'V';
    else
        jobz = 'N';

    char uplo = 'U';

    int lwork, liwork;

    if(calc_eigenvectors)
    {
        lwork = 6*dim+1+2*dim*dim;
        liwork = 3+5*dim;
    } else
    {
        lwork = 2*dim+1;
        liwork = 1;
    }

    double *work = new double[lwork];

    int *iwork = new int[liwork];

    // initialized to avoid uninitialized value errors
    int info = 1;

    dsyevd_(&jobz, &uplo, &dim, mat, &dim, eigs, work, &lwork, iwork, &liwork, &info);

    if(info != 0)
        std::cerr << "Calculating eigenvalues failed..." << std::endl;

    delete [] work;
    delete [] iwork;
}