function up=pcsmulsw(au)
% PCSMULSW
% symmetric multiplicative Schwarz as a precondtioner on an m x m mesh
%
% function up=pcsmulsw(au)
%
global ph pv overlap
[m1,m2]=size(au); m=sqrt(m1);
au2=zeros(m,m); 
au2(:)=au;
%ph=2; pv=2; overlap=max(3,floor(m/5));
up2=mutswp(m,au2,overlap,ph,pv);
up=up2(:);
%
% This is a modification of the solver. Divide the domain into ph * pv
% subdomains, ph splits in horozontal and pv splits in vertical
%
function u2d=mutswp(m,rhs2d,overlap,ph,pv)
[l,r]=partition(m,2,overlap);
[lv,rv]=partition(m,pv,overlap);
[lh,rh]=partition(m,ph,overlap);
h=1/(m+1);
u2d=zeros(m,m);
maxit=1;
for i=1:maxit
	%
	% forward sweep
	%
	for ih=1:ph
		for iv=1:pv
			%
			% residual update
			%
			r2d=res2d(u2d,rhs2d);
			lbh=lh(ih); ubh=rh(ih); lbv=lv(iv); ubv=rv(iv);
			vd=ubv-lbv+1; hd=ubh-lbh+1;
			rsub=r2d(lbh:ubh,lbv:ubv);
			ux=subdomain(rsub,hd,vd,h);
			u2d(lbh:ubh,lbv:ubv)=u2d(lbh:ubh,lbv:ubv)+ux;
		end
	end
	%
	% backward sweep
	%
	for ih=ph:-1:1
		for iv=pv:-1:1
			%
			% residual update
			%
			r2d=res2d(u2d,rhs2d);
			lbh=lh(ih); ubh=rh(ih); lbv=lv(iv); ubv=rv(iv);
			vd=ubv-lbv+1; hd=ubh-lbh+1;
			rsub=r2d(lbh:ubh,lbv:ubv);
			ux=subdomain(rsub,hd,vd,h);
			u2d(lbh:ubh,lbv:ubv)=u2d(lbh:ubh,lbv:ubv)+ux;
		end
	end
end
function r2=res2d(uin,rhs2d)
[m1,m2]=size(rhs2d);
if m1 ~= m2
	disp('error in res2d');
	return;
end
r2=zeros(m1,m1);
uin1d=uin(:); rhs1d=rhs2d(:);
res1d=rhs1d-lapmf(uin1d);
r2(:)=res1d;