% Test accuracy of svd() by a hard problem
%
% With svd_driver('gesvd')
% svd(A) is good upto N==14, with max relative error 1.4e-6.
% at N==19, svd(A) max relative error 0.34, start to return wrong answer.
% MATLAB behaviours exactly the same as 'gesdd'.
%sumsq = @(x) sum(x.*x);  % for MATLAB^TM

svd_driver('gesvd')
%svd_driver('gesdd')

N = 26
A = compan (fliplr ([1, 1 ./ cumprod(1:N)]));

[u, compan_svd, v] = svd(A);
compan_svd = diag(compan_svd);

a = A(1,:);
% Characteristic polynomial of A*A', note svd(A)=eig(A*A')=roots(c).
c = [1, -(1-sumsq(a)), -sumsq(a(1:end-1)), zeros(1, length(a)-2)];

% Numerical stable answer.
compan_svd_formula = sort(1-[roots(c(1:3)); zeros(length(a)-2,1)], 'descend') .^0.5;

% Solve svd(A) by eig(A*A').
compan_svd_eig = sort(eig(A*A'), 'descend') .^0.5;

disp('[compan_svd, compan_svd_eig, compan_svd_formula]');
fprintf('%14.6g, %14.6g, %14.6g\n', ...
   [compan_svd, compan_svd_eig, compan_svd_formula]');

max_relative_error = max(abs(compan_svd - compan_svd_formula)./compan_svd_formula)

% Running result
% Octave-4.0.3 (lapack-3.7.0, openblas-0.2.19)
%{
N =  26
[compan_svd, compan_svd_eig, compan_svd_formula]
   6.08901e+26,    6.08901e+26,    6.08901e+26
   3.76175e+09,        1.03585,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
             1,              1,              1
      0.662326,       0.956064,              1
   2.65834e-10,    8.48585e-06,       0.662326
max_relative_error =    3.7617e+09
%}
