Many bugs and incorrect documentation in {s,d}lasd{1,2,3}, {s,d}lasd{6,7,8}, and {s,d}laed2

[poulson]

The documentation for the deflation process in {s,d}lasd2 claims that the deflated singular values are sorted in increasing order, but they are in fact more-or-less sorted in decreasing order, with this constraint broken in certain cases when there was a small update vector entry deflation before a deflation of index JPREV due to a small diagonal difference deflation. As can be seen from {s,d}lasd1's call to {s,d}lamrg, the deflated portion of D is claimed to be in decreasing order, which turns out not to be necessarily true for the mentioned reason. I would therefore recommend calling DLASRT rather than DLAMRG to avoid breaking the assumption that the two subsets are already sorted.

Also, there is no d(j)-d(0) <= tol check within the initial deflation loop starting at line 425 of dlasd2. One could easily fix this by checking if d(j) is less than the tolerance rather than after checking if |z(j)| is, but this requires a Givens rotation to mix the first and j'th columns of U (and V), which destroys the fact that the first column of U2 would only have a single nonzero entry (in general, said column would become dense), so the packing logic would need to be changed to incorporate a rank-one update rather than copying a single row of the secular left singular vectors into the updated left singular vectors.

Further, the writes to DSIGMA during the deflation loop in {s,d}lasd2 are superfluous given the subsequent filling of DSIGMA at line 555 of dlasd2.

Further, no absolute value should be required on line 457 of dlasd2 (and the analogue in slasd2) since the diagonal entries are already sorted.

Further, the check on line 568 of dlasd2 of whether abs(DSIGMA(2)) <= tol/2 is strange for two reasons:

1. DSIGMA(2) >= 0, so the absolute value is superfluous
2. If DSIGMA(2) < tol, then it must have been deflated due to the small diagonal difference condition on line 457 . But then DSIGMA(2) is essentially the largest singular value due to the fact that the deflated components in DSIGMA are more or less stored in descending order. But this appears to be a contradiction given the definition of the deflation tolerance unless the matrix was of a sufficiently small dimension for the perturbations of a descending order to significantly break this argument (e.g., a pathological 3x3 matrix). Either way, some sort of comment is warranted for this check given the many subtle errors in the current implementation. Fixing the above-mentioned missing deflation criterion for the d(j) <= tol might obsolete this check.

Also, the documentation for U and V in {s,d}lasd3 incorrectly suggests that only the last N-K vectors are relevant, when, on exit, they provide the singular vectors for the merged bidiagonal problem.

[poulson]

I will submit a patch after I figure out how to appropriately handle the copyright notice given constraints imposed by my current employer.

[poulson]

It's worth noting that these deflation comments apply equally well to {s,d}lasd{6,7,8}.

[poulson]

{s,d}laed2 do not incorporate the norm of the rank-one update when choosing the deflation tolerance. In particular, they only look at the max-norm of z, where || z ||_2 = sqrt(2), and at the maximum absolute eigenvalue of the two subproblems. Thus, if a split happens to occur where the off-diagonal entry dominates the two subproblems, then the tolerance will be unnecessarily tight.