%TEST GMRES using "cheb_fun" -- a Chebyshev polynomial preconditioner 
%
%Usage:  x = chebg(A,b,m,k,tol,maxit,x0)    
%            for solving Ax=b by GMRES(m) maxit steps to reach tol,
%             starting from x0 and using Chebyshev p_k(A) preconditioner
%
%% e.g. USE:  cheb_fun.m (plus Mfile internal GMRES)   
%    if ~exist('A'), A = delsq( numgrid('S',20) ); end
%    b = sum(A,2);        % the true solution is the vector of all ones
%    x = chebg(A,b,10, 3,1.0e-4,2,[]);
%    disp('-------------------------------------------------------------')
%    xg = gmres(A,b,10,  1.0e-4,4);  %%% No preconditioner by more steps
%%---------------------------------------------------------------------
function y=c_gmres(B,b,m,kk,tol,maxit,x0)    
if nargin<6, help chebg, return;end

global A Theta Delta k_cheby
   A=B; help chebg
   r=sum(abs(A-diag(diag(A)))')';
   r=[-r+diag(A) r+diag(A)];
   alp=double(min(r(:,1))); bet=double(max(r(:,2)));
   alp=alp(1); bet=bet(1); % de-sparse to a normal scalar
       if alp<=0 %% Trigger a more exact estimate of \lambda_min
          OPTS.disp=0;
          alp=eigs(A,1,'SA', OPTS);
       end % bet=eigs(A,1,'LR'),
       if alp<=0 %% Issue Warnings 
          disp('*** Warning: Develop a complex version or A unsuitable')
          alp=0.001; % Re-set
       end % bet=eigs(A,1,'LR'),
   Theta=(bet+alp)/2; Delta=(bet-alp)/2; k_cheby=kk;
   fprintf('\n%s: Eigs = [%5.3f, %5.2f], \\theta=%5.2f, \\delta=%5.2f\n',...
            mfilename,alp,bet,Theta,Delta)

y = my_gmres('cheb_fun',b,m,tol,maxit,x0);
return  %% 
%===================================== END ==================

return %=====================================================
%% gmres_k.m --- adapted
function x=my_gmres(myfun,b,k,tol,maxit,x0)
   m=length(b); n=m;  % matrix-free now
   global k_cheby
if length(x0)==0,x0=0*b; end
fprintf(' INPUT to my_gmres:  restart = %d,  maxit = %d,  tol=%e\n',...
         k,maxit,tol)
x=x0;
           r = b - feval(myfun,x,0);  norm_0=norm(r); iter=0; ok=0;

while iter<maxit

% Step (1)
          x0 = x;  r = b - feval(myfun,x,0);  r2=norm(r); Q=[];R=[];c=[];s=[];
% Step (2)
          q1 = r/r2; Q=q1;   rhs=zeros(m,1);rhs(1)=r2;
   for i=1:k
% Step (3)
          w = feval(myfun,Q(:,i),1);
% Step (4)
   for L=1:i
     R(L,i) = w'*Q(:,L);  w=w-R(L,i)*Q(:,L);
   end 
% Step (5)
   R(i+1,i) = norm(w);   Q(:,i+1)=w/R(i+1,i);
% Step (6) 
   for L=1:(i-1)
     t=c(L)*R(L,i) + s(L)*R(L+1,i);
     R(L+1,i)=s(L)*R(L,i) - c(L)*R(L+1,i);
     R(L,i)=t;
   end
% Step (7)
   if abs(R(i+1,i))<1.0e-16, c(i)=1; s(i)=0;
   elseif abs(R(i+1,i))>abs(R(i,i))
      t=R(i,i)/R(i+1,i); s(i)=1/sqrt(1+t^2);  c(i)=s(i)*t; else
      t=R(i+1,i)/R(i,i); c(i)=1/sqrt(1+t^2);  s(i)=c(i)*t; 
   end
% Step (8)
   t=c(i)*rhs(i);  rhs(i+1)=s(i)*rhs(i);  rhs(i)=t;
   R(i,i)=c(i)*R(i,i) + s(i)*R(i+1,i);
   R(i+1,i) = 0; 
% Step (9)
   y = R(1:i,1:i) \ rhs(1:i,1);
% Step (10)
    y=Q(:,1:i)*y; y = feval(myfun,y,2);  % For right GMRES
    x = x0 + y;
% Step (11)
   r=b-feval(myfun,x,0);   r2=norm(r)/norm_0;
%  fprintf('%s iter=%3d, res=%e\n',mfilename,iter*k+i,r2)
   if r2 <= tol, ok=1; break; disp('gmres_k done...'); end
% Step (12)
end % i
   if ok==1, break; end
   iter=iter+1;
end %while
   fprintf('%s: Total GMRES(%d) Steps = %d to achieve %e [N=%d]\n',...
            mfilename, k, iter*k+i, norm(r)/norm_0, m);
   fprintf('%s: Order k_Chebyshev = %d used\n\n',mfilename,k_cheby)
