Skip to content

lax::UPLO references column-major layout #394

Open
@bdura

Description

@bdura

Hello everyone, thanks for the great work.

Disclaimer: I'm new to LAPACK & scientific computing in general, apologies if what follows is not perfectly rigorous.

I've noticed that the UPLO flag does not necessarily do what we expect, in the sense that it references the column-major alignment of the data. In particular, I'm not sure this aligns with newcomers' expectations.

To wit, in the following example:

extern crate ndarray;
extern crate ndarray_linalg;

use ndarray::*;
use ndarray_linalg::*;

fn main() {
    // A very clearly upper-diagonal matrix
    let a = arr2(&[[3.0, 1.0, 1.0], [0.0, 3.0, 1.0], [0.0, 0.0, 3.0]]);
    // Works with UPLO::Lower, not UPLO::Upper...
    let (e, vecs) = a.eigh(UPLO::Lower).unwrap();
    println!("eigenvalues = \n{:?}", e);
    println!("V = \n{:?}", vecs);
    let av = a.dot(&vecs);
    println!("AV = \n{:?}", av);
}

Using UPLO::Upper does not work, even though the matrix is clearly upper-diagonal in ndarray-world.

It looks like this stems from this method, perhaps we should transpose uplo in that case? Note that in the case of a tuple of arrays, we should probably assert that both arrays have the same layout, otherwise it does not make sense to call uplo.t().

It might make sense to decide that UPLO is column-major only, and leave the responsibility to the caller. But in that case I think we should document that behaviour in UPLO's documentation.

I'm happy to propose a PR, let me know which strategy seems best to you.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions