% eigen9.m  demonstrate eigenvalues in MatLab special cases
function eigen5
  diary('eigen9_m.out')
  format compact
  n = 1
  A=zeros(n);
  A(1,1)=1.0; % build matrix A
  disp('A=')
  disp(A)
  disp('det(A)=')
  disp(det(A))

  [v e]=eig(A) % compute eigenvectors and eigenvalues
  I=eye(n);    % identity matrix

  disp('length of each v should be 1')
  for i=1:n
    l = 0.0;
    for j=1:n
      l = l + v(j,i)*v(j,i);
    end
    l = sqrt(l)
  end
  disp(' ')

  disp('det(A-eI) should be near zero')  
  for i=1:n
    z1=det(A-I.*e(i,i)) % should be near zero
  end
  disp(' ')

  disp('A*v-e*v should be near zero')
  for i=1:n
    z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp(' ')
  disp(' ')

  disp('trace(A)-trace(e) should be near zero')
  z3=trace(A)-trace(e) % should be near zero
  disp(' ')

  n = 3
  A(1,1)=1.0; % build matrix A 
  A(1,2)=1.0; %
  A(1,3)=1.0; % | 1.0  1.0  1.0 |
  A(2,1)=1.0; % | 1.0  1.0  1.0 |
  A(2,2)=1.0; % | 1.0  1.0  1.0 |
  A(2,3)=1.0;
  A(3,1)=1.0;
  A(3,2)=1.0;
  A(3,3)=1.0;
  disp('A=')
  disp(A)
  disp('det(A)=')
  disp(det(A))

  [v e]=eig(A) % compute eigenvectors and eigenvalues
  I=eye(n);    % identity matrix

  disp('length of each v should be 1')
  for i=1:n
    l = 0.0;
    for j=1:n
      l = l + v(j,i)*v(j,i);
    end
    l = sqrt(l)
  end
  disp(' ')

  disp('det(A-eI) should be near zero')  
  for i=1:n
    z1=det(A-I.*e(i,i)) % should be near zero
  end
  disp(' ')

  disp('A*v-e*v should be near zero')
  for i=1:n
    z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp(' ')

  disp('trace(A)-trace(e) should be near zero')
  z3=trace(A)-trace(e) % should be near zero
  disp(' ')
  disp(' ')


  n = 4
  A(1,1)=2.0; % build matrix A 
  A(1,2)=2.0; %
  A(1,3)=2.0; % | 2.0  2.0  2.0  2.0 |
  A(1,4)=2.0; % | 2.0  2.0  2.0  2.0 |
  A(2,1)=2.0; % | 2.0  2.0  2.0  2.0 |
  A(2,2)=2.0; % | 2.0  2.0  2.0  2.0 |
  A(2,3)=2.0;
  A(2,4)=2.0;
  A(3,1)=2.0;
  A(3,2)=2.0;
  A(3,3)=2.0;
  A(3,4)=2.0;
  A(4,1)=2.0;
  A(4,2)=2.0;
  A(4,3)=2.0;
  A(4,4)=2.0;
  disp('A=')
  disp(A)
  disp('det(A)=')
  disp(det(A))

  [v e]=eig(A) % compute eigenvectors and eigenvalues
  I=eye(n);    % identity matrix

  disp('length of each v should be 1')
  for i=1:n
    l = 0.0;
    for j=1:n
      l = l + v(j,i)*v(j,i);
    end
    l = sqrt(l)
  end
  disp(' ')

  disp('det(A-eI) should be near zero')  
  for i=1:n
    z1=det(A-I.*e(i,i)) % should be near zero
  end
  disp(' ')

  disp('A*v-e*v should be near zero')
  for i=1:n
    z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp(' ')

  disp('trace(A)-trace(e) should be near zero')
  z3=trace(A)-trace(e) % should be near zero
  disp(' ')

  n = 4
  A(1,1)=1.0/(1+1); % build matrix A  hermitian
  A(1,2)=1.0/(1+2); %
  A(1,3)=1.0/(1+3); % | 0.5    0.333  0.25  0.2   |
  A(1,4)=1.0/(1+4); % | 0.333  0.25   0.2   0.166 |
  A(2,1)=1.0/(2+1); % | 0.25   0.2    0.166 0.143 |
  A(2,2)=1.0/(2+2); % | 0.2    0.166  0.143 0.125 |
  A(2,3)=1.0/(2+3);
  A(2,4)=1.0/(2+4);
  A(3,1)=1.0/(3+1);
  A(3,2)=1.0/(3+2);
  A(3,3)=1.0/(3+3);
  A(3,4)=1.0/(3+4);
  A(4,1)=1.0/(4+1);
  A(4,2)=1.0/(4+2);
  A(4,3)=1.0/(4+3);
  A(4,4)=1.0/(4+4);
  disp('A=')
  disp(A)
  disp('det(A)=')
  disp(det(A))

  [v e]=eig(A) % compute eigenvectors and eigenvalues
  I=eye(n);    % identity matrix

  disp('length of each v should be 1')
  for i=1:n
    l = 0.0;
    for j=1:n
      l = l + v(j,i)*v(j,i);
    end
    l = sqrt(l)
  end
  disp(' ')

  disp('det(A-eI) should be near zero')  
  for i=1:n
    z1=det(A-I.*e(i,i)) % should be near zero
  end
  disp(' ')

  disp('A*v-e*v should be near zero')
  for i=1:n
    z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp(' ')

  disp('trace(A)-trace(e) should be near zero')
  z3=trace(A)-trace(e) % should be near zero
  disp(' ')

  Y=[10, 20, 30, 40]
  X=[4, 3, 2, 1]
  YY=A*X'
  z4=YY-Y'  % should be near zero
  ap=pinv(A)
  X=ap*Y'
  YYY=A*X % should be Y
  z5=YYY-Y' % should be near zero
  diary off
end
