Demonstraton of bug #401 in ARPACK-NG 3.9.0

arpack-ng-bug-410.c

// See:
// https://github.com/opencollab/arpack-ng/issues/410
// https://github.com/opencollab/arpack-ng/issues/401
//
// Build ARPACK with ISO_C_BINDING to make arpack.h available, then compile this file as:
// 
// cc prog.c -o prog -larpack

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include <arpack.h>

// Represents a sparse square matrix of zeros and ones.
typedef struct {
    int n; // this is an n-by-n matrix
    int nelems; // no of non-zero elements
    int *rowidx, *colidx; // vectors of row and column indices of non-zero elements
} binmatrix_t;

// Init an n-by-n binary matrix.
void binmatrix_init(binmatrix_t *mat, int n) {
    assert(n >= 0);

    mat->n = n;
    mat->nelems = 0;

    // Allocate dummy arrays of size 1 suitable for later realloc()
    // through binmatrix_enter().
    mat->rowidx = calloc(1, sizeof(int));
    mat->colidx = calloc(1, sizeof(int));
}

// Deallocate a binary matrix.
void binmatrix_destroy(binmatrix_t *mat) {
    free(mat->rowidx);
    free(mat->colidx);
}

// Add a single non-zero element, mat[row, col] = 1.
// Inefficient, but simple.
void binmatrix_enter(binmatrix_t *mat, int row, int col) {
    assert(0 <= row && row < mat->n);
    assert(0 <= col && col < mat->n);

    mat->nelems++;

    mat->rowidx = realloc(mat->rowidx, sizeof(int) * mat->nelems);
    mat->colidx = realloc(mat->colidx, sizeof(int) * mat->nelems);

    mat->rowidx[mat->nelems - 1] = row;
    mat->colidx[mat->nelems - 1] = col;
}

// Print a length-n vector of doubles.
// Provided for debugging purposes.
void print_vec(const double *vec, int n) {
    for (int i=0; i < n; i++) {
        printf("%g ", vec[i]);
    }
    printf("\n");
}

// Multiplies invec by mat and stores the result in outvec.
// Invec and outvec are assumed to be allocated and have size mat->n.
void binmatrix_mulvec(binmatrix_t *mat, const double *invec, double *outvec) {
    memset(outvec, 0, sizeof(double) * mat->n);
    for (int i=0; i < mat->nelems; i++) {
        outvec[mat->rowidx[i]] += invec[mat->colidx[i]];
    }
}

// Uniform random real in [0, 1).
double uniform_random_real() {
    return (double) rand() / (RAND_MAX + 1.0);
}

// Run the eigensolver for a specific predefined symmetric problem.
// Returns true if ARPACK indicated an error, false otherwise.
bool run_solver(void) {
    const int N = 10; // matrix size, code below assumes N >= 2
    
    // N-by-N symmetric matrix with two 1 elements
    // The bug in ARPACK 3.9.0 is reproduced with N=10 or N=20.
    // Only reproducible with certain BLAS/LAPACK implementations,
    // such as Apple's vecLib or several OpenBLAS versions.
    binmatrix_t mat;
    binmatrix_init(&mat, N);
    binmatrix_enter(&mat, 0, 1);
    binmatrix_enter(&mat, 1, 0);

    // Variables for ARPACK reverse communication interface.
    // These are the parameters of DSAUPD().
    a_int ido;
    const char *bmat = "I";
    const char *which = "LA"; // largest (algebraic) eigenvalue
    const a_int nev = 1; // find one eigenvalue
    const double tol = 0; // use default tolerance
    double *resid;
    const a_int ncv = 5; // manually set to suit N=10
    double *v;
    const a_int ldv = N;
    a_int iparam[11];
    a_int ipntr[11];
    double *workd;
    double *workl;
    a_int lworkl = ncv * (8 + ncv);
    a_int info;

    // Allocate work storage for ARPACK
    v = calloc(ncv, sizeof(double) * ldv * ncv);
    workl = calloc(lworkl, sizeof(double));
    workd = calloc(3*N, sizeof(double));
    resid = calloc(N, sizeof(double));

    // Problem parameters.
    // When referencing the docs, keep in mind that it uses Fortran-style
    // 1-based indexing, while these below are 0-based indices.
    iparam[0] = 1; // ISHIFT
    iparam[1] = 0; // not referenced
    iparam[2] = 3000; // MXITER, max no of iterations
    iparam[3] = 1; // NB, "The code currently works only for NB = 1."
    iparam[4] = 0; // NCONV, only for output
    iparam[5] = 0; // not referenced
    iparam[6] = 1; // MODE, mode 1
    iparam[7] = 0; // NP, only for output
    iparam[8] = 0; // NUMOP, only for output
    iparam[9] = 0; // NUMOPB, only for output
    iparam[10] = 0; // NUMREO, only for output

    info = 1; // we will provide our own initial resid
    for (int i=0; i < N; i++) {
        resid[i] = uniform_random_real();
    }

    ido = 0; // must be set to 0 for first DSAUPD() call
    while (true) {
        dsaupd_c(&ido, bmat, N, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);

        // On error (i.e. info != 0), ARPACK should request termination with ido == 99
        assert(ido == 99 || info == 0);

        // Handle ido values
        switch (ido) {
        case -1:
        case 1:
        {
            // perform multiplication
            double *invec  = workd + ipntr[0] - 1;
            double *outvec = workd + ipntr[1] - 1;
            binmatrix_mulvec(&mat, invec, outvec); // outvec = mat * invec
            break;
        }
        case 99:
            // finished
            goto done;
        default:
            // should never reach here
            assert(! "Unexpected IDO value!");
        }
    }
done:

    printf("Finished.\n");
    printf("INFO = %d\n", (int) info);

    // No need to do post-processing with DSEUPD() for this example

    // free storage allocated for ARPACK
    free(resid);
    free(workd);
    free(workl);
    free(v);

    binmatrix_destroy(&mat);

    return info != 0;
}

int main(void) {
    // Seed the RNG for reproducibility, but keep in mind that the
    // standard library's RNG, which we use here, is platform-specific.
    srand(42);

    // With an appropriate problem setup, the error is reproducible with high probability,
    // 100 trials are more than sufficient. Note that in each of these trials, RESID will
    // be initialized to a different random value.
    for (int i=0; i < 100; i++) {
        printf("\nTrial %d:\n", i);
        if (run_solver()) {
            // ARPACK indicated an error, stop.
            break;
        }
    }

    return 0;
}