function svdzero(n)
% Tim Mitchell
% mitchell@mpi-magdeburg.mpg.de
% Date: February 13, 2019
%
% svdzero(n): 
% This demo shows that svd can sometimes return -0 as a singular value,
% when it should return +0 instead.  The demo creates an upper triangular
% randn matrix of order n and then sets the top left entry to zero.  It
% then computes the minimum singular value to show that when singular
% vectors are also requested, -0 can sometimes be returned.  This appears
% to be an issue in GESVD.  On R2017b, even 1/svd(-0) will return -inf (but
% not R2018a).
% 
% If n <= 25, MATLAB will only call GESVD for all tests. 
% If n > 25, MATLAB will instead call GESDD when the following singular 
% vectors are requested:
%   - right AND left, i.e. [U,S,V] = svd, on R2017b and earlier
%   - right OR left, i.e. [U,S] = svd or [U,S,V] = svd, on R2018a and later
%
% This behavior has also been confirmed on Octave for [U,S,V] = svd calls.
% Note that Octave does not support [U,S] = svd calls so if you are running
% this test on Octave, you will need to comment out the two [U,S] = svd 
% calls below, along with their associated print statements.
%
% 
% 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 = randn(n);
    A = triu(A);
    
    A(1) = 0;
    fprintf('\nA = randn(n); A = triu(A); A(1) = 0\n')
    fprintf('Reciprocal of smallest singular value:\n');
    s = svd(A);
    fprintf('s = svd(A):       %g\n',s(end));
    [U,S] = svd(A);
    fprintf('[U,S] = svd(A):   %g\n',S(end));
    [U,S,V] = svd(A);
    fprintf('[U,S,V] = svd(A): %g\n\n',S(end));

    A(1) = -0;
    fprintf('Same matrix but now A(1) = -0\n');
    fprintf('Reciprocal of smallest singular value:\n');
    s = svd(A);
    fprintf('s = svd(A):       %g\n',s(end));
    [U,S] = svd(A);
    fprintf('[U,S] = svd(A):   %g\n',S(end));
    [U,S,V] = svd(A);
    fprintf('[U,S,V] = svd(A): %g\n\n',S(end));

end