function x = SemiSepLsolve(d,p,q,b)
% d, p, q, and b are column n-vectors
% The matrix L = diag(d) + tril(p*q') is nonsingular
% x solve Lx = b

n = length(b);
nmin = 2;
if n<=nmin
    % Solve conventionally...
    x = (diag(d) + tril(p*q'))\b;
    %x = b/(d(1) + p(1)*q(1));
else
    % Divide and conquer...
    m = floor(n/2);
    top = 1:m;
    bot = m+1:n;
    % Solve the L11 system recursively...
    xT = SemiSepLsolve(d(top),p(top),q(top),b(top));
    % Update the bottom rhs...
    btilde = b(bot) - (q(top)'*xT)*p(bot);
    % Solve the L22 system recursively...
    xB = SemiSepLsolve(d(bot),p(bot),q(bot),btilde);
    x = [xT;xB];
end

% -3 if you form 2d arrays