% eigen3q.m  test for finding A given e and v or 1/e and v
%            the eigenvalues of inv(A) are 1/eigenvalues of A
%            the eigenvectors of inv(A) are eigenvectors of A
%            the computed order may be different and possibly sign change  
function eigen3q
 diary('eigen3q_m.out')
 format compact
 e1=1; % desired eigenvalues, complex
 e2=2;
 e3=3;
 P=poly([e1 e2 e3]) % make polynomial from roots
 r=roots(P) % check that roots come back
 A=zeros(3);
 A(1,1)=-P(2); % build matrix A
 A(1,2)=-P(3);
 A(1,3)=-P(4);
 A(2,1)=1;
 A(3,2)=1;
 A
 da=det(A)
 [v e]=eig(A) % compute eigenvectors and eigenvalues
 av1 = A*v(:,1)
 ev1 = e(1,1)*v(:,1)
 av2 = A*v(:,2)
 ev2 = e(2,2)*v(:,2)
 av3 = A*v(:,3)
 ev3 = e(3,3)*v(:,3)
 ee=e % expect e1, e2, e3
 I=eye(3);    % identity matrix
 for i=1:3
   z1=det(A-I.*e(i,i)) % should be near zero
 end
 for i=1:3
   z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
 end
 z3=trace(A)-trace(e) % should be near zero
 for i=1:3
   z3=det(A-I.*ee(i,i)) % should be near zero, not in MatLab
 end
 z4=det(e) % should be zero
 AI=inv(A)
 dai=det(AI)
 [v e]=eig(AI) % compute eigenvectors and eigenvalues
 av1 = AI*v(:,1)
 ev1 = e(1,1)*v(:,1)
 av2 = AI*v(:,2)
 ev2 = e(2,2)*v(:,2)
 av3 = AI*v(:,3)
 ev3 = e(3,3)*v(:,3)

 diary off
end
