I am trying to get invert shift mode working on ARPACK, and I am interfacing the library through rust.
I am using this laplacian matrix for testing:
2, -1, 0, 0, 0, -1, 0, 0, 0, 0
-1, 3, -1, 0, 0, 0, -1, 0, 0, 0
0, -1, 3, -1, 0, 0, 0, -1, 0, 0
0, 0, -1, 3, -1, 0, 0, 0, -1, 0
0, 0, 0, -1, 3, -1, 0, 0, 0, -1
-1, 0, 0, 0, -1, 3, -1, 0, 0, 0
0, -1, 0, 0, 0, -1, 3, -1, 0, 0
0, 0, -1, 0, 0, 0, -1, 3, -1, 0
0, 0, 0, -1, 0, 0, 0, -1, 3, -1
0, 0, 0, 0, -1, 0, 0, 0, -1, 2
For the regular mode I am getting the correct eigenvalues according to scipy. However if I try switching to inverse shift mode I get incorrect values. I tried setting it up in rust as follows:
fn faer_arpack_symmetric_f64(
mat: &SparseColMat<usize, f64>,
dimension: usize,
eigenvalue_kind: &str,
number_of_eigenvalues: usize,
number_of_lanczos_vectors: usize,
maxiter: usize,
vectors: bool,
) -> (Col<f64>, Mat<f64>)
{
let mut ido = 0;
let mut residual: Col<f64> = Col::zeros(dimension);
let mut arnoldi_data = vec![0.0; dimension * number_of_lanczos_vectors];
let mut iparam = [0; 11];
iparam[0] = 1;
iparam[2] = maxiter as i32;
iparam[6] = 3; // code 3 means inverse shift according to the book
let mut ipntr = [0; 11];
let mut workd = Col::zeros(3 * dimension);
let lworkl = number_of_lanczos_vectors * (number_of_lanczos_vectors + 8);
let mut workl = Col::<f64>::zeros(lworkl);
let mut info = 0;
let sigma = 0.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_lanczos_vectors as i32,
arnoldi_data.as_mut_ptr() 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,
);
match info
{
0 =>
{}
1 =>
{
panic!("Maximum number of Lanczos iterations reached.")
}
3 =>
{
panic!(
"No shifts could be applied during implicit Arnoldi update, try increasing NCV."
)
}
-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 = mat * workd.subrows(ipntr[0] as usize - 1, dimension);
let y = faer_shift_invert_solve(mat, sigma, &res);
workd
.subrows_mut(ipntr[1] as usize - 1, dimension)
.copy_from(&y);
}
else if ido == 99
{
break;
}
else
{
panic!("ido code: {}", ido);
}
}
}
let select = vec![false as i32; number_of_lanczos_vectors];
let mut d: Col<f64> = Col::zeros(number_of_eigenvalues + 1);
let mut z_storage = vec![0.0; dimension * number_of_lanczos_vectors];
unsafe {
dseupd_c(
vectors as i32,
"A".as_ptr() as *const i8,
select.as_ptr(),
d.as_ptr_mut() as *mut f64,
z_storage.as_mut_ptr() as *mut f64,
dimension as i32,
sigma,
"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_lanczos_vectors as i32,
arnoldi_data.as_mut_ptr() 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,
);
}
match info {
0 => { /* Normal exit, do nothing */ }
1 => panic!("INFO = 1: The Schur form computed by LAPACK routine dsytrd failed to reorder or to compute eigenvalues and eigenvectors to high relative accuracy."),
3 => panic!("INFO = 3: No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge."),
-1 => panic!("INFO = -1: N must be greater than 0."),
-2 => panic!("INFO = -2: NEV must be greater than 0 and less than N."),
-3 => panic!("INFO = -3: NCV must be greater than NEV and less than or equal to N."),
-5 => panic!("INFO = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'."),
-6 => panic!("INFO = -6: BMAT must be one of 'I' or 'G'."),
-7 => panic!("INFO = -7: The part of the spectrum specified by WHICH is not supported."),
-8 => panic!("INFO = -8: The spectrum range parameters are inconsistent."),
-9 => panic!("INFO = -9: WORKD[1:N] is not correctly allocated."),
-10 => panic!("INFO = -10: The starting vector is zero."),
-11 => panic!("INFO = -11: The standard Lanczos iteration (for generalized symmetric eigenvalue problem) failed to converge."),
-12 => panic!("INFO = -12: The shifts could not be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge."),
_ => panic!("Unknown INFO value: {}", info),
}
let mut eigs = d.subrows(0, number_of_eigenvalues).to_owned();
(
eigs,
faer::mat::from_column_major_slice_mut(
&mut z_storage,
dimension,
number_of_lanczos_vectors,
)
.to_owned(),
)
}
fn faer_shift_invert_solve(
a: &SparseColMat<usize, f64>,
sigma: f64,
b: &Col<f64>,
) -> Col<f64>
{
let triplets = (0..a.nrows()).map(|i| (i, i, 1.)).collect::<Vec<_>>();
let identity =
SparseColMat::try_new_from_triplets(a.nrows(), a.nrows(), &triplets).unwrap();
use faer::prelude::SpSolver;
let qr = (a - sigma * identity).sp_qr().unwrap();
qr.solve(b)
}
```
