void fcn(int n, d_complex x[], int itask, d_complex y[])
{
    switch (itask) {
        /*
         * Define responses to different tasks for the generalized
         * eigenvalue problem
         *        A*x = lambda * B * x,
         * which includes the ordinary case B = I.
         */
    case IMSL_ARPACK_PREPARE:
        /* 
         * Take initial steps to prepare for the operations
         * that follow. Note that arpack_general internally
         * always calls fcn with this enum value, even if it is
         * not required by the user. 
         */
        break;
    case IMSL_ARPACK_A_X:
        /*
         * Compute matrix-vector product y = A * x
         */
        break;
    case IMSL_ARPACK_B_X:
        /*
         * Compute matrix-vector product y = B * x
         */
        break;
    case IMSL_ARPACK_INV_SHIFT_X:
        /*
         *  Compute the matrix-vector product
         *     z = inv(A - sigma * B) * x,
         *  and return
         *     y = z.
         *
         *  Usually, matrix A - sigma * B is not directly inverted.
         *  Instead, a factorization of A - sigma * B is determined,
         *  and the factors are used to compute z via backsolves.
         *
         *  Example:
         *  If an LU factorization of A - sigma * B exists, then
         *      A - sigma * B = P * L * U,
         *  P a permutation matrix. Vector z can then be determined
         *  as solution of the linear system
         *      L * U * z = trans(P) * x.
         *  The LU factorization only has to be computed once, for
         *  example outside of fcn or within IMSL_ARPACK_PREPARE.
         */
        break;
    case IMSL_ARPACK_INV_B_X:
        /*
         *  Compute matrix-vector product
         *     y = inv(B) * x.
         *  Usually, matrix B is not directly inverted.
         *  Instead, a factorization of B is determined, and the
         *  factors are used to compute y via backsolves.
         *  
         *  Example:
         *  If matrix B is positive definite, then a Cholesky
         *  factorization B = L * ctrans(L) exists. Vector y can then
         *  be determined by solving the linear system
         *     L * ctrans(L) * y = x
         *  via the two backsolves
         *     L * z = x  and  ctrans(L) * y = z.
         *  The Cholesky factorization only has to be computed once,
         *  for example outside of fcn or within IMSL_ARPACK_PREPARE.
         */
        break;
    default:
        /*
         *  Define error conditions, if necessary.
         */
        break;
    }
}

#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>

static void fcn(int, d_complex[], int, d_complex[]);

int main()
{
    int n = 4, nev = 2, n_acc;
    d_complex *eval = NULL;

    /* Compute eigenvalues */
    eval = imsl_z_arpack_general(fcn, n, nev,
        IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
        0);

    printf("Number of requested eigenvalues : %d\n", nev);
    printf("Number of computed (converged) eigenvalues : %d\n", n_acc);
    
    /* Print eigenvalues of largest magnitude */
    imsl_z_write_matrix("Largest eigenvalues", 1, n_acc, eval, 0);

    if (eval)
        imsl_free(eval);
}

static void fcn(int n, d_complex x[], int itask, d_complex y[])
{
    d_complex a[] = {
        {5.0, 9.0}, {5.0, 5.0}, {-6.0, -6.0}, {-7.0, -7.0},
        {3.0, 3.0}, {6.0, 10.0}, {-5.0, -5.0}, {-6.0, -6.0},
        {2.0, 2.0}, {3.0, 3.0}, {-1.0, 3.0}, {-5.0, -5.0},
        {1.0, 1.0}, {2.0, 2.0}, {-3.0, -3.0}, {0.0, 4.0}
    };

    switch (itask) {
    case IMSL_ARPACK_PREPARE:
        /* Do nothing here */
        break;
    case IMSL_ARPACK_A_X:
        imsl_z_mat_mul_rect("A*x",
            IMSL_A_MATRIX, n, n, a,
            IMSL_X_VECTOR, n, x,
            IMSL_RETURN_USER, y,
            0);
        break;
    default:
        imsl_set_user_fcn_return_flag(1);
        break;
    }
}

Number of requested eigenvalues : 2
Number of accurate (converged) eigenvalues : 2
 
                 Largest eigenvalues
                        1                          2
(          4,          8)  (          3,          7)

#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>

static void fcn(int, d_complex[], int, d_complex[]);
static void fcn_w_data(int, d_complex[], int, d_complex[], void *);
static void ax(int, double, d_complex[], d_complex[], d_complex[]);
static void tx(int, d_complex[], d_complex[], d_complex[]);
static void compute_residuals_fcn(
    void(*fcn_w_data) (int n, d_complex x[], int itask, d_complex y[],
        void *data),
    int n, int n_acc, d_complex eval[], d_complex evecu[],
    int evecu_col_dim, d_complex cwork[], double res_norm[], void *data);
static d_complex l_mul(d_complex, d_complex);
static void l_zaxpy(int, d_complex, d_complex[], d_complex[]);

typedef struct {
    int nx;
    double hsq;
    d_complex rho;
    d_complex *dt;
} imsl_arpack_data;

int main() {
    int nx = 10, nev = 6, ncv = 25, n_acc, n, i;
    d_complex *eval = NULL, dt[3];
    d_complex *evecu = NULL, *cwork = NULL;
    d_complex rho = { 100.0, 1.0 };
    double *res_norm = NULL;
    imsl_arpack_data usr_data;

    usr_data.nx = nx;
    usr_data.rho = rho;
    usr_data.dt = dt;

    n = nx * nx;

    evecu = (d_complex *)malloc(n * nev * sizeof(d_complex));
    cwork = (d_complex *)malloc(2 * n * sizeof(d_complex));
    res_norm = (double *)malloc(nev * sizeof(double));

    /* Compute eigenvalues */
    eval = imsl_z_arpack_general(NULL, n, nev,
        IMSL_NUM_ACCURATE_EIGVALS, &n_acc,
        IMSL_NUM_ARNOLDI_VECTORS, ncv,
        IMSL_VECTORS_USER, evecu,
        IMSL_FCN_W_DATA, fcn_w_data, &usr_data,
        0);

    compute_residuals_fcn(fcn_w_data, n, n_acc, eval, evecu, nev, cwork,
        res_norm, &usr_data);

    printf("\nNumber of requested eigenvalues : %d\n", nev);
    printf("Number of computed (converged) eigenvalues : %d\n",
        n_acc);
    printf("\n%20sLargest Magnitude Operator eigenvalues\n", " ");
    printf("%15sRitz values (Real, Imag)%18srelative "
        "residuals\n", " ", " ");
    for (i = 0; i < n_acc; i++) {
        printf("%d  (%20.13lf, %20.13lf)    %25.17lf\n",
            i, eval[i].re, eval[i].im, res_norm[i]);
    }

    if (eval)
        imsl_free(eval);
    if (evecu)
        free(evecu);
    if (cwork)
        free(cwork);
    if (res_norm)
        free(res_norm);
}

static void fcn_w_data(int n, d_complex x[], int itask,
    d_complex y[], void *data)
{
    int nx;
    double h, hsq;
    d_complex dd, dl, du, rho;
    d_complex *dt;
    imsl_arpack_data *usr_data = (imsl_arpack_data *)data;
    dt = usr_data->dt;
    rho = usr_data->rho;
    nx = usr_data->nx;

    switch (itask) {
    case IMSL_ARPACK_PREPARE:
        /* Create tridiagonal matrix in band storage format */
        h = 1.0 / (nx + 1);
        hsq = h * h;
        dd.re = 4.0 / hsq;
        dd.im = 0.0;
        dl.re = -1.0 / hsq - 0.5 * rho.re / h;
        dl.im = -0.5 * rho.im / h;
        du.re = -1.0 / hsq + 0.5 * rho.re / h;
        du.im = 0.5 * rho.im / h;

        dt[0] = dl;
        dt[1] = dd;
        dt[2] = du;
        usr_data->hsq = hsq;
        break;
    case IMSL_ARPACK_A_X:
        hsq = usr_data->hsq;
        ax(nx, hsq, dt, x, y);
        break;
    default:
        imsl_set_user_fcn_return_flag(1);
        break;
    }
}

/*
 *     Matrix-vector function.
 *
 *     Computes y <- A*x, where A is the nx^2 by nx^2 block
 *     tridiagonal matrix
 *
 *                 | T -I          |
 *                 |-I  T -I       |
 *             A = |   -I  T       |
 *                 |        ...  -I|
 *                 |           -I T|
 *
 *    derived from the standard central difference discretization
 *    of the convection-diffusion operator (Laplacian u) + rho*(du/dx)
 *    with zero boundary condition.
 *
 *    Function tx() is called to compute y <- T*x.
 */
static void ax(int nx, double hsq, d_complex dt[], d_complex x[],
    d_complex y[]) {
    int i, j;

    tx(nx, x, dt, y);
    for (i = 1; i <= nx; i++) {
        y[i - 1].re -= x[nx + i - 1].re / hsq;
        y[i - 1].im -= x[nx + i - 1].im / hsq;
    }

    for (j = nx + 1; j <= nx * (nx - 1); j += nx) {
        tx(nx, &x[j - 1], dt, &y[j - 1]);
        for (i = j; i <= j + nx - 1; i++) {
            y[i - 1].re -= (x[i - nx - 1].re + x[i + nx - 1].re) / hsq;
            y[i - 1].im -= (x[i - nx - 1].im + x[i + nx - 1].im) / hsq;
        }

    }

    tx(nx, &x[nx * (nx - 1)], dt, &y[nx * (nx - 1)]);
    for (i = (nx - 1) * nx + 1; i <= nx * nx; i++) {
        y[i - 1].re -= x[i - nx - 1].re / hsq;
        y[i - 1].im -= x[i - nx - 1].im / hsq;
    }
}

/*
 *    Compute the matrix vector multiplication y <- T * x
 *    where T is an nx by nx tridiagonal matrix with dd on the
 *    diagonal, dl on the subdiagonal, and du on the superdiagonal.
 */
static void tx(int nx, d_complex x[], d_complex dt[],
    d_complex v[]) {
    int j;
    d_complex dl, dd, du, temp1, temp2, temp3;

    dl = dt[0];
    dd = dt[1];
    du = dt[2];

    temp1 = l_mul(dd, x[0]);
    temp2 = l_mul(du, x[1]);
    v[0].re = temp1.re + temp2.re;
    v[0].im = temp1.im + temp2.im;

    for (j = 2; j <= nx - 1; j++) {
        temp1 = l_mul(dl, x[j - 2]);
        temp2 = l_mul(dd, x[j - 1]);
        temp3 = l_mul(du, x[j]);

        v[j - 1].re = temp1.re + temp2.re + temp3.re;
        v[j - 1].im = temp1.im + temp2.im + temp3.im;
    }

    temp1 = l_mul(dl, x[nx - 2]);
    temp2 = l_mul(dd, x[nx - 1]);
    v[nx - 1].re = temp1.re + temp2.re;
    v[nx - 1].im = temp1.im + temp2.im;
}

/*
 *    Compute the residual norm
 *
 *     || A * x - lambda * x || / |lambda|
 *
 *    for the n_acc accurately computed eigenvalues and
 *    eigenvectors.
 */
static void compute_residuals_fcn(
    void(*fcn_w_data) (int n, d_complex x[], int itask, d_complex y[],
        void *data),
    int n, int n_acc, d_complex eval[], d_complex evecu[],
    int evecu_col_dim, d_complex cwork[], double res_norm[], void *data) {
    int i, j;
    double eval_abs;
    d_complex temp;

    for (j = 0; j < n_acc; j++) {
        /* Copy j-th eigenvector into cwork */
        for (i = 0; i < n; i++) {
            cwork[i] = evecu[i * evecu_col_dim + j];
        }
        /*  Compute ax <- A * x  */
        fcn_w_data(n, cwork, IMSL_ARPACK_A_X, &cwork[n], data);
        /*  Compute ax <- ax - lambda * x */
        temp.re = -eval[j].re;
        temp.im = -eval[j].im;
        l_zaxpy(n, temp, cwork, &cwork[n]);
        /* Save residual norms */
        res_norm[j] = imsl_z_vector_norm(n, &cwork[n], 0);
        eval_abs = imsl_z_abs(eval[j]);
        if (eval_abs > 0.0) {
            res_norm[j] /= eval_abs;
        }
    }
}

/*
 *  Compute complex scalar product x * y.
 *
 *  This is faster than imsl_z_mul since the IMSL error handler is not
 *  involved.
 */
static d_complex l_mul(d_complex x, d_complex y) {
    d_complex t;

    t.re = x.re * y.re - x.im * y.im;
    t.im = x.re * y.im + x.im * y.re;

    return t;
}

/*
 *  Compute y <- alpha * x + y.
 */
static void l_zaxpy(int n, d_complex alpha, d_complex x[],
    d_complex y[]) {
    int i;
    d_complex temp;

    for (i = 0; i < n; i++) {
        temp = l_mul(alpha, x[i]);
        y[i].re += temp.re;
        y[i].im += temp.im;
    }
}

Number of requested eigenvalues : 6
Number of computed (converged) eigenvalues : 6

                    Largest Magnitude Operator eigenvalues
               Ritz values (Real, Imag)             relative residuals
0  (   727.0167725961721,  -1029.5865512774324)     0.00000000000000454
1  (   705.3778266332353,   1029.5865512774278)     0.00000000000000788
2  (   698.4028279266130,  -1029.5865512774401)     0.00000000000000541
3  (   676.7638819636808,   1029.5865512774371)     0.00000000000000512
4  (   653.2957705962293,  -1029.5865512774360)     0.00000000000000389
5  (   631.6568246332987,   1029.5865512774408)     0.00000000000000820