function svdbug(n)
% Tim Mitchell
% mitchell@mpi-magdeburg.mpg.de
% Date: January 22, 2019
%
% Using the companion_demo from EigTool, with dimension n >= 26, svd
% returns very different values for the singular values when singular
% vectors are also requested.  The issue is due to the gesdd routine from
% LAPACK, which is called when singular vectors are requested.
% Interestingly, gesdd computes the correct result when the matrix
% dimension is 25 or smaller.
%
% I have observed this bug on both macOS (10.12.6) and Ubuntu 16.04.5 LTS.
% On Ubuntu, I have tested all of the MATLAB releases below.  On the Mac, I
% have only tested R2016b, R2017b, and R2018a.
%
% R2013b, R2014b, R2015b, R2016b, R2017a, R2017b: 
% - s = svd         WORKS
% - [U,S] = svd     WORKS
% - [U,S,V] = svd   FAILS
%
% R2018a: 
% - s = svd         WORKS
% - [U,S] = svd     FAILS
% - [U,S,V] = svd   FAILS
%
% The problem also exists in Octave 4.4.1 but apparently not in 4.2.2.
%
% THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
% IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
% THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
% (AT YOUR OPTION) ANY LATER VERSION.
%
% THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL,
% BUT WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
% MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  SEE THE
% GNU GENERAL PUBLIC LICENSE FOR MORE DETAILS.
%
% YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
% ALONG WITH THIS PROGRAM; IF NOT, SEE <HTTP://WWW.GNU.ORG/LICENSES/>.
%
% COPYRIGHT (C) Tim Mitchell, 2019

    A = companion_demo(n);
    
    % Note that the bug still occurs if either:
    % - all entries of A are perturbed by small additive noise, e.g. on the 
    %   order of 1e-12
    % - the entries of A are randomly permuted, e.g. 
    %   R = randperm(n); A = A(R,R);
        
    % Works just fine
    fprintf('Singular values via: s = svd\n');
    s       = svd(A)
    
    % Works just fine on R2017a but not R2018a
    fprintf('Singular value via: [U,S] = svd\n');
    [U,S]   = svd(A);
    s       = diag(S)
    
    % Produces very large garbage for n >= 26 on all tested MATLAB versions
    fprintf('Singular value via: [U,S,V] = svd\n');
    [U,S,V] = svd(A);
    s       = diag(S)  
    
    % The largest singular value seems to be the only one that is computed
    % consistently.  As the other computed singular values are all at least
    % 1e16 times smaller than the largest, checking whether the SVD matches
    % the original matrix to machine precision is essentially pointless.
    
end