How to properly use ARPACK's dsaupd and dseupd?

Question

In Rust, I am trying to solve an eigendecomposition problem through ARPACK. I made the following subroutine for this purpose:

fn faer_arpack_symmetric_f64<F>(
    mut transform: F,
    dimension: usize,
    eigenvalue_kind: &str,
    number_of_eigenvalues: usize,
    number_of_basis_vectors: usize,
    maxiter: usize,
    vectors: bool,
) -> (Col<f64>, Mat<f64>)
where
    F: FnMut(ColRef<f64>) -> Col<f64>,
{
    let mut ido = 0;
    let mut residual: Col<f64> = Col::zeros(dimension);
    let mut eigenvectors = Mat::<f64>::zeros(dimension, number_of_basis_vectors);
    let mut iparam = [0; 11];
    iparam[0] = 1;
    iparam[2] = maxiter as i32;
    iparam[6] = 1;
    let mut ipntr = [0; 14];
    let mut workd = Col::zeros(3 * dimension);
    let lworkl = 3 * number_of_basis_vectors.pow(2) + 6 * number_of_basis_vectors;
    let mut workl = Col::<f64>::zeros(lworkl);
    let mut info = 0;

        loop
    {
        unsafe {
            dsaupd_c(
                &mut ido,
                "I".as_ptr() as *const i8,
                dimension as i32,
                eigenvalue_kind.as_ptr() as *const i8,
                number_of_eigenvalues as i32,
                f64::EPSILON,
                residual.as_ptr_mut() as *mut f64,
                number_of_basis_vectors as i32,
                arnoldi_vectors.as_ptr_mut() as *mut f64,
                dimension as i32,
                iparam.as_mut_ptr(), // no idea what this does
                ipntr.as_mut_ptr(),  // no idea what this does
                workd.as_ptr_mut() as *mut f64,
                workl.as_ptr_mut() as *mut f64,
                lworkl as i32,
                &mut info,
            );

            match info
            {
                0 | 1 | 2 =>
                {
                    println!("a {}", ido);
                }
                -1 => panic!("N must be positive."),
                -2 => panic!("NEV must be positive."),
                -3 => panic!("NCV-NEV >= 2 and less than or equal to N."),
                -4 => panic!("Maximum iterations must be greater than 0."),
                -5 => panic!("Maximum iterations must be greater than 0."),
                i => panic!("dsaupd_c returned error code {}", i),
            }

            if (ido == -1) || (ido == 1)
            {
                let res = transform(workd.subrows(ipntr[0] as usize - 1, dimension));
                workd
                    .subrows_mut(ipntr[1] as usize - 1, dimension)
                    .copy_from(&res);
                break;
            }
        }
    }

    let select = vec![false as i32; number_of_basis_vectors];
    let mut d: Col<f64> = Col::zeros(number_of_eigenvalues + 1);
    let mut z: Mat<f64> = Mat::zeros(dimension, number_of_basis_vectors);
    unsafe {
        dseupd_c(
            vectors as i32,
            "A".as_ptr() as *const i8,
            select.as_ptr(),
            d.as_ptr_mut() as *mut f64,
            z.as_ptr_mut() as *mut f64,
            dimension as i32,
            0.0 as f64,
            "I".as_ptr() as *const i8,
            dimension as i32,
            "LR".as_ptr() as *const i8,
            number_of_eigenvalues as i32,
            f64::EPSILON,
            residual.as_ptr_mut() as *mut f64,
            number_of_basis_vectors as i32,
            eigenvectors.as_ptr_mut() as *mut f64,
            dimension as i32,
            iparam.as_mut_ptr(),
            ipntr.as_mut_ptr(),
            workd.as_ptr_mut() as *mut f64,
            workl.as_ptr_mut() as *mut f64,
            lworkl as i32,
            &mut info,
        );
    }

    (d.subrows(0, number_of_eigenvalues).to_owned(), z)
}

And I am calling this this way:

faer_arpack_symmetric_f64(
        |v| mat * v,
        mat.nrows(),
        Eigenvaluekind::LargestMagnitude.as_str(),
        2,
        4,
        100,
        true,
    )

With a matrix that I know has an interesting eigendecomposition. However, the result I get at the end is 0 for the eigenvalues and eigenvectors that I requested. The input matrix is not 0, I have checked. And the operation

            let res = transform(workd.subrows(ipntr[0] as usize - 1, dimension));
            workd
                .subrows_mut(ipntr[1] as usize - 1, dimension)
                .copy_from(&res);

Is returning non-zero vectors on each iteration. I am also certain that at the end of all iterations, workd contains many non-zero elements. Although some of its subarrays are contiguous rows of 0.

I suspect I am probably making a mistake in the ways I am copying the data or passing the pointers, but I am not sure where.

Answer 1

The problem was in my use of dseupd, I was giving it the parameter:

"LR".as_ptr() as *const i8,

What I should have used is:

"LM".as_ptr() as *const i8,

Additionally, one should be careful that the memory layout of the high level matrix structures is continous.