  function ShowSHiftedSolver()
% The Schur decomposition can be handy in certain linear system settings
n = 300;
m = 2000;
b = randn(n,1);
mu = rand(m,1);
A = randn(n,n); A = A'*A;
I = eye(n,n);
tic
X = ShiftedSolver(A,mu,b);
t1 = toc;
tic
for k=1:m
    X(:,k) = (A+mu(k)*I)\b;
end
t2 = toc;
clc
fprintf('n = %1d\nm = %1d\n',n,m)
fprintf('Time(ShiftedSolver) = %5.2f\n',t1)
fprintf('Time(m Cholesky''s)  = %5.2f\n',t2)





  function X = ShiftedSolver(A,mu,b)
% A is an nxn symmetric positive definite matrix
% mu is a length-m vector of positive values
% (A + mu(k))X(:,k) = b, k=1M

% Idea:
% Q'(A + mu(k)I)Q(Q'x) = Q'b, i.e., (D + mu(k)I)xtilde = btilde 
% where
[Q,D] = eig(A);
[n,n] = size(A);
m = length(mu);
X = zeros(n,m);
% Get all the xtilde's...
btilde = Q'*b;
d = diag(D);
for k=1:m
    X(:,k) = btilde./(d + mu(k));
end
% Since Q'x = xtilde, x = Q*xtilde so
X = Q*X;
