luhenry / netlib

An high-performance, hardware-accelerated implementation of Netlib in Java
Other
60 stars 12 forks source link

AbstractLAPACK dgeev argument checking #12

Closed uhoefel closed 2 years ago

uhoefel commented 2 years ago

There seems to be an issue in AbstractLAPACK.java in dgeev. My test input works with netlib.fommil and looks like this:

int n = 2;
double[] a = {1.0, 0.75, 0.5, 1.0};
double[] wr = new double[2];
double[] wi = new double[2];
int lwork = 68; // optimal size from previous query
double[] work = new double[lwork];
intW info = new intW(2);

LAPACK.getInstance().dgeev("N", "N", n, a, n, wr, wi, new double[1], 1, new double[1], 1, work, lwork, info);

To my understanding, that should give the following eigenvalues in wr:

$$\lambda_{1,2} = \frac{\pm\sqrt{6}+4}{4}$$

My suspicion is that lines 338-341

if (lsame("N", jobvl))
    checkIndex(offsetvl + n * ldvl - 1, vl.length);
if (lsame("N", jobvr))
    checkIndex(offsetvr + n * ldvr - 1, vr.length);

should read

if (lsame("V", jobvl))
    checkIndex(offsetvl + n * ldvl - 1, vl.length);
if (lsame("V", jobvr))
    checkIndex(offsetvr + n * ldvr - 1, vr.length);

The documentation of LAPACK states for vl that (similar for vr)

If JOBVL = 'N', VL is not referenced.

so I would at least not intuit that its length matters (or I may very well miss something).

PS: Thanks for fixing the split packages!

luhenry commented 2 years ago

That's right, let me submit your fix.

luhenry commented 2 years ago

Could you please verify https://github.com/luhenry/netlib/pull/15. If it seems right to you, I'll do a 3.1.0 release with that fix.

uhoefel commented 2 years ago

Seems correct as far as I can judge. I think the checks some lines above (e.g. 4962-4965, but also in dgeev) probably should be flipped to "V" as well:

    if (lsame("N", jobvl))
      requireNonNull(vl);
    if (lsame("N", jobvr))
      requireNonNull(vr);

-->

    if (lsame("V", jobvl))
      requireNonNull(vl);
    if (lsame("V", jobvr))
      requireNonNull(vr);
luhenry commented 2 years ago

It's been released in v3.0.2