function B = jacobi_B( A, p, q )
    u = A(p, q);
    v = .5 * (A(p,p) - A(q,q));
    w = sqrt(u^2 + v^2);
    c = sqrt( (w + abs(v))/(2*w) );
    if v == 0
        signum = 1;
    else
        signum = sign(v);
    end
    s = signum*u/(2*w*c);

    S = eye(size(A));
    S(p,p) = c;
    S(q,q) = c;
    S(q,p) = s;
    S(p,q) = -s;

    B = S * A * S';
end

