% 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