rand('state',31415926535);   % initialize random number generator -> reproducible results
V = rand(3);                 % create random eigenvector matrix
L = diag([3,3,1]);           % eigenvalues: {3, 3, 1}
A = V*L*inv(V);              % A has eigenvectors from V, eigenvalues from L

fprintf(1,'eigenvalues: %23.15e\n', diag(L));

v = zeros(3,21);
v(1:3,1) = [1;1;1];          % starting vector
l(1) = norm(v(:,1),'inf'); 
v(:,1) = v(:,1)/l(1);        % normalize

for k = 2:20,
  v(:,k) = A*v(:,k-1);
  l(k) = norm(v(:,k),'inf');
  v(:,k) = v(:,k)/l(k);
  fprintf(1,'iteration = %3d, approx. abs. eig. val. = %23.15e\n', k-1, l(k));
  pause
end

k = 21;
v(:,k) = A*v(:,k-1);
[maxv,idx] = max(abs(v(:,21)));
lambda = v(idx,21)/v(idx,20);
keyboard
fprintf(1,'approx. eig. val. = %23.15e\n', lambda);
fprintf(1,'norm(A*v(:,20)-lambda*v(:,20)) = %23.15e\n', ...
  norm(v(:,21)-lambda*v(:,20)));
fprintf(1,'give eigenvector relative to eigenvector basis V(:,1), V(:,2), V(:,3)\n');
fprintf(1,'%23.15e\n',(V'*V)\(V'*v(:,21)) );

keyboard

L = diag([3+4*i, 3-4*i,1]);         % eigenvalues: {3+4i, 3-4i, 1}
A = V*L*inv(V);                     % A has eigenvectors from V, eigenvalues from L

fprintf(1,'eigenvalues:\n');
fprintf(1,'(%23.15e , %23.15e)\n', real(L(1,1)), imag(L(1,1)) );
fprintf(1,'(%23.15e , %23.15e)\n', real(L(2,2)), imag(L(2,2)) );
fprintf(1,'(%23.15e , %23.15e)\n', real(L(3,3)), imag(L(3,3)) );

v = zeros(3,21);
v(1:3,1) = [1;1;1];
l(1) = norm(v(:,1),'inf'); 
v(:,1) = v(:,1)/l(1);               % normalize

for k = 2:20,
  v(:,k) = A*v(:,k-1);
  l(k) = norm(v(:,k),'inf');
  v(:,k) = v(:,k)/l(k);
  fprintf(1,'iteration = %3d, approx. abs. eig. val. = %23.15e\n', k-1, l(k));
  pause
end

k = 21;
v(:,k) = A*v(:,k-1);
[maxv,idx] = max(abs(v(:,21)));
lambda = v(idx,21)/v(idx,20);
keyboard
fprintf(1,'approx. eig. val. = %23.15e\n', lambda);
fprintf(1,'norm(A*v(:,20)-lambda*v(:,20)) = %23.15e\n', ...
  norm(v(:,21)-lambda*v(:,20)));
fprintf(1,'give eigenvector relative to eigenvector basis V(:,1), V(:,2), V(:,3)\n');
fprintf(1,'%23.15e\n',(V'*V)\(V'*v(:,21)) );