m			= 3;
n_tests  = 100000;

S			= randn([m m n_tests], 'single');
S			= S + permute(S, [2 1 3]);

order		= 'u-';

%% compute with matlab
fprintf('Computing eigenvalues/vectors with matlab ... \n');
tic;
V_matlab = zeros([m m n_tests], 'single');
D_matlab = zeros([m n_tests], 'single');
for n = 1 : n_tests
	[V,D]					= eig(S(:,:,n));
	if strcmp(order, 'u+')
		D					= diag(D);
		[~,idx]			= sort(abs(D), 'ascend');
		D					= D(idx);
	elseif strcmp(order, 'u-')
		D					= diag(D);
		[~,idx]			= sort(abs(D), 'descend');
		D					= D(idx);
	elseif strcmp(order, 's+')
		[D,idx]			= sort(diag(D), 'ascend');
	else
		[D,idx]			= sort(diag(D), 'descend');
	end
	V_matlab(:,:,n)	= V(:, idx);
	D_matlab(:,n)		= D;
end
t_matlab = toc;
fprintf('Took %0.3fs in MATLAB (%0.1f mat/s)\n', t_matlab, n_tests/t_matlab);

%% compute with ndsyev
tic;
[D_syev, V_syev] = ndsyev(S, order);
t_syev	= toc;
fprintf('Took %0.3fs with syev (%0.1f mat/s)\n', t_syev, n_tests/t_syev);

%% compare
D_err			= D_matlab - D_syev;
fprintf('RMS error in eigenvalues: %0.3f\n', rms(D_err,2));
fprintf('RMS eigenvalue:           %0.3f\n', rms(D_matlab,2));

V_err			= zeros(1, n_tests);
for n = 1 : n_tests
	I_err		= V_matlab(:,:,n) * V_syev(:,:,n)' - eye(3);
	V_err(n) = norm(I_err);
end
fprintf('Mean norm of V*V'' - I:   %0.3f\n', mean(V_err));