% test gejsv svd sovler vsign = @(u) sign(sum(u .* abs(u))); maxabs = @(x) max(abs(x(:))); %normalize = @(x) x ./ norm(x, 2, 'cols'); %vsign = @(u) conj(normalize(sum(u .* abs(u)))); function max_norm_err = check_vec(u1, u2) % for checking (complex valued) unitary matrixes fres = @(u, u0) norm(u - (u0 \ u) * u0); norm_err = zeros(1, size(u1, 2)); for k = 1:size(u1, 2) norm_err(k) = fres(u1(:,k), u2(:,k)); end max_norm_err = max(norm_err); end s_ty = {'D', 'S', 'C', 'Z'}; s_is_complex = [0 0 1 1]; typename.D = 'double'; typename.S = 'single'; typename.C = 'single'; typename.Z = 'double'; s_mn = {[4, 4], [4, 3], [3, 4], [1, 5], [5, 1], [100, 200]}; for k_mtype = 1:numel(s_ty) for k_size = 1:numel(s_mn) m = s_mn{k_size}(1); n = s_mn{k_size}(2); p = min(m, n); ty = s_ty{k_mtype}; fprintf('svd: %d x %d ty = %s\n', m, n, ty); if s_is_complex(k_mtype) a = rand(m, n, typename.(ty)) + 1i * rand(m, n, typename.(ty)); else a = rand(m, n, typename.(ty)); end if typename.(ty)(1) == 'd' feps = 1e-11; else feps = 1e-3; end svd_driver('gesvd'); [u1, s1, v1] = svd(a); svd_driver('gejsv'); assert( strcmp(svd_driver(), 'gejsv') == 1 ); [u2, s2, v2] = svd(a); u1 = u1 .* vsign(u1); u2 = u2 .* vsign(u2); v1 = v1 .* vsign(v1); v2 = v2 .* vsign(v2); assert(maxabs(s2 - s1) < feps) assert(check_vec(u2(:,1:p), u1(:,1:p)) < feps) assert(check_vec(v2(:,1:p), v1(:,1:p)) < feps) [u2, s2] = svd(a); u2 = u2 .* vsign(u2); assert(maxabs(s2 - s1) < feps) assert(check_vec(u2(:,1:p), u1(:,1:p)) < feps) [~, s2, v2] = svd(a); v2 = v2 .* vsign(v2); assert(maxabs(s2 - s1) < feps) assert(check_vec(v2(:,1:p), v1(:,1:p)) < feps) s2 = svd(a); assert(maxabs(s2 - s1(eye(size(s1))==1)) < feps) end end