function M = get_spd(n)
% M = get_spd(n)
%
% Make a random n-by-n SPD matrix. The
% eigenvalues of M are
%   logspace(0, 1, n)
%
% In:  n : integer > 0
% Out: M : random n-by-n SPD matrix.

M = diag( logspace(0, 1, n) );

r = 2 * pi * rand(n,1);
for i = 2 : n
   c = cos(r(i));
   s = sin(r(i));
   R = [c, s; -s, c];
   M([i,1], :) = R' * M([i,1], :);
   M(:, [i,1]) = M(:, [i,1]) * R;
end

r = 2 * pi * rand(n,1);
for i = 2 : n
   j = i + ceil( (n-i) * rand );
   c = cos(r(i));
   s = sin(r(i));
   R = [c, s; -s, c];
   M([i,j], :) = R' * M([i,j], :);
   M(:, [i,j]) = M(:, [i,j]) * R;
end
