function B = sqrtm_schurparlett(A)
[Q, U] = schur(A, 'complex');
n = size(A, 1);
S = zeros(n, n);
for j = 1:n
    S(j,j) = sqrt(U(j,j));
    for i = j-1:-1:1
        num = U(i,j) - S(i,i+1:j-1)*S(i+1:j-1,j);
        S(i,j) = num / (S(i,i)+S(j,j));
    end
end
B = Q*S*Q';
if isreal(A)
    B = real(B);
end