function [U]=block_Krylov(A,B,m)
%function [U]=block_Krylov(A,B,m)
%
%block arnoldi process, with modified gram orthogon to build basis for block standard Krylov space
%

[n,s]=size(B); 
[U(:,1:s),beta]=qr(B,0);
H=zeros((m+1)*s,m*s);
for j=1:m,

    jms=(j-1)*s+1;j1s=(j+1)*s;js=j*s;js1=js+1;
    Up = A*U(:,jms:js);

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

    if (j<=m)
     [U(1:n,js1:j1s),h]=qr(Up,0);
     H(js1:j1s,jms:js)=h;
    end

end;

