% 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

