nmatrix
nmatrix copied to clipboard
SciRuby
LAPACK getrf/getri function and spec broken (and causes #det to be wrong)
getri is supposed to compute the inverse of a matrix. In the spec we test getri on
a = NMatrix.new(:dense, 3, [1,0,4,1,1,6,-3,0,-10], dtype)
and compare the result with
b = NMatrix.new(:dense, 3, [-5,0,-2,-4,1,-1,1.5,0,0.5], dtype)
We can check that a and b are not inverses of each other.
This is another case of feeding row-major arrays into LAPACK functions that expect column-major arrays.
I suspect that this is not the only LAPACK function that is broken. At some point, I will check them all and fix what needs to be fixed. For now, this probably isn't urgent because if anybody had been using this function they would have complained already.
Wait, I take that back, a and b actually are inverses. I definitely remember finding something wrong with these functions, I will have to check exactly what it was.
OK, it's easier to see if you just look getrf by itself. getrf is supposed to do an LU decomposition with output as described here (it's also the precursor to getri). If you look at our spec for getrf, you can check the example and see that it doesn't give a proper LU decomposition.
As in the spec is wrong, or the result doesn't match the spec but it's still passing?
OK, I was wrong about this. I was misinterpreting the output from getrf. According to the ATLAS convention for row-major matrices, the upper matrix is supposed to have unit diagonal elements, but I was assuming it was the lower matrix. I think the current behavior is correct. The det thing is a different issue, I'll submit a PR for that.
@wlevine please update the documentation accordingly. If you thought it was wrong, probably others will too.
Good point. The comments in the C code are accurate (which is how I figured out what was going on), but they don't seem to be reflected in the documentation. The documentation from the ruby code seems a little confused, and I can fix that. It looks like there is a way (http://rdoc.rubyforge.org/RDoc/Parser/C.html) to get rdoc to process C files, do we use this currently?
I think that we should use YARD and Rubydoc from now on. You can find NMatrix documentation on Rubydoc. (I probably should take the job of converting the docs to YARD...)
I'm not sure why the comments in C source aren't being reflected in the docs, by the way. I will need to investigate.