function [M,J] = jordan(A,epp,str)
%
%  Function JORDAN finds a Normal Jordan Form of the matrix A,
%  A = M*J*inv(M). If A can be ....
%
% Does not always give you the Normal Jordan Form.
%
%

   if nargin < 3
      str = '';
   end
   if nargin < 2
      epp = 1e-12;
   end
   
   if (strcmp(str, 'balance') == 1) | (strcmp(str,'') == 1)
       str = 'nobalance';
       [Dsc,B] = balance(A);
       [M,J] = jordan(B,epp,str);
       M = Dsc*M; J = J;
       return;
   end
      
   
   [U1,T1] = schur(A);
   T1 = zerm(T1,1e-12);
   [U2,T2] = rsf2csf(U1,T1);

   P = diag(T2);
   Nx = max(size(P));
   PR = sorte(P);
   Pinfo = []; PD = [];
   
   while isempty(PR) == 0
      ipm = find( abs(PR-PR(1)) < epp );
      mp = max(size(ipm));
      Pinfo = [Pinfo [PR(1); mp]];
      for i = 1:mp
         Ph(i,1) = PR(1);
      end
      PR = rmvecel(PR,ipm);
      PD = [PD; Ph];
   end
   Pinfo = Pinfo
   %
   %  Omrokker egenverdiene. 
   %
   PD = sorte(PD);
   [U,T] = qrsort(U2,T2,PD);
   T = zerm(T,1e-12);
   sjnrm = norm(zerm(U*T*U'-A));
   if sjnrm > epp*10
      disp('Warning: reordered schur decomposition is inaccurate');
      [U,T] = reordschr(U2,T2,PD,epp);
      sjnrmII = norm(zerm(U*T*U'-A));
      if sjnrmII > epp*10
         disp('Warning: reordered schur decomposition is inaccurate');
      else
         disp('Corrected: reordered schur decomposition is OK');
      end
   end
   % Now we are ready to start computing eigenvalue problems.
   % Consentrates on the right eigenvalue problem.
   
   if str ~= ''
      [VRh,DR] = eig(T,str);                   PR = diag(DR);
%      % Solve the left eigenvalue problem.
%      [VLh,DL] = eig(T',str); VLh = conj(VLh); PL = diag(DL);
   else
      [VRh,DR] = eig(T);     PR = diag(DR);
%      % Solve the left eigenvalue problem.
%      [VLh,DL] = eig(T'); VLh = conj(VLh); PL = diag(DL);  
   end


   NP = size(Pinfo,2);
   if NP == Nx
      M = U*VRh;
      J = DR;
      return
   end
   
   VR = []; ON = 1; 
   ERRmat = [];
   for i=1:NP
      Rinz = find( abs( Pinfo(1,i)-PR ) < epp );
      
      if Pinfo(2,i) > 1
         EV = zerm(VRh(:,Rinz(:,1)));
	 IDPEV  = rank(EV, epp);
	 
	 if IDPEV == Pinfo(2,i)         % Independent eigenvectors.
	    M = [M zerm(EV)];
	    
	 else
	    MZ = T-Pinfo(1,i)*eye(Nx);
	 
	    VNUM = []; MH = [];
	    REMVR = Pinfo(2,i) - IDPEV + 1; 
	    for j = 1:IDPEV
	       if j > 1
	          REMVR = Pinfo(2,i)-IDPEV;
	       end
	       [sev,sid] = fidpev(M,EV);
	       if isempty(sev)
	          disp('Warning: Cannot find the missing eigenvector.');
		  ERRmat = [ERRmat [Pinfo(1,i);Pinfo(2,i); REMVR; size(M,2)+1]]; 
	       else
	          EV = rmmatcd(EV,sid);
	          VR = itprjr(MZ,sev,REMVR,epp);
	          VNUM(j) = size(VR,2);
	          MH = [MH zerm(VR)];
	       end
	    end

            if size(MH,2) > Pinfo(2,i)	
               MHV = []; nstr = 1;
               for blkc=1:IDPEV 
                  MTMP = jorana(MH(:,nstr:nstr+VNUM(blkc)-1),T,Pinfo(1,i),epp);
                  MHV = [MHV MTMP];
                  nstr = nstr+VNUM(blkc);
               end
               M = [M MHV];
            elseif size(MH,2) == Pinfo(2,i)
               M = [M MH];
            else
	       disp('Warning: to few generalized vectors.')
            end

	 end
      else
         M = [M VRh(:,Rinz(:,1))];
      end
   end
   
   M = zerm(U*M);
   if isempty(ERRmat) == 0
      for i = 1:size(ERRmat,2);
         ERRmat = ERRmat;
         [EVA,DA] = eig(A); PDA = diag(DA);
         Rinz = find( abs( ERRmat(1,i)-PDA ) < epp );
      
         EV = zerm(EVA(:,Rinz(:,1)));
         ML = M(:, 1:ERRmat(4,i)-1);
  	 MR = M(:,ERRmat(4,i):size(M,2));
         [sev,sid] = fidpev(M,EV);
	 MZ = A-eye(Nx)*ERRmat(1,i);
	 VR = itprjr(MZ,sev,ERRmat(3,i));
	 nvr = size(VR,2);
	 M = zerm([ML,VR,MR]);
	 for k=i+1:size(ERRmat,2)
	    ERRmat(4,k) = ERRmat(4,k)+nvr;
	 end
      end
   end
   
   RKM = rank(M);
   if RKM < Nx
      disp(['Error: could not find ',num2str(Nx), ' independent generalized eigenvectors']);
      J = zerm(pinv(M)*A*M);
   else
      J = zerm(inv(M)*A*M);
   end   
end





