Allow gesdd to only compute one of U and V

[chtobler]

Hello,

In the gesvd function, it was possible to select jobu and jobvt separately, allowing us to specify that we want to compute U but not VT, or VT but no U. With gesdd, there is just one option jobz that will compute either none of U and VT, or both U and VT.

This causes some difficult decisions of which function to use, since gesdd is generally faster, but the additional wasted work of computing an output just to discard it can make it slower to use after all.

In MATLAB, for a while we were using gesvd for the 2-output syntax [U, S] = svd(A) and gesdd for all other cases. However, this was causing confusion with users, since U and S would change between the calls [U, S] = svd(A) and [U, S, V] = svd(A). We have changed svd to always use gesdd to avoid this confusion, but this comes at a cost of performance, most noticeably when calling orth(A) to find an orthonormal basis of the columns of A.

It would be ideal for us if LAPACK provided an interface to request that gesdd only computes U, or only computes V.

Some ideas for a new syntax:

• New function (gesdds with the "s" standing for separate U/V computation?) which takes inputs jobu and jobvt instead of only jobz. This could be called internally from gesdd, which would set jobu=jobz, jobvt=jobz, to keep only one implementation.

• Make the jobz input to gesdd additionally accept a length-3 character array, the first character could be 'M' for "mixed", followed by one of 'A', 'S', 'O', 'N' for U and another character of the same set for V. I added a new first character, because a C-style API would need some first character that currently errors to tell it that it's save to access additional characters after the jobz pointer.

Note: I imagine some of the reason for gesdd to only have the one jobz option was that some combinations of the 'A', 'S' and 'O' options for U and VT could be difficult to compute and were unlikely to be used. I would be perfectly happy if the mixed case requires one of jobu and jobvt to be 'N'. I think it would also be fine to have additional restrictions that 'O', in the mixed case, can only be used for either 'U' or 'V' depending on the matrix size.

Such a new interface would be extremely helpful for us, allowing for better performance in two-output calls to svd, and for simpler code where we call gesdd.

Thank you,
Christine

[mgates3]

I think bdsdc (bidiagonal divide and conquer) necessarily always computes both U and V. My understanding is it is required by the D&C algorithm. An optimization is not computing U or V vectors at the root node of the D&C tree, only in the lower levels. This could potentially save considerable time in bdsdc, since the root node is the most expensive node.

Outside of bdsdc, if only U is desired, gesdd could skip the back-transformation of V to save some time (and similarly for V, skipping back-transform of U).