function uplus = ras(u, coeffs, nsd, msd, lap, f)
%
% Applies the restricted additive Schwarz method to an input u for a
% convection-diffusion system 
%
%        Lu = Laplacian(u) + c*u + d*u_x + e*u_y = f 
%
% The discretization is by centered differences on an (m x m)
% square grid in [0,1]x[0,1].
% 
% Restricted additive Schwarz is applied on an (nsd x nsd) array 
% of subdomains with an (msd x msd) square grid in each subdomain 
% (so m = nsd*msd) and lap grid lines of overlap. 
%
% Inputs (required): 
%   u       = (nsd*msd)*(nsd*msd) input array.  
%   coeffs  = row vector of coefficients [c d e] in the operator L.  
%   nsd     = number of subdomains per axis => (nsd x nsd) array 
%             of subdomains.  
%   msd     = number of subdomain gridpoints per axis 
%             => (msd x msd) grid per subdomain.  
%   lap     = number of gridlines of overlap in additive Schwarz.  
%   f       = (nsd*msd)*(nsd*msd) right-hand side array. 
%
% Output: 
%   uplus   = (nsd*msd)x(nsd*msd) output array.

if nargin~=6, 
  error('ras: Number of arguments must be exactly 6.'); 
end
if norm(double(size(u) ~= size(f))), 
  error('ras: u and f must be the same size.'); 
end
if nsd < 1, error('ras: nsd must be >= 1.'); end 
if msd < 1, error('ras: msd must be >= 1.'); end
if lap < 0, error('ras: lap must be >= 0.'); end
if nsd == 1 & lap > 1, 
  error('ras: lap must be 0 or 1 when nsd = 1.'); 
end
if lap < 1 | lap > msd+1, error('ras: Must have 1 <= lap <= msd+1.'); end 

global h_GL;

m = nsd*msd; 
h_GL = (m+1)\1; 

d = coeffs(1,2); 
e = coeffs(1,3); 

% Set up the necessary matrices. 

if nsd == 1
  L = condif_DD(coeffs,m,m);
else
  Lcorner = condif_DD(coeffs,msd+lap-1,msd+lap-1);
  if nsd >= 3
    Lside = condif_DD(coeffs,msd+2*lap-2,msd+lap-1);
    Ltopbot = condif_DD(coeffs,msd+lap-1,msd+2*lap-2); 
    Lint  = condif_DD(coeffs,msd+2*lap-2,msd+2*lap-2);
  end
end

if nsd == 1
  % Treat nsd = 1 as a special case. 
  rhs = reshape(f,m*m,1);
  uplus = reshape(L\rhs,m,m);
%  return;
else
  % Embed u in a larger array that includes the boundary
  % values. For the convection-diffusion problem with Dirichlet BCs,
  % these are all zero.
  mp2 = m+2;
  uu = zeros(mp2,mp2);
  uu(2:m+1,2:m+1) = u; 

  % Cycle through the subdomain solves. 

  for p = 1:nsd
    if p == 1 
      mb = 1; mt = msd+lap+1; 
      mbf = 1; mtf = msd+lap-1;
      for q = 1:nsd
        if q == 1, ml = 1; mlf = 1; 
        else ml = (q-1)*msd-lap+2; mlf = (q-1)*msd-lap+2; end
        if q == nsd, mr = mp2; mrf = m; 
        else mr = q*msd+lap+1; mrf = q*msd+lap-1; end
        rhssd_vec = reshape(f(mbf:mtf,mlf:mrf)',(mtf-mbf+1)*(mrf-mlf+1),1); 
        [bc_Lap, bc_dx, bc_dy] = bc_vectors_DD(uu(mb:mt,ml:mr)); 
        if q == 1 
          uuxsd = Lcorner\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(1:msd,1:msd) = uuxsd(1:msd,1:msd);
        elseif 1 < q & q < nsd
          uuxsd = Ltopbot\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(1:msd,(q-1)*msd+1:q*msd) = uuxsd(1:msd,lap:lap+msd-1);
        else
          uuxsd = Lcorner\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(1:msd,m-msd+1:m) = uuxsd(1:msd,lap:lap+msd-1);
        end
      end
    elseif 1 < p & p < nsd
      mb = (p-1)*msd-lap+2; mt = p*msd+lap+1; 
      mbf = (p-1)*msd-lap+2; mtf = p*msd+lap-1; 
      for q = 1:nsd
        if q == 1, ml = 1; mlf = 1; 
        else ml = (q-1)*msd-lap+2; mlf = (q-1)*msd-lap+2; end
        if q == nsd, mr=mp2; mrf = m; 
        else mr = q*msd+lap+1; mrf = q*msd+lap-1; end
        rhssd_vec = reshape(f(mbf:mtf,mlf:mrf)',(mtf-mbf+1)*(mrf-mlf+1),1);
        [bc_Lap bc_dx, bc_dy] = bc_vectors_DD(uu(mb:mt,ml:mr)); 
        if q == 1 
          uuxsd = Lside\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus((p-1)*msd+1:p*msd,1:msd) = uuxsd(lap:lap+msd-1,1:msd);
        elseif 1 < q & q < nsd
          uuxsd = Lint\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus((p-1)*msd+1:p*msd,(q-1)*msd+1:q*msd) = uuxsd(lap:lap+msd-1,lap:lap+msd-1);
        else
          uuxsd = Lside\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus((p-1)*msd+1:p*msd,m-msd+1:m) = uuxsd(lap:lap+msd-1,lap:lap+msd-1);
        end
      end
    else
      mb = m-msd-lap+2; mt = mp2; mbf = m-msd-lap+2; mtf = m; 
      for q = 1:nsd
        if q == 1, ml = 1; mlf = 1; 
        else ml = (q-1)*msd-lap+2; mlf = (q-1)*msd-lap+2; end
        if q == nsd, mr = mp2; mrf = m; 
        else mr = q*msd+lap+1; mrf = q*msd+lap-1; end
        rhssd_vec = reshape(f(mbf:mtf,mlf:mrf)',(mtf-mbf+1)*(mrf-mlf+1),1);
        [bc_Lap, bc_dx, bc_dy] = bc_vectors_DD(uu(mb:mt,ml:mr)); 
        if q == 1 
          uuxsd = Lcorner\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(m-msd+1:m,1:msd) = uuxsd(lap:lap+msd-1,1:msd);
        elseif 1 < q & q < nsd
          uuxsd = Ltopbot\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(m-msd+1:m,(q-1)*msd+1:q*msd) = uuxsd(lap:lap+msd-1,lap:lap+msd-1);
        else
          uuxsd = Lcorner\(rhssd_vec - (h_GL*h_GL)\bc_Lap -((2*h_GL)\d)*bc_dx-((2*h_GL)\e)*bc_dy);
          uuxsd = reshape(uuxsd,mr-ml-1,mt-mb-1)';
          uplus(m-msd+1:m,m-msd+1:m) = uuxsd(lap:lap+msd-1,lap:lap+msd-1);
        end
      end
    end
  end
end
