The point is: GESVD is also returning wrong result for the original example. i.e. you only gain extra precision of one or two digits after the switch to GESVD. see (4) below, and my uploaded file #46268.
@Tim Mitchell
1. Sorry for my misleading words. Yes, I mean I'm suggesting to use GESDD as the default driver, and keep GESVD as an alternative, and a proper documentation in svd().
The reason is, as I said and as you understood: (i) error=1e15 is not too bad; (ii) GESVD too slow for large scale problem.
The speed is alread discussed elsewhere, but let me post one more sample.
a = randn(2000); tic; [u,s,v]=svd(a); toc;
octave, GESVD: 54.822 sec
octave, GESDD: 3.6 sec
matlab: 3.2 sec
We gain extra (say 2 digits) precision occasionally at the cost of 10 times slower, every time, despite of the user need it or not (assume average users keep using the default driver). Most users will feel that Octave is slow, instead of been more accurate.
2. For the dgesvj or dgejsv drivers. I will file a new bug.
3. Thanks for linking the SIAM Review. It helps a lot. Thanks for correcting the relative error computations.
4. Right, GESDD returns wrong result there, but GESVD is also wrong for the same example. The second singular value of A:
N = 26
formula solution = 1
GESVD = 3.76175e+09
GESDD = 6.07965e+10
Clearly, this example can't be used to demonstrate that GESVD is better than GESDD. They are all wrong, well GESVD is less wrong by one digit.
Compared to the first singular value s_max = 6.08901e+26, they are all around or under machine epsilon.
In the SIAM review, fig.29 and fig.30 compared the accuracy for us. The improvement of QR (GESVD) compare to D&C (GESDD) range from zero to several digits. I think it is not worth to switch for this small and uncertain improvement.

@count: Go ahead and file a different bug report about adding additional drivers to svd_drivers() as a feature request. If you have software skills it would also be good to submit a patch. This is a relatively esoteric request which may not attract anyone else to write such a patch.
Beyond that, no algorithm is perfect. Octave is prioritizing accuracy over speed, butunlike Matlabwe give you the choice of using whatever driver you want. I think the documentation, quoted below, is clear on that.

@count:
1. I don't think anyone is proposing to remove any of the different SVD drivers. We are just discussing to not use GESDD by default (because of its worse accuracy than GESVD).
2. Yes, in my #15, I already mentioned that both dgesvj and dgejsv are available in LAPACK. Perhaps it would be nice if Octave could also add support for them to be chosen as other alternatives if the user wishes.
3. Regarding demo2(), sorry about changing the order of the log space and yes, now I can replicate that. However, again no one was saying GESVD is perfect, just that it is much better than GESDD; again, see the SIAM Review article. Also, your update completely broke the relative error computations. Below is a corrected version. For N=26, while both methods indeed struggle with the small singular values, GESDD is still much worse (1e23 relative error versus 1e11). Such an example would be a reason to try the Jacobi routines.
4. GESDD does not fail at around 1e15; it just seems to return any singular value below norm(A)*eps as about norm(A)*eps, regardless of its true value. This means such singular values can be completely wrong (as in no digits accurate) and this can even happen to large singular values if norm(A) is very large. That is demonstrated in the original example.
If risking such an accuracy loss for some speedup is fine for your needs, great. But that tradeoff doesn't seem like a good default choice for everyone, particularly since there are lots of applications where the higher accuracy is needed and expected.

...In comment #17, the typo fix is for my #15
@Tim Mitchell
For demo1(), though 'gdesdd' start to fail around 1e15. I don't feel too bad about it. It is around machine epsilon.
For demo2(), don't change the order of the logspace(), my code is for demo that 'dgesvd' is not better than 'dgesdd' too much.
It should looks like this (note the logspace() range and the transpose of s):
My Matlab version is 9.3.0.713579 (R2017b) under Linux, using CPU i76700HQ (with avx2 and fma). No idea why can't be reproduced.

I forgot to mention that I vote for keep using the 'gesdd' driver.
There is one more SVD solver in LAPACK: 'dgejsv', which uses Jacobi SVD algorithm. It should be the most accurate one in LAPACK. If possible, I wish to add that as an addition svd_driver. Should I open another ticket/issue for this?
Fix for some typos in comment #16 ...
"demo 'gesdd' is better than 'gesvd'"
it is a demo of 'gesvd' is better than 'gesdd':
In the attachement (file #46268), The comment
% Characteristic polynomial of A*A', note svd(A)=eig(A*A')=roots(c).
should be
% Characteristic polynomial of IA*A', note svd(A)^2=eig(A*A')=1roots(c).

@count: for me, your two demos show GESVD being much more accurate than GESDD. Perhaps someone else here can test too.
In your first demo, you accidentally reversed s1 and s2. On R2017b, for N=26, it actually shows that GESVD returns the B matrix (exactly, not surprisingly) while GESDD incurs large errors in the singular values (as in 1e14 relative errors bad). Here is a corrected version, which outputs the singular values and then the relative errors:
On R2017b, I can't replicate the bad results you report for N=11. For N=26, GESVD computes each singular value to at least 11 digits while GESDD has relative errors up to 10e23. Here is the updated code, with same output updates I made for demo1:
In all the cases and tests I've seen so far, GESVD is much better than GESDD in terms of relative accuracy of the singular values, even though it too is not perfect. See Figs. 29 and 30 of the recent SIAM Review article on the SVD by Dongarra et al, which precisely demonstrate how both can be very bad but GESDD is generally much worse accuracywise on poorly scaled matrices. Meanwhile, the figures also show how the Jacobi and QRP alternatives do handle such poorly scaled problems to high relative accuracy, but from Figs. 2728, are generally slower than GESVD. By the way, LAPACK already has at least two Jacobibased SVD routines: dgesvj and dgejsv (both are onesided Jacobi, the latter with preconditioning).

Just wait a minute...
I found the following:
 1) 'gesvd' did not give the "correct" result either, for N=26.
I worked out an explicity formula for svd(A) in this case. So we can measure the relative error now. See attachment for code. A point is, the second singular value should be 1.0, but we get 3.76175e+09 from 'gesvd'.
A straightforward example that demo 'gesdd' is better than 'gesvd':
Here svd(B) is trivial to compute, yet, 'gesdd' is not that satisfactory. And indeed, we see that 'gesvd' is carefully programmed.
However, here is a frustrating example that convince you that 'gesvd' is no much better:
We know exactly the SVD of C by definition: U=I, S=diag(scaling), V=A0'. Yet, 'gesvd' is disappointed, in an idealism scene.
 2) Matlab(2017b) also returns wrong results, the numbers (both the correct and the wrong) are the same as octave. Indicating that Matlab is also using 'gesdd'.
 3) The error of 'gesdd' is acceptable.
Note that the ratio between the second and the first svd(A) is about 1.0e16, which is coincede with the machine epsilon. Similary for svd(B) and svd(C) above. This is the precision we can expect, but no more. As 'gesvd' could also produce wrong results around 1.0e17 (in the case of svd(C) above and svd(A) in attachement), there is no urgent need to pursuit the occasional extra precision.
 4) Speed is important, given that 'gesdd' is not terribly wrong.
So, I suggest document the difference between 'gesvd' and 'gesdd', let the careful user pick the "correct" one.
PS: I heard that the "Jacobi" method is an even more accurate method (ref), but I don't have a good code in hand that can handle any of above extreme case.
(file #46268)

I reopened bug #49940 and linked it to this report.

I made "gesvd" the default driver in Octave (https://hg.savannah.gnu.org/hgweb/octave/rev/1212568010a8). I updated the documentation to explain the issue and with a URL reference to this bug report.

@Mike: Given that the inaccuracies of GESDD's results can be arbitrarily bad, I'd personally be wary to add any documentation encouraging users to try it, even if speed is valued more than accuracy. If anything, until this is better resolved on the LAPACK side, I think it would be prudent to warn users that using GESDD could be trading off all accuracy for speed!

@Mike: At the moment, no other solution than documentation. We could look to see if there is anything within Octave, such as unnecessarily duplicating matrices, that we might tweak. But I believe the difference in performance is really located with LAPACK.

From my understanding, bug #49940 is about the performance of the GESVD function, which we are now considering making the default again. Is there any solution for that other than documenting it and suggesting that users who care more about speed than accuracy should use GESDD?

I generally agree. I guess my point was that we switched to gesdd
to resolve bug 49940. So if we go back (and I agree that at this
moment we should), we need to reopen 49940.
Also the benchmark ratio seems got worse.
Dmitri.

We can look at bug #49940 again. But, I still believe we need to prioritize correctness over speed. The issue with this bug is that most people will be unaware that they have incorrect results. Even the pessimists, who verify the results using the acceptance test in the documentation, will believe that the results are good. It is only if you calculate using GESVD that you will find the problem. Since any new acceptance test is really going to require running GESVD, we might as well do that in the first place.
I reread the discussion in bug #49940, and I'm pretty sure it isn't Octave's problem. My guess is that the GESDD "Divide and Conquer" algorithm is able to spread the problem over multiple threads/cores, but arrives at the wrong answer. And the GESVD algorithm is probably old and singlethreaded with few opportunities for parallelism.

May be we should revisit bug 49940.
It looks to me that gesvd mostly runs singlethread
even with multithreaded blas.
I run a simple benchmark and watched top. With
gesdd cpuload was 400% (4 threads), with
gesvd it started ~400% but dropped quickly to 100%
similar results with liopenblasp (pthread interface) and atlas.
I am not sure if this is octave problem (posted results in 49940
did not show such a big difference), or lapack's.
Dmitri.

Sure, seems prudent to me. The default was last changed in 2016 for the 4.4.0 release.
I get identical test suite results with both drivers, and confirm switching to GESVD fixes the problem reported here.

@Mike: I think we should favor correctness over speed. What do you think about changing the default driver in Octave back to GESVD for the 5.1 release?

@Mike: Thanks. I renamed the function to tst_svd on my platform. Here is a sample session which shows changing the backend SVD function to GESVD fixes the issue.

Copy/paste error, this should be a little clearer:

Here is what I've distilled out of the original problem statement on https://github.com/ReferenceLAPACK/lapack/issues/316:

Is the original problem matrix available in a file (.mat or otherwise)? Could that be uploaded?

I recently opened an issue on LAPACK's Github page about how the GESDD routine can produce vastly different singular values if the singular vectors are also requested. At the moment, it is still unclear if this behavior is a coding bug or just an inherent deficiency of the "divide and conquer" approach that GESDD implements. You will need to follow the LAPACK issue for developments on that and for complete details.
As Octave 4.4 is using GESDD by default (like MATLAB and others), I just wanted to raise your attention of how bad this problem can be in practice. For future Octave releases, you might consider only using GESVD for all svd calls, unless the user specifically enables GESDD (at their own risk). ;)
In any case, GESDD's inaccuracy in singular values can basically be arbitrarily inaccurate. To be fair, note that this may only happen rarely and indeed the example I found was designed to be numerically challenging (companion_demo(n) from EigTool, with n>=26). Nevertheless, even though this matrix is nearly singular (GESVD returns 2.6e10 as the smallest singular value), GESDD instead returns 2.3e+9 as the smallest singular value! For this problem, only the largest singular value is computed consistently; all others don't even agree to a single digit between GESVD and GESDD.
Lastly, note that using norm(U*S*V'A)/norm(A) to verify the computation for this example is not a sufficient test. The reason is that the largest singular value is 6e+25 and all the others are at least 1e+16 times smaller. So, to machine precision, the matrix seems to be rank one. In other words, one can randomly choose any value in [0,1e9] for the other n1 singular values and still have norm(U*S*V'A)/norm(A) look good numerically! A better test is to compare the singular values to those computed by GESVD and make sure each singular value has a high number of digits in agreement.
