function [OUT,nrm] = trfomgrid(A, lu, rd, cols, rows, degree, m, maxiter, tol, b, flag)
%TRFOMGRID: Pseudospectra with GRID, Transfer Functions and resrarted FOM
%
%        INPUT
%
%A       :input matrix
%lu      :Left up corner of mesh (i.e. 1+1*i)    
%rd      :Right bottom corner of mesh
%cols    :Number of columns of mesh
%rows    :Number of rows of mesh 
%degree  :Dimension of transfer function Krylov subspace 
%m       :Dimension of FOM Krylov space
%maxiter :Maximum number of FOM restarts
%tol     :Convergence tolerance for FOM
%b       :Starting vector for Transfer Functions. ||b||=1 (i.e. b=rand(size(A,1),1);b=b/norm(b) ) 
%flag    :if flag==1, then compute only norms of residuals of FOM
%
%        OUTPUT
%OUT     :Estimation of \sigma_{\min}(A-zI) for all mesh nodes z
%nrm     :Norms of residuals for FOM
%
%		 Constantine Bekas. 

%Start timming
tic;

[XX,YY] = meshgrid(linspace(real(lu),real(rd),cols), linspace(imag(lu),imag(rd),rows));
points_z = reshape(XX+sqrt(-1)*YY,(cols)*(rows),1);

%Perform Initial Arnoldi
%-----------------------
N = max(size(A));
b=b/norm(b);
[H,V] = reorth_arnoldi(A,b,degree);

%Initialize solution
%-------------------
fz = zeros(degree+1,length(points_z));
par_res = zeros(m,length(points_z));

%Find the current seed remainder
%-------------------------------
r = V(:,degree+1);

%tm = 0;
%tic;

%Do arnoldi for the seed system
[H1,V1]=reorth_arnoldi(A,r,m);

%tm = tm + toc;

%Outer restated iteration
%------------------------
active_points_z = 1:length(points_z);
nrm = ones(length(points_z),1);

for k = 1 : maxiter    
    
    if (k > 1)
        r = V1(:,m+1);        
        %Do arnoldi for the seed system
        [H1,V1]=reorth_arnoldi(A,r,m);                               
    end
    
    
    %Inner point iteration
    for l = 1 : length(active_points_z)
              
        par_res(:,active_points_z(l))=(H1(1:m,1:m)-points_z(active_points_z(l))*eye(m))\(nrm(active_points_z(l))*[1;zeros(m-1,1)]);

        nrm(active_points_z(l)) = - H1(m+1,m)*par_res(m,active_points_z(l));        
        
    end
    
    fz(:,active_points_z) = fz(:,active_points_z) + (V(1:N,1:degree+1)'*V1(:,1:m))*par_res(1:m,active_points_z);    

    sprintf('Number of remaining mesh points: %d \n',length(active_points_z))
    
    active_points_z = find(abs(nrm) > tol);
    if (isempty(active_points_z))
        break
    end
end

if flag
    hold on
    plot(eig(H1(1:m,1:m)),'r.');
    plot(eig(H(1:degree,1:degree)),'b.');
    return
end

tm1 = toc;

tic;
 for k=1:length(points_z)
     %Compute Transfer
     RES = eye(degree+1,degree)- H(degree+1,degree)*[zeros(degree+1,degree-1), fz(1:degree+1,k)]; 
     RES = [((H(1:degree,1:degree)-points_z(k)*eye(degree))'\RES')',fz(:,k)];
     OUT(k)=norm(RES,2);
 end
 
 OUT = (1./(reshape(OUT,rows,cols)));
 nrm = abs(reshape(nrm,rows,cols));
 tm2 = toc;
disp(' ');
s= sprintf('Restarted FOM ellapsed time (secs) %f:\n ', tm1);
disp(s);

tm2 = toc;
s=sprintf('Time on mesh (secs) %f: \n', tm2 );
disp(s)
tm3 = tm1+tm2;
s=sprintf('Total ellapsed time (secs) %f: \',tm3);
disp(s);