%  [D, S, VL,VR] = eigpr(A,epp,str);
%
%  Function to solve both left and right eigenvectors simultanoulsly.
%  Based on the use of MATLABS EIG function.
%  Solves:          A*VR = VR*D       and    VL'*A = D*VL';
%  Scaling:     VR(:,i)'*VR(:,i) = 1  and  VL(:,i)'*VL(:,i) = 1;
%  In addition: S(i,i) = VL(:,i)'*VR(:,i);
%  Distinct eigenvalues:
%     VR*inv(S)*D*VL' = A  and  VR*D*inv(S)*VL' = A
%
%  Output arguments:
%     D   - Diagonal matrix containing the eigenvalues.
%     S   - Diagonal matrix containing the scalings as described above.
%     VL  - Left eigenvectors.
%     VR  - Right eigenvectors.
%  Input arguments:
%     A   - Matrix to find eigenvalues and eigenvectors.
%     epp - Tolerance, default is 1e-12.
%     str - 'balance' or 'nobalance', default is 'nobalance'.
%
%  Written by: Kjetil Havre, 12/9-1996.
%

function [D,S,VL,VR] = eigpr(A,epp,str)
   D = []; S=[]; VL =[]; VR=[];
   if nargin < 3
      str = '';
   end
   if nargin < 2
      epp = 1e-12;
   end
   
   % Solve the standard eigenvalue problem.
   if str ~= ''
      [VRh,DR] = eig(A,str);                   PR = diag(DR);
      % Solve the left eigenvalue problem.
      [VLh,DL] = eig(A',str); VLh = conj(VLh); PL = diag(DL);
   else
      [VRh,DR] = eig(A);     PR = diag(DR);
      % Solve the left eigenvalue problem.
      [VLh,DL] = eig(A'); VLh = conj(VLh); PL = diag(DL);  
   end
   Nx = max(size(diag(DR))); S = zeros(Nx,Nx);
   D = zeros(Nx,Nx);
   
   VL = []; VR = [];
   i = 1; 
   while isempty(PR) == 0
      Linz = find( abs( PL-PR(1) ) < epp );
      Rinz = find( abs( PR-PR(1) ) < epp );
      ML   = max(size(Linz));
      MR   = max(size(Rinz));
      if ML ~= MR
         disp('Error: multiplicity of left and right eigenvalues does not agree');
	 return;
      end
	       
      if MR > 1
         [U,S,V] = svd(A-PR(1)*eye(Nx));
	 isvz = find(diag(S)<epp);
	 nisvz = max(size(isvz));
	 for j=1:MR
	    D(i+j-1,i+j-1) = PR(1);
	    VR = [VR VRh(:,Rinz(j))];
	    VL = [VL VLh(:,Linz(j))];
    	    S(i+j-1,i+j-1)  = VL(:,i+j-1)'*VR(:,i+j-1); 
         end
      elseif MR == 1
         D(i,i) = PR(1);
	 VR = [VR VRh(:,Rinz(j))];
	 VL = [VL VLh(:,Linz(j))]; 
      else
         disp('Should never occure');
	 return;
      end
	 
      PR = rmvecel(PR,Rinz);
      PL = rmvecel(PL,Linz);
      VRh = rmmatcd(VRh,Rinz);
      VLh = rmmatcd(VLh,Linz);
      i = i + MR;
   end
end
