Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lapack integration still buggy for complex eigenvalues #1105

Open
geoeo opened this issue May 1, 2022 · 6 comments
Open

Lapack integration still buggy for complex eigenvalues #1105

geoeo opened this issue May 1, 2022 · 6 comments

Comments

@geoeo
Copy link
Contributor

geoeo commented May 1, 2022

In the notes of nalgebra v. 31 and lapack v.22 it states improved behavior when dealing with the computation of eigenvalue problem.
While I can confirm this for nalgebra, the lapack variant still crashes.

This works: action_matrix.complex_eigenvalues()
This crashes: nalgebra_lapack::Eigen::complex_eigenvalues(action_matrix)

The lapack integration works since I use it use compute SVD etc.

The error is

 ** On entry to DGEEV  parameter number  1 had an illegal value
The application panicked (crashed).
Message:  Lapack error.
Location: /home/marc/.cargo/registry/src/github.com-1ecc6299db9ec823/nalgebra-lapack-0.22.0/src/eigen.rs:272

Tested on Ubuntu 20.04 (WSL)

@geoeo
Copy link
Contributor Author

geoeo commented May 2, 2022

The parameters passed in the unsafe block are invalid.
Looking at the source once call is e.g.

        let lwork = T::xgeev_work_size(
            b'T',
            b'T',
            n as i32,
            m.as_mut_slice(),
            lda,
            wr.as_mut_slice(),
            wi.as_mut_slice(),
            &mut placeholder1,
            n as i32,
            &mut placeholder2,
            n as i32,
            &mut info,
        );

        lapack_panic!(info);

        let mut work = vec![T::zero(); lwork as usize];

        T::xgeev(
            b'T',
            b'T',
            n as i32,
            m.as_mut_slice(),
            lda,
            wr.as_mut_slice(),
            wi.as_mut_slice(),
            &mut placeholder1,
            1 as i32,
            &mut placeholder2,
            1 as i32,
            &mut work,
            lwork,
            &mut info,
        );
        lapack_panic!(info); 

According to the docs http://www.netlib.org/lapack/explore-html/d9/d8e/group__double_g_eeigen_ga66e19253344358f5dee1e60502b9e96f.html the only valid chars are 'N' or 'V'.
When I change it to e.g. 'V' I get

 ** On entry to DGEEV  parameter number  9 had an illegal value
The application panicked (crashed).
Message:  Lapack error.
Location: /home/marc/Rust/nalgebra/nalgebra-lapack/src/eigen.rs:292

A different error in the second call! Progress!

@geoeo
Copy link
Contributor Author

geoeo commented May 2, 2022

It seem all the work is already done expect the complex part is discarded in the normal eigenvalue() call.
I changed the complex eigenvalue method signature to

    pub fn complex_eigenvalues(mut m: OMatrix<T, D, D>, left_eigenvectors: bool, eigenvectors: bool) 
        -> (OVector<Complex<T>, D>,  Option<OMatrix<Complex<T>, D, D>>,  Option<OMatrix<Complex<T>, D, D>>)

But apparently the check macro doesnt like that

macro_rules! lapack_check(
    ($info: expr) => (
        // TODO: return a richer error.
        if $info != 0 {
            return None;
        }
        // if $info < 0 {
        //     return Err(Error::from(ErrorKind::LapackIllegalArgument(-$info)));
        // } else if $info > 0 {
        //     return Err(Error::from(ErrorKind::LapackFailure($info)));
        // }
    );
);

Any suggestions @Andlon ?

@geoeo
Copy link
Contributor Author

geoeo commented May 7, 2022

The accuracy of this method does not seem to be very high when compared to matlab/octave This is not really true.

@geoeo geoeo closed this as completed May 7, 2022
@Andlon
Copy link
Sponsor Collaborator

Andlon commented May 9, 2022

@geoeo: did you mean to close this? It looks like you've discovered an important deficiency in nalgebra-lapack here.

I'm otherwise unfortunately not familiar with the codebase of nalgebra-lapack and I'm sorry to say that I don't have the bandwidth to dig into it at the moment. Hopefully someone else is able to step up, however.

@geoeo
Copy link
Contributor Author

geoeo commented May 12, 2022

@Andlon I made a PR. maybe its better to discuss this there. #1106

@Andlon
Copy link
Sponsor Collaborator

Andlon commented Sep 12, 2022

I've reopened this because the linked PR (#1106) does not directly address this issue.

@Andlon Andlon reopened this Sep 12, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants