function [alpha,beta] = MakeTridiag(d,q)
% d is a column n-vector with distinct entries.
% q is a column n-vector with unit 2-norm and positive entries.
% alpha is a column n-vector and beta is a column (n-1)-vector with the property that
% there is an orthogonal Q with Q(:,n) = q so that T(alpha,beta) =
% Q'*diag(d)*Q.

% Idea: Run Lanczos with A = diag(d1,...,dn)...
n = length(d); alpha = zeros(n,1); beta = zeros(n-1,1);
Q = zeros(n,n); Q(:,1) = q;
for k=1:n
y = d.*Q(:,k);
alpha(k) = Q(:,k)'*y;
if k>1
rk = y - alpha(k)*Q(:,k) - beta(k-1)*Q(:,k-1);
else
rk = y - alpha(k)*Q(:,k);
end
if k<n
% Project out the unwanted stuff...
rk = rk - Q(:,1:k)*(Q(:,1:k)'*rk);
beta(k) = norm(rk);
Q(:,k+1) = rk/beta(k);
end
end
% We now have Q'*diag(d)*Q = T but Q(:,n) = q.
% Let E = I(:,n:-1:1). We then have (QE)'*diag(d)*(QE) = E'T*E = T1
% Note that q is the last column of QE and that T1 is tridiagonal.
alpha = alpha(n:-1:1);
beta = beta(n-1:-1:1);