function [U]=eksm(A,rhs,m)
%
% Generate basis for Extended Krylov space for   (A, rhs)
%
%

n=length(rhs);
maxm=m;
LM=speye(n);UM=LM;
     [LA,UA]=lu(A);
%    [LM,UM]=lu(M);

s=2;
rhs1=UM*(UA\(LA\(LM*rhs)));
%rhs1=M*(UA\(LA\(rhs)));
U=[rhs,rhs1];

rho=0;
H=sparse((m+1)*s,m*s);T=H;
l=sparse((m+1)*s,m);

%[ibeta,U]=gram_sh(U);
[U,bb]=qr(U,0); ibeta=inv(b);
j=0;

   while j<m 
    j=j+1;

    jms=(j-1)*s+1;j1s=(j+1)*s;js=j*s;js1=js+1;
    Up(:,1) = LM\(A*(UM\U(:,jms))); %Up(:,1) = A*(UM\(LM\U(:,jms))); 
    Up(:,2) = UM*(UA\(LA\(LM*U(:,js)))); %Up(:,2) = M*(UA\(LA\U(:,js)));

%new bases block (modified gram)
    for il=1:2
        k_min=1;  %max(1,j-k_max);
        for kk=k_min:j
            k1=(kk-1)*s+1; k2=kk*s;
            coef= U(1:n,k1:k2)'*Up;
            H(k1:k2,jms:js) = H(k1:k2,jms:js)+ coef; 
            Up = Up - U(:,k1:k2)*coef; 
        end
    end

    if (j<=m)
    %  [hinv,U(1:n,js1:j1s)]=gram_sh(Up);
    %  H(js1:j1s,jms:js)=inv(hinv);
       [U(:,js1:j1s),hinv]=qr(Up,0);
       H(js1:j1s,jms:js)=hinv;
       hinv = inv(hinv);
    end
  
    I = speye(j*s+s); 

    if (j==1),
      l(1:j*s+1,j) =...
       [ H(1:s+1,1), speye(s+1,1)]*(ibeta(1:2,2)/ibeta(1,1));
    else
      l(1:j*s+s,j) = l(1:j*s+s,j) + H(1:j*s+s,js-1)*rho;
    end

    T(1:js+s,js-1)=H(1:js+s,js-1); T(1:js+s,js)=l(1:js+s,j);
    l(1:j*s+2,j+1) = ( I(1:j*s+2,j*s)- T(1:js+s,1:js)*H(1:js,j*s))*hinv(2,2);

    rho = hinv(1,2)/hinv(1,1);

 end
      
