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.
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
UPLOflag 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:
Using
UPLO::Upperdoes not work, even though the matrix is clearly upper-diagonal inndarray-world.It looks like this stems from this method, perhaps we should transpose
uploin 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 calluplo.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.