On 4/26/11 11:49 AM, Pauli Virtanen wrote:
Tue, 26 Apr 2011 11:36:19 -0500, Jason Grout wrote: [clip]
Okay, just one more data point. Our people that are seeing the problem with numpy returning a non-unitary V also see a non-unitary V being returned by the test C call to zgesdd. In other words, it really appears that zgesdd follows the Lapack docs, and if rows<columns, the returned V is not necessarily unitary, but may contain a zero row. This makes numpy's assumptions in using zgesdd false.
You can see this report at http://trac.sagemath.org/sage_trac/ticket/11248#comment:25
What LAPACK promises is not very clear. Earlier on the the man page:
JOBZ (input) CHARACTER*1 Specifies options for computing all or part of the matrix U: = 'A': all M columns of U and all N rows of V**H are returned in the arrays U and VT; = 'S': the first min(M,N) columns of U and the first min(M,N) rows of V**H are returned in the arrays U and VT; = 'O': If M>= N, the first N columns of U are over‐ written in the array A and all rows of V**H are returned in the array VT; otherwise, all columns of U are returned in the array U and the first M rows of V**H are overwritten in the array A; = 'N': no columns of U or rows of V**H are computed.
Looks to me the statement should be parsed as:
return_n_rows = (jobz == 'A') or (jobz == 'O' and m>= n)
So the current usage should be OK. (At least, as long as jobz == 'A', in the other cases, I don't think the wrapper does the right thing.)
But apparently either there's a bug, or the LAPACK man page needs to be understood as you say.
Ah, you're right that it makes sense to parse their statement that way too, so I'm not so sure what Lapack really is saying anymore either. If it's parsed the way you propose (which makes sense given the JOBZ description), I think it points to a bug in the Accelerate Lapack on the affected platforms, as we get the same zero row when we call the function directly from C, without numpy, python, or Sage in the middle. The updated rwork calculation makes no difference with a 3x4 matrix (both the old calculation and the new calculation give 66 in the 3x4 case), so I don't think that is affecting anything. Jason