% Basic svd test

%for k=1:1000
%  a=randn(100); [u,s,v] = svd(a); svd_std_self_test(a,u,s,v)
%end

% a: m x n
% u: m x m  or  m x p
% s: m x n  or  p x p
% v: n x n  or  n x p

function svd_std_self_test (a, u, s, v, tol)

if ~exist('tol', 'var') || isempty(tol)
  tol = 10; % unit ulp
end

if isempty(a)
  assert (rows(a), rows(u));
  assert (columns(a), rows(v));

  if max(size(s)) == 0
    % 'econ'
    assert (columns (u), 0);
    assert (columns (v), 0);
  else
    % full
    assert (diff (size (u)) == 0);
    assert (diff (size (v)) == 0);
    assert (u, eye (size (u), class (u)));
    assert (v, eye (size (v), class (v)));
  end
  return;
end

m = size (u, 1);
n = size (v, 1);

assert (isdiag (s));
sv = s(1 == eye (size (s)));
assert ( issorted (sv, 'descend') );
assert (~ any (sv(min (m, n)+1:end)) );

ulp = eps(class(s));
anorm = @(x) norm(x, 1);

% See http://www.netlib.org/lapack/explore-html/d1/d3c/group__double__eig_ga154b2caeffc955bb3dede4c310ca3cd6.html#ga154b2caeffc955bb3dede4c310ca3cd6
% All results should be in O(1).
assert ( anorm(a - u * s * v') / (anorm(a) * max(m,n) * ulp) < tol);
assert ( anorm(eye(size(u,2)) - u' * u) / (m * ulp) < tol);
assert ( anorm(eye(size(v,2)) - v' * v) / (n * ulp) < tol);

end
