Test of different LAPACK functions for computing eigenvalues of a symmetric matrix (corresponding to the routines used by numpy.linalg.eigh and scipy.linalg.eigh, and numpy.linalg.eig)

testcase.cc

// Compile with: `clang++ -o test testcase.cc -framework Accelerate`
// run with ./test <N> <seed>

// Copyright [2015] Robert T. McGibbon
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
// The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#include "Accelerate/Accelerate.h"
#include <cstdlib>
#include <cstdio>
#include <algorithm>  // std::sort
#include <vector>

using std::vector;

extern "C" int
dsyevd(char *jobz, char *uplo, int *n,
        double a[], int *lda, double w[], double work[],
        int *lwork, int iwork[], int *liwork,
        int *info);

extern "C" int
dsyevr(char *jobz, char *range, char *uplo, int *n,
       double *a, int *lda, double *vl, double *vu,
       int *il, int *iu, double *abstol, int *m, double *w,
       double *z, int *ldz, int *isuppz, double *work,
       int *lwork, int *iwork, int *liwork, int *info);

extern "C" int
dgeev(char *jobvl, char *jobvr, int *n, double *
      a, int *lda, double *wr, double *wi, double *vl,
      int *ldvl, double *vr, int *ldvr, double *work,
      int *lwork, int *info);


typedef struct eigh_params_struct {
    char JOBZ;
    char UPLO;
    double* work;
    int lwork;
    int* iwork;
    int liwork;
    int info;

    char RANGE;
    double vl;
    double vu;
    int il;
    int iu;
    double abstol;
} EIGH_PARAMS_t;

void free_params(EIGH_PARAMS_t* params) {
    free(params->work);
    free(params->iwork);
}

void init_dsyevd(EIGH_PARAMS_t* params, char JOBZ, char UPLO, int N) {
    double* null = NULL;
    double* work = NULL;
    int lwork  = -1;
    int liwork = -1;
    int info;

    double query_work_size = 0;
    int query_iwork_size = 0;

    dsyevd(&JOBZ, &UPLO, &N,
           null, &N, null,
           &query_work_size, &lwork,
           &query_iwork_size, &liwork,
           &info);

    params->lwork = (int) query_work_size;
    params->work = (double*) malloc(params->lwork*sizeof(double));
    params->liwork = query_iwork_size;
    params->iwork = (int*) malloc(params->liwork*sizeof(double));
    params->info = info;
    params->JOBZ = JOBZ;
    params->UPLO = UPLO;
}

// Compute the eigenvalues of a symmetric matrix A
// (of shape N x N in Fortran order) using dsyevd, the LAPACK routine
// used by numpy.linalg.eigh
vector<double> dsyevd_eigenvalues(double* A, int N) {
    char JOBZ = 'V';
    char UPLO = 'L';
    EIGH_PARAMS_t params;
    init_dsyevd(&params, JOBZ, UPLO, N);

    vector<double> a(A, A+N*N);
    vector<double> w(N);

    dsyevd(&params.JOBZ, &params.UPLO, &N,
           &a[0], &N, &w[0],
           params.work, &params.lwork,
           params.iwork, &params.liwork,
           &params.info);

    free_params(&params);
    std::sort(w.begin(), w.end());
    return w;
}


void init_dsyevr(EIGH_PARAMS_t* params, char JOBZ, char UPLO, int N) {
    char RANGE = 'A';
    double vl = 0;
    double vu = 0;
    int il = 0;
    int iu = 0;
    double abstol = 0;
    double* null = NULL;

    int M;
    double query_work_size = 0;
    int query_iwork_size = 0;
    int lwork = -1;
    int liwork = -1;
    int info;
    dsyevr(&JOBZ, &RANGE, &UPLO, &N, null,
           &N, &vl, &vu, &il, &iu, &abstol,
           &M, null, null, &N, (int*) null,
           &query_work_size, &lwork,
           &query_iwork_size, &liwork,
           &info);

   params->lwork = (int) query_work_size;
   params->work = (double*) malloc(params->lwork*sizeof(double));
   params->liwork = query_iwork_size;
   params->iwork = (int*) malloc(params->liwork*sizeof(int));
   params->JOBZ = JOBZ;
   params->UPLO = UPLO;
   params->RANGE = RANGE;
   params->vl = vl;
   params->vu = vu;
   params->il = il;
   params->iu = iu;
   params->abstol = abstol;
}


// Compute the eigenvalues of a symmetric matrix A
// (of shape N x N in Fortran order) using dsyevr, the LAPACK routine
// used by scipy.linalg.eigh
vector<double> dsyevr_eigenvalues(double* A, int N) {
    char JOBZ = 'V';
    char UPLO = 'L';

    EIGH_PARAMS_t params;
    init_dsyevr(&params, JOBZ, UPLO, N);

    int M;
    vector<int> isuppz(2*N);
    vector<double> w(N);
    vector<double> z(N*N);
    vector<double> a(A, A+N*N);

    dsyevr(&params.JOBZ, &params.RANGE, &params.UPLO, &N, &a[0],
           &N, &params.vl, &params.vu, &params.il, &params.iu, &params.abstol,
           &M, &w[0], &z[0], &N, &isuppz[0],
           params.work, &params.lwork,
           params.iwork, &params.liwork,
           &params.info);

    free_params(&params);
    std::sort(w.begin(), w.end());
    return w;
}


void init_dgeev(EIGH_PARAMS_t* params, char JOBVL, char JOBVR, int N) {
    // dgeev(char *jobvl, char *jobvr, int *n, double *
    //       a, int *lda, double *wr, double *wi, double *vl,
    //       int *ldvl, double *vr, int *ldvr, double *work,
    //       int *lwork, int *info);

    double* A = NULL;
    double* WR = NULL;
    double* WI = NULL;
    double* VL = NULL;
    double* VR = NULL;
    double query_work_size;
    int lwork = -1;
    int info;

    dgeev(&JOBVL, &JOBVR, &N, A, &N, WR, WI, VL, &N, VR, &N, &query_work_size,
          &lwork, &info);

    params->lwork = (int) query_work_size;
    params->work = (double*) malloc(params->lwork*sizeof(double));
    params->liwork = 1;
    params->iwork = (int*) malloc(params->liwork*sizeof(int));
}

vector<double> dgeev_eigenvalues(double* A, int N) {
    char JOBVL = 'V';
    char JOBVR = 'V';
    EIGH_PARAMS_t params;

    init_dgeev(&params, JOBVL, JOBVR, N);

    vector<double> a(A, A+N*N);
    vector<double> wr(N);
    vector<double> wi(N);
    vector<double> vl(N*N);
    vector<double> vr(N*N);
    int info;

    dgeev(&JOBVL, &JOBVR, &N, &a[0], &N, &wr[0], &wi[0],
          &vl[0], &N, &vr[0], &N, params.work,
          &params.lwork, &info);

    free_params(&params);
    std::sort(wr.begin(), wr.end());
    return wr;
}


int main(int argc, char** argv) {
    if (argc < 3) {
        fprintf(stderr, "usage: %s N seed\n", argv[0]);
        exit(1);
    }

    int N = atoi(argv[1]);  // size of matrix
    int seed = atoi(argv[2]);
    srand48(seed);

    printf("N = %d\n", N);
    printf("seed = %d\n", seed);

    vector<double> A(N*N);
    for (int i = 0; i < N; i++) {
        for (int j = i; j < N; j++) {
            double v = drand48();
            A[i*N+j] = v;
            A[j*N+i] = v;
        }
    }

    std::vector<double> w1 = dsyevr_eigenvalues(&A[0], N);
    std::vector<double> w2 = dsyevd_eigenvalues(&A[0], N);
    std::vector<double> w3 = dgeev_eigenvalues(&A[0], N);

    printf("dsyevr Eigenvalues: [");
    for (int i = N-5; i < N; i++) {
        printf("%.8f, ", w1[i]);
    }
    printf("]\n");

    printf("dsyevd Eigenvalues: [");
    for (int i = N-5; i < N; i++) {
        printf("%.8f, ", w2[i]);
    }
    printf("]\n");

    printf("dgeev  Eigenvalues: [");
    for (int i = N-5; i < N; i++) {
        printf("%.8f, ", w3[i]);
    }
    printf("]\n");

    return 0;
}

I need a help: I am using "scipy.linalg.eig". What is the C++ Library and function name?

Please let me know: How to find the C++ library/function , which is used by scipy.linalg in Python.