% SCHWARZ_TEST
%
% Run a test of Schwarz methods as preconditioners. 
%
global ph pv overlap ac zm lc uc lf uf at global_levels
%
% Set up the problem by table number.
%
table_number=5.3;
%
%
maxits=200; splits=4; tol=1.d-2; ngrids=4;
%
% Tables 4.1-4.4 are the one level methods; 5.1-5.3 the two level.
%
if table_number < 5
    levels=1;
    if table_number < 4.3
    ngrids=3; % Using the finest grid with some one-level methods 
              % takes too long. You'd be bored.
    end
else
    levels=2;
end
global_levels=levels;
%
% Now it's time to figure out the overlap. Your choices are 0, 1, 
% 1% , and 2%
%
% I've picked 2%. You could loop over all the possibilities, but
% that would eat up some time.
%
for p=1:ngrids
    mp=p+4;
%   olap(p) = 0;              % none
%   olap(p) = 1;              % 1
%   olap(p) = ceil(2^mp*.01); % one percent
    olap(p) = ceil(2^mp*.02); % two percent
end
%
%
%
itcount=zeros(splits,ngrids);
splim=splits*ones(10,1);
for p=1:ngrids
	mp=p+4;
        overlap=olap(p);
	m=2^mp+overlap; h=1/(m+1); x=1:m; x=x'*h;
        if levels == 2
              at=pmat2(m,m,h);
        end
	%
	% Fix the right side so that the solution is
	%
	% u(x,y) = sin(pi x) sin(pi y) e^x
	%
	wx=sin(pi*x); wxp=pi*cos(pi*x); wxpp=-pi*pi*wx;
	wy=wx; wypp=wxpp;
	u=(wx.*exp(x))*wy';
	wxx=(wx + 2*wxp + wxpp).*exp(x); uxx=wxx*wy';
	wyy=wypp; uyy=(wx.*exp(x))*wyy';
	rhs2d=-(uxx+uyy); records=zeros(maxits+1,splits); rlen=0;
%
%       Now loop over the decompositions
%
	for id=1:splits
		ph=2^(id+1); pv=ph;	
        if levels == 2
%
% form and factor the coarse mesh matrix
%
		[ac,zm]=acform(m,ph,pv,overlap); uc=chol(ac); lc=uc';
%
       end
%
% Factor the matrix af for this mesh size
%
%
 		af=afform(m,ph,pv,overlap); uf=chol(af); lf=uf';
%
%               Use an LU factorization if the problem is non-symmetric
%
%		af=afform(m,ph,pv,overlap); [lf,uf]=lu(af);
		rhs1d=rhs2d(:); u0=zeros(m*m,1); 
%
% Here are the methods.
%
switch table_number
%
  case 4.1
% One-level symmetric multiplicative with CG
	[ug,errors,tits]=pcgsol(u0,rhs1d,'lapmf',[tol,maxits],'pcsmulsw');
%
  case 4.2
% One-level multiplicative with GMRES; Table 4.2
	[ug,errors,tits]=gmres2(u0,rhs1d,'mullap',[tol,maxits]);
%
  case 4.3
% One level additive with CG; Table 4.3
	[ug,errors,tits]=pcgsol(u0,rhs1d,'lapmf',[tol,maxits],'pcaddsw');
%
  case 4.4
% One level additive with GMRES; Table 4.4
        [ug,errors,tits]=gmres2(u0,rhs1d,'addlap',[tol,maxits]);
%
  case 5.1
% Two-level additive with CG; Table 5.1
	[ug,errors,tits]=pcgsol(u0,rhs1d,'lapmf',[tol,maxits],'pcaddsw');
%
  case 5.2
% Two-level additive with GMRES; Table 5.2
        [ug,errors,tits]=gmres2(u0,rhs1d,'addlap',[tol,maxits]);
%
  case 5.3
% Two-level hybrid with GMRES; Table 5.3
	[ug,errors,tits]=gmres2(u0,rhs1d,'hyblap',[tol,maxits]);
%
end % end of switch
%
		errors=errors/errors(1);
		records(1:tits+1,id)=errors';
		rlen=max(rlen,tits+1);
		delta=norm(lapmf(ug)-rhs1d,inf);
		itcount(id,p)=tits;
	end
	records(1:rlen,:)
end
%
% itcount is the array for Tables 4.1, 4.2, 4.3, 4.4, 5.1, 5.2, 5.3
%
itcount