Thiết kế website giá rẻ

Question

I am trying to use arpack from within rust and I am getting weird results.

In order to test if I was correctly interfacing with arpack, I decided to make a sparse symmetric tri-diagonal Toeplitz matrix by setting all diagonal values to 6 and all off diagonal values to 2.

My understanding is that for an nxn such matrix the eigenvalues will be 6 – 4 cos(i * PI / (n + 1)) for i in [1, n]. So I wrote a test to verify if that was my output, and I am getting (for 5 eigenvalues), this result:

// Expected:
5.594633205988414
5.189266411976828
4.7838996179652415
4.378532823953655
3.973166029942069
3.5677992359304826
3.1624324419188965
2.7570656479073103
2.351698853895724
1.9463320598841376
1.5409652658725523
1.1355984718609653
0.7302316778493791
0.324864883837793
-0.08050191017379227
-0.4858687041853793
-0.8912354981969663
-1.2966022922085516
-1.7019690862201378
-2.107335880231725
-2.512702674243311
-2.9180694682548953
-3.3234362622664815
-3.7288030562780694
-4.134169850289656
-4.539536644301242
-4.944903438312828
-5.350270232324414
-5.755637026336
-6.1610038203475845
// Obtained
[
[9.497386464578327],
[9.675831246480941],
[9.816557025600208],
[9.918119765009969],
[9.97947729356757],
]

So the results are not even close to what they should be.

This is the smallest MVE I could make:

use faer::sparse::SparseColMat;
use faer::*;

extern "C" {
    pub fn dsaupd_c(
        ido: *mut ::std::os::raw::c_int,
        bmat: *const ::std::os::raw::c_char,
        n: ::std::os::raw::c_int,
        which: *const ::std::os::raw::c_char,
        nev: ::std::os::raw::c_int,
        tol: f64,
        resid: *mut f64,
        ncv: ::std::os::raw::c_int,
        v: *mut f64,
        ldv: ::std::os::raw::c_int,
        iparam: *mut ::std::os::raw::c_int,
        ipntr: *mut ::std::os::raw::c_int,
        workd: *mut f64,
        workl: *mut f64,
        lworkl: ::std::os::raw::c_int,
        info: *mut ::std::os::raw::c_int,
    );
}

extern "C" {
    pub fn dseupd_c(
        rvec: ::std::os::raw::c_int,
        howmny: *const ::std::os::raw::c_char,
        select: *const ::std::os::raw::c_int,
        d: *mut f64,
        z: *mut f64,
        ldz: ::std::os::raw::c_int,
        sigma: f64,
        bmat: *const ::std::os::raw::c_char,
        n: ::std::os::raw::c_int,
        which: *const ::std::os::raw::c_char,
        nev: ::std::os::raw::c_int,
        tol: f64,
        resid: *mut f64,
        ncv: ::std::os::raw::c_int,
        v: *mut f64,
        ldv: ::std::os::raw::c_int,
        iparam: *mut ::std::os::raw::c_int,
        ipntr: *mut ::std::os::raw::c_int,
        workd: *mut f64,
        workl: *mut f64,
        lworkl: ::std::os::raw::c_int,
        info: *mut ::std::os::raw::c_int,
    );
}

pub fn eigen_decompose_sparse_matrix(
    mat: &SparseColMat<usize, f64>,
) -> (Col<f64>, Mat<f64>)
{
    debug_assert!(mat.nrows() == mat.ncols());
    faer_arpack_symmetric_f64(
        |v| mat * v,
        mat.nrows(),
        Eigenvaluekind::LargestMagnitude.as_str(),
        5,
        10,
        3000,
        true,
    )
}

pub enum Eigenvaluekind
{
    LargestMagnitude,
    SmallestMagnitude,
    LargestRealPart,
    SmallestRealPart,
    LargestImaginaryPart,
    SmallestImaginaryPart,
}

impl Eigenvaluekind
{
    fn as_str(&self) -> &'static str
    {
        match self
        {
            Eigenvaluekind::LargestMagnitude => "LM",
            Eigenvaluekind::SmallestMagnitude => "SM",
            Eigenvaluekind::LargestRealPart => "LR",
            Eigenvaluekind::SmallestRealPart => "SR",
            Eigenvaluekind::LargestImaginaryPart => "LI",
            Eigenvaluekind::SmallestImaginaryPart => "SI",
        }
    }
}

// http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf
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 arnoldi_vectors = 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 = number_of_basis_vectors * (number_of_basis_vectors + 8);
    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 =>
                {}
                -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);
            }
            else if ido == 99
            {
                break;
            }
            else
            {
                panic!("ido code: {}", ido);
            }
        }
    }

    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,
            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(),
            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),
    }

    use crate::sparse_vector::stringify_mat;
    // println!("{:?}", residual);
    // println!("{}", stringify_mat(&z, 2));
    (d.subrows(0, number_of_eigenvalues).to_owned(), z)
}

fn faer_shift_and_invert(mat: &SparseColMat<usize, f64>, col: &Col<f64>, sigma: f64)
{
    // let shifted_mat = mat - sigma * ;
}

// +| Tests |+ =======================================================
#[cfg(test)]
mod tests
{
    use super::*;

    #[test]
    fn test_sparse_eigen_decomposition()
    {
        // Symmetric tri diagonal toeplitz matrix.
        // https://en.wikipedia.org/wiki/Tridiagonal_matrix#Eigenvalues
        let mut triplets = Vec::new();
        let n = 30;
        let a = 6.;
        let b = 2.;
        for i in 0..n
        {
            triplets.push((i, i, a));
            if i < n - 1
            {
                triplets.push((i, i + 1, b));
            }
            if i > 0
            {
                triplets.push((i, i - 1, b));
            }
        }
        let mat = SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap();

        let (eigs, eigvs) = eigen_decompose_sparse_matrix(&mat);

        for i in 0..n
        {
            use std::f64::consts::PI;
            let toeplitz_val = a - 2. * b * ((i + 1) as f64 * PI / (n as f64 + 1.));
            println!("{}", toeplitz_val);
        }
        use crate::sparse_vector::stringify_sparse_mat;
        println!("{:?}", eigs);
    }
}

Thiết kế website giá rẻ

Danh mục

ARPACK producing nonsense eigenvalues?