[Numpy-discussion] Zero row in SVD's unitary matrix on some Mac's

[Pauli Virtanen]

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