function [V]=rksm(A,res,m,tol,emin,emax);
%function [V]=rksm(A,res,m,tol,emin,emax);
%      
% Basis of rational Krylov space
%
% WARNING: A < 0
% choose emin, emax as approximations to lambda_min(-A), lambda_max(-A)

s1=emin;
ch=0;
iter=0;
ttime=cputime;
tet = 10.; hinorm=0.;
[n,n]=size(A);
E=speye(n);EL=E;EU=E;
Lres=EL\res;
nrmb=Lres'*Lres;
V(:,1) = Lres/sqrt(nrmb);             % Include rhs in the basis.
beta2=(V(:,1)'*Lres)^2;

H=sparse(m+2,m+1);
ind=[]; err=[];
nrmrestot=[];
nrma=norm(A,'fro');
singE=condest(E)/norm(E,'fro');

if (norm(A-E-(A-E)',1)<1e-12), symm=1; else symm=0;end

 newAv=EL\(A*(EU\V(:,1)));
 K=full(V'*newAv);   % K=V'*(L\(A*(L'\V)));
 s=s1;
 eH=eig(K);
 eHpoints = sort([s1,emax]);
 snew=newpolei(eHpoints,eH,s); 
 if real(snew)<0, snew=-real(snew)+sqrt(-1)*imag(snew);end
 s=[s,snew];

cmplxflag=0;
itsinner=0;
i=0;
% additional steps
while i < m
  i=i+1;

    i1=i+1; it = 0; t=0.;
    w=EL*V(:,i);
    if (iter)  
      if symm
        [UA]=cholinc(-(A-snew*E),1e-2); wrk=pcg1(A-snew*E,w,0*w,500,tol1,UA',-UA);
      else
        [LA,UA]=luinc((A-(snew)*E),1e-2);      Prec.L=LA; Prec.U=UA;
        [wrk,resinner,itsinner] = idrs_new(A-snew*E,w,4,tol1,n,Prec);
				%wrk = bicgstab(A-snew*E,w,tol1,n,LA,UA,0*w);
      end
    else 
%disp(full([i,snew]))
      wrk = (A-snew*E)\w;   %V(:,1)=L'*((A-s1*E)\(L*V(:,1)));
    end
    wrk= EU*wrk;              %V(:,i1)=L'*((A-snew*E)\(L*V(:,i)));
    for it=1:2,
      for j=1:i
        t = V(:,j)'*wrk;
        hinorm = hinorm + abs(t^2); H(j,i)=H(j,i)+t;
        wrk=wrk-t*V(:,j);
      end;
      t = norm(wrk);
    end;
    H(i1,i)=t;
    if (t ~= 0.)
      t=1.0/t; V(:,i1)=wrk*t; 
    % if (itp==2), V(:,i:i1)=orth([real(V(:,i)),imag(V(:,i1))]);end
    end;

%rhs(i1,1)=V(:,i1)'*Lres;
 ih1=i1; ih=i;
 newAv=EL\(A*(EU\V(:,ih1)));
 D= diag(s(2:ih1)); 
 g = V(:,1:ih)'*newAv;
 K=full( ( H(1:ih,1:ih)*D + speye(ih)-g*H(ih1,ih)*[sparse(1,ih-1),1])/H(1:ih,1:ih) );

Htrick = [H(1:ih,1:ih-1),V(:,1:ih)'*(A*V(:,ih))];
Dtrick = diag([1./s(2:ih).';0]);
Ktrick = full(Htrick/(eye(ih)+Htrick*Dtrick));
%K,Ktrick
K=V(:,1:ih)'*A*V(:,1:ih);
%pause

 ind=[ind,ih];


% New poles and zeros

 eH=sort(eig(K(1:ih,1:ih))); eHorig=eH;

 if (ch)                     % Complex poles. Compute set for next complex pole of r_m

    if (sum(abs(imag(eH))>1e-4)>0 & length(eH)>2),    % Roots lambdas come from convex hull too
      ij=convhull(real(eH),imag(eH)); eH=eH(ij);
      ieH=length(eH); missing=ih-ieH;
      while missing>0,                         % include enough points from the border
        neweH=(eH(1:ieH-1)+eH(2:ieH))/2;missing=ih-length(eH);
        eH=[eH;neweH];
      end
      eH=eH(1:ih);
      eHpoints=-eH;
      eH=eHorig;
    else                                  % if all real eigs, no convex hull possible
      eHpoints = sort([s1; emax.';-real(eH)]);
    end


 else   % Real poles s from real set. Compute complex roots of r_m via Ritz convex hull
    if (sum(abs(imag(eH))>1e-2)>0 & length(eH)>2),    % Roots lambdas come from convex hull too
       eH=[eH;-s1;-emax.'];
       ij=convhull(real(eH),imag(eH)); eH=eH(ij);
      %ij=convhull(real(eH),imag(eH),{'Qt','Pp'}); eH=eH(ij);
       ieH=length(eH); missing=ih-ieH;
       while missing>0, % include enough points from the border
         neweH=(eH(1:ieH-1)+eH(2:ieH))/2;
         eH=[eH;neweH];
         missing=ih-length(eH);
       end
       eH=eH(1:ih);
     end
      eHpoints = sort([s1; emax.';-real(eH)]);
      eH=eHorig;
%     eHpoints = sort([s1; emax.';s']);
 end

 snew = newpolei(eHpoints,eH,s(2:ih1).');%pause
 if real(snew)<0, snew=-real(snew)+sqrt(-1)*imag(snew);end

% If pole is complex, include its conjugate
 if (imag(snew) ~=0), cmplxflag=1;end
 s=[s,snew];

end;

%s.'

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function r=ratfun(x,eH,s)

%c1=poly(eH); c1=c1/c1(1);
%c2=poly(s); c2=c2/c2(1);

for j=1:length(x)
r(j)=abs(prod( (x(j)-s)./(x(j)-eH) ));
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function r=ratfuni(x,eH,s)

for j=1:length(x)
r(j)=1./abs(prod( (x(j)-s)./(x(j)-eH) ));
end

return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%function snew=newpole(eHpoints,eH,s)
%
%for j=1:length(eHpoints)-1
%    sval=linspace(eHpoints(j),eHpoints(j+1),100); 
%    [sf,jx] = max (abs(ratfun(sval,eH,s)));
%    snew(j)=sval(jx);
%end
%[sn,jx]=max(abs(ratfun(snew,eH,s)));
%snew=snew(jx);
%
%return
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function snew=newpolei(eHpoints,eH,s)

eHpoints=sort(eHpoints);
for j=1:length(eHpoints)-1

    snew(j) = fminbnd( @(x) ratfuni(x,eH,s), eHpoints(j),eHpoints(j+1),optimset('TolX',1e-3));
end
[sn,jx]=max(abs(ratfun(snew,eH,s)));
snew=snew(jx);
return

