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. Pauli