% eigen3r.m  test for finding A given e and v or given 1/e and v
%            the eigenvalues of inv(A) are 1/eigenvalues of A
%            the eigenvectors of inv(A) are the eigenvectors of A
%            the eigenvalues may differ in order, eigenvectors sign may flip   
function eigen3r
 diary('eigen3r_m.out')
 format compact
 sprintf('random matrix')
 A=rand(3,3); %real
 B=rand(3,3); %imag
 for j=1:3
   for k=1:3
     A(j,k)=A(j,k)+(B(j,k)*1.0i);
   end
 end
 A = A
 sprintf('determinant of complex A')
 da=det(A)
 [v e]=eig(A) % compute eigenvectors and eigenvalues
 sprintf('av1 = A*v(i) should equal  evi = e(i)*v(i) known')
 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)
 
 I=eye(3);    % identity matrix
 sprintf('det(A-Ie) error')  
 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
 sprintf('trace(A)-trace(e) = ')
 z3=trace(A)-trace(e) % should be near zero
 sprintf('det(A-Ie) = ')
 for i=1:3
   z3=det(A-I.*e(i,i)) % should be near zero
 end
 sprintf('det(e) =')  
 z4=det(e) % should be zero
 sprintf('inverse of A')
 AI=inv(A)
 sprintf('det(inverse of A)')
 dai=det(AI)
 [vi ei]=eig(AI) % compute eigenvectors and eigenvalues
 sprintf('check ei = 1/e')
 e1inv = 1.0/e(1,1)
 e2inv = 1.0/e(2,2)
 e3inv = 1.0/e(3,3)
	 
 sprintf('av1 = AI*vi(i) should equal  evi = ei(i)*vi(i) known')
 av1 = AI*vi(:,1)
 ev1 = ei(1,1)*vi(:,1)
 av2 = AI*vi(:,2)
 ev2 = ei(2,2)*vi(:,2)
 av3 = AI*vi(:,3)
 ev3 = ei(3,3)*vi(:,3)

 diary off
end
