% eigen10.m  demonstrate more eigenvalue testing
function eigen10
  diary('eigen10_m.out')
  format compact
  disp('eigen10.m running')
  disp('A =')
  A=zeros(5);
  A(1,1)= 1.0; % build matrix A
  A(1,2)= 1.0;
  A(1,3)= 1.0;
  A(1,4)= 1.0;
  A(1,5)= 0.0;
  A(2,1)=-1.0;
  A(2,2)= 1.0;
  A(2,3)= 1.0;
  A(2,4)= 1.0;
  A(2,5)= 1.0;
  A(3,1)= 0.0;
  A(3,2)=-1.0;
  A(3,3)= 1.0;
  A(3,4)= 1.0;
  A(3,5)= 1.0;
  A(4,1)= 0.0;
  A(4,2)= 0.0;
  A(4,3)=-1.0;
  A(4,4)= 1.0;
  A(4,5)= 1.0;
  A(5,1)= 0.0;
  A(5,2)= 0.0;
  A(5,3)= 0.0;
  A(5,4)=-1.0;
  A(5,5)= 1.0;
  disp(A)
  disp(' ')
  AT = transpose(A); % should have same eigenvalues
  disp('AT =')
  disp(AT)
  disp(' ')
  disp('check A')
  [v e]=eig(A) % compute eigenvectors and eigenvalues

  disp('A-eI=0')
  I=eye(5);    % identity matrix
  for i=1:5
    z1=det(A-I.*e(i,i)) % should be near zero
  end
  disp('A*v-e*v=0')
  for i=1:5
    z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(A)=trace(e)=0')
  z3=trace(A)-trace(e) % should be near zero

  disp('orthogonal col')
  for i=1:5
    disp(i)
    for j=i+1:5
      disp(j)
      dot=0.0+0j;
      for k=1:5
	dot = dot + v(k,i)*v(k,j);
      end
      disp(dot)
    end
  end

  disp('orthogonal row')
  for i=1:5
    disp(i)
    for j=i+1:5
      disp(j)
      dot=0.0+0j;
      for k=1:5
	dot = dot + v(i,k)*v(j,k);
      end
      disp(dot)
    end
  end

  disp(' ')
  disp('check AT')
  disp('AT =')
  disp(AT)
  [v e]=eig(AT) % compute eigenvectors and eigenvalues
  disp('AT-eI=0')
  I=eye(5);    % identity matrix
  for i=1:5
    z1=det(AT-I.*e(i,i)) % should be near zero
  end
  disp('AT*v-e*v=0')
  for i=1:5
    z2=AT*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(AT)=trace(e)=0')
  z3=trace(AT)-trace(e) % should be near zero

  disp(' ')
  disp(' ')
  disp('B =')
  B=zeros(4);
  B(1,1)= 4.0; % build matrix B
  B(1,2)= 4.0;
  B(1,3)= 4.0;
  B(1,4)= 4.0;
  B(2,1)= 4.0;
  B(2,2)= 4.0;
  B(2,3)= 4.0;
  B(2,4)= 4.0;
  B(3,1)= 4.0;
  B(3,2)= 4.0;
  B(3,3)= 4.0;
  B(3,4)= 4.0;
  B(4,1)= 4.0;
  B(4,2)= 4.0;
  B(4,3)= 4.0;
  B(4,4)= 4.0;
  disp(B)
  disp(' ')
  disp('check B')
  [v e]=eig(B) % compute eigenvectors and eigenvalues

  disp('B-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(B-I.*e(i,i)) % should be near zero
  end
  disp('B*v-e*v=0')
  for i=1:4
    z2=B*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(B)=trace(e)=0')
  z3=trace(B)-trace(e) % should be near zero
  disp(' ')
  % BI = inv(B); % no singular

  disp(' ')
  disp(' ')
  disp('C =')
  C=zeros(4);
  C(1,1)= 1.0; % build matrix C
  C(2,2)= 2.0;
  C(3,3)= 3.0;
  C(4,4)= 4.0;
  disp(C)
  disp(' ')
  disp('check C')
  [v e]=eig(C) % compute eigenvectors and eigenvalues

  disp('C-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(C-I.*e(i,i)) % should be near zero
  end
  disp('C*v-e*v=0')
  for i=1:4
    z2=C*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(C)=trace(e)=0')
  z3=trace(C)-trace(e) % should be near zero
  disp(' ')
  CI = inv(C); % should have recipricol eigenvalues
  disp('CI =')
  disp(CI)
  disp(' ')
  disp('check CI')
  [v e]=eig(CI) % compute eigenvectors and eigenvalues

  disp('CI-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(CI-I.*e(i,i)) % should be near zero
  end
  disp('CI*v-e*v=0')
  for i=1:4
    z2=CI*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(CI)=trace(e)=0')
  z3=trace(CI)-trace(e) % should be near zero
  disp(' ')

  disp(' ')
  disp(' ')
  disp('D =')
  D=zeros(4);
  D(1,1)= 0.5; % build matrix D
  D(1,2)= 0.333333;
  D(1,3)= 0.25;
  D(1,4)= 0.2;
  D(2,1)= 0.333333;
  D(2,2)= 0.25;
  D(2,3)= 0.2;
  D(2,4)= 0.166666;
  D(3,1)= 0.25;
  D(3,2)= 0.2;
  D(3,3)= 0.166666;
  D(3,4)= 0.142857;
  D(4,1)= 0.2;
  D(4,2)= 0.166666;
  D(4,3)= 0.142857;
  D(4,4)= 0.125;
  disp(D)
  disp(' ')
  disp('check D')
  [v e]=eig(D) % compute eigenvectors and eigenvalues

  disp('D-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(D-I.*e(i,i)) % should be near zero
  end
  disp('D*v-e*v=0')
  for i=1:4
    z2=D*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(D)=trace(e)=0')
  z3=trace(D)-trace(e) % should be near zero
  disp(' ')
  DI = inv(D); % should have recipricol eigenvalues
  disp('DI =')
  disp(DI)
  disp(' ')
  disp('check DI')
  [v e]=eig(DI) % compute eigenvectors and eigenvalues

  disp('DI-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(DI-I.*e(i,i)) % should be near zero
  end
  disp('DI*v-e*v=0')
  for i=1:4
    z2=DI*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(DI)=trace(e)=0')
  z3=trace(DI)-trace(e) % should be near zero
  disp(' ')

  disp(' ')
  disp(' ')
  disp('E =')
  E=zeros(4);
  E(1,1)= 2.0; % build matrix D
  E(1,2)= 3.0;
  E(1,3)= 4.0;
  E(1,4)= 5.0;
  E(2,1)= 3.0;
  E(2,2)= 4.0;
  E(2,3)= 5.0;
  E(2,4)= 6.0;
  E(3,1)= 4.0;
  E(3,2)= 5.0;
  E(3,3)= 6.0;
  E(3,4)= 7.0;
  E(4,1)= 5.0;
  E(4,2)= 6.0;
  E(4,3)= 7.0;
  E(4,4)= 8.0;
  disp(E)
  disp(' ')
  disp('check E')
  [v e]=eig(E) % compute eigenvectors and eigenvalues

  disp('E-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(E-I.*e(i,i)) % should be near zero
  end
  disp('E*v-e*v=0')
  for i=1:4
    z2=E*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp('trace(E)=trace(e)=0')
  z3=trace(E)-trace(e) % should be near zero
  disp(' ')
  % EI = inv(E); % E is singular

  disp(' ')
  disp(' ')
  disp('F =')
  F=zeros(4);
  F(1,1)= 1.0/2.0; % build matrix D
  F(1,2)= 1.0/3.0;
  F(1,3)= 1.0/4.0;
  F(1,4)= 1.0/5.0;
  F(2,1)= 1.0/3.0;
  F(2,2)= 1.0/4.0;
  F(2,3)= 1.0/5.0;
  F(2,4)= 1.0/6.0;
  F(3,1)= 1.0/4.0;
  F(3,2)= 1.0/5.0;
  F(3,3)= 1.0/6.0;
  F(3,4)= 1.0/7.0;
  F(4,1)= 1.0/5.0;
  F(4,2)= 1.0/6.0;
  F(4,3)= 1.0/7.0;
  F(4,4)= 1.0/8.0;
  disp(F)
  disp(' ')
  disp('check F')
  [v e]=eig(F) % compute eigenvectors and eigenvalues

  disp('F-eI=0')
  I=eye(4);    % identity matrix
  for i=1:4
    z1=det(F-I.*e(i,i)) % should be near zero
  end
  disp('F*v-e*v=0')
  for i=1:4
    z2=F*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
  end
  disp(' ')
  disp('trace(F)=trace(e)=0')
  z3=trace(F)-trace(e) % should be near zero
  disp(' ')

  disp('vector length=')
  for i=1:4
    len = 0.0;
    for j=1:4
    len = len + v(i,j);
    end
    disp(sqrt(len))
  end
  disp(' ')
	    
  disp('orthogonal col')
  for i=1:3
    disp(i)
    for j=i+1:4
      dot=0.0+0j;
      for k=1:4
	dot = dot + v(k,i)*v(k,j);
      end
      disp('dot =')
      disp(dot)
    end
  end
  disp(' ')
	
  diary off
end
