%mgm_2s.m --- Linear Multigrid method for a 2D model PDE on [0,1]^2
%         in the same usage as mgm_2d.m:         (short / compact file)
%See mgm_2d for details   (c) K Chen (2004) University of Liverpool, UK
function [u,res]=mgm(m,n,levi, gam1,gam2,nu, maxit,tol,iout,fmg,u_init)
         global iprob, if nargin<1, help mgm_2s, return; end, iprob=2; 
         global hx hy nnx nny levs  iter  U BC COF RHS R 
         if nargin<7, maxit=2; end, if nargin<8,  tol=1.0e-3; end
         if nargin<9,  iout=0; end, if nargin<10, fmg=1; end
         levs = levi; t0=cputime;   if nargin<11, u_init=zeros(m,n); end
        U=[];RHS=[];COF=[];R=[];BC=[]; hx=[];hy=[]; nnx=[];nny=[];
disp(['|___ ' mfilename ' __Linear Multigrid Method___|'...
      ' T1 to T' num2str(levs) ' ___| with ' num2str(m) ' X ' num2str(n)])
fprintf('\t\tUsing %d levels, gamma_1=%d,gamma_2=%d,%d-cycling\n',...
        levs,gam1,gam2,nu), maxi=floor(log(min(m,n))/log(10)); nus=[]; 
for i=1:levs  
   nnx(levs-i+1)=(m+1)/2^(i-1) - 1; nny(levs-i+1)=(n+1)/2^(i-1) - 1; 
end, space='                                                      ';
for i=1:levs  
 hx(i)=1/(nnx(i)+1);  hy(i)=1/(nny(i)+1); id = floor(log(i)/log(10)); 
      s=sprintf('U%d%s',i,space(1:maxi-id));     U=[U;s];
      s=sprintf('RHS%d%s',i,space(1:maxi-id)); RHS=[RHS;s];
      s=sprintf('COF%d%s',i,space(1:maxi-id)); COF=[COF;s];
      s=sprintf('R%d%s',i,space(1:maxi-id));     R=[R;s];
      s=sprintf('BC%d%s',i,space(1:maxi-id));   BC=[BC;s];
end 
for i=1:levs 
 eval(['global ' U(i,:) ' ' BC(i,:) ' ' COF(i,:) ' ' RHS(i,:) ' ' R(i,:) ])
end 
  Get_A_rhs_R_BC; iter=1; eval(['U' num2str(levs) '=u_init;']); 
  if fmg==1, alev=1; Goto_finest(1,alev,gam1,gam2,nu, iout); end
r_fine=sprintf('R%d=residual(COF%d,U%d,RHS%d); err2=norm(R%d,''fro'');',...
         levs,levs,levs,levs,levs);
  eval(r_fine); erratio=-1; res(iter)=err2;
while (err2>tol & erratio<1.00 & iter<maxit) 
  alev=levs; Goto_finest(0,alev,gam1,gam2,nu, iout); 
  err0=err2; eval(r_fine) 
  erratio=err2/max(err0,eps); 
  iter=iter+1; res(iter)=err2;
end 
fprintf('%s: Total Iterations=%2d Err=%5.1e (final rate %5.2f) CPU=%4.1f\n',...
   mfilename,iter, err2, erratio,cputime-t0), eval(['u=U' num2str(levs) ';']); 
function uc = restrict(u,lev); %==========================================[1]
        global hx hy nnx nny levs  
   range_x=2:2:nnx(lev); range_y=2:2:nny(lev); uc=[];
  for j=range_y, jc=floor(j/2); 
      for i=range_x, ic=floor(i/2);
      uc(ic,jc) = ( ...
          u(i-1,j-1)+u(i+1,j+1)+u(i-1,j+1)+u(i+1,j-1) +...
          2*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) +...
          4*u(i,j) )/16;
  end, end %i
function v = interpo(uc, Type,BC); % BC.4 pieces ==========================[2]
   [m n] = size(uc);  
    if Type==1 
     uc = [BC.bottom(:) uc BC.top(:)]; 
      left = [given_bc(0,0) BC.left given_bc(0,1)];
     right = [given_bc(1,0) BC.right given_bc(1,1)];
     uc = [left; uc; right];  
    else       
     uc = [zeros(m,1) uc zeros(m,1)];
     uc = [zeros(1,n+2); uc; zeros(1,n+2)];
    end 
for i=2:m+2, i2 = 2*(i-1); 
   for j=1:n+1,  j2 = 2*j;     
     v(i2-1,j2-1) = uc(i-1,j); v(i2-1,j2+1) = uc(i-1,j+1);
     v(i2+1,j2-1) = uc(i,j);   v(i2+1,j2+1) = uc(i,j+1);
     v(i2  ,j2) = ( uc(i,j)+uc(i,j+1) + uc(i-1,j) ...
                   +uc(i-1,j+1) )/4; 
     v(i2-1,j2  ) = (uc(i-1,j)+uc(i-1,j+1))/2; 
     v(i2-1,j2-1) = (uc(i-1,j)+uc(i  ,j  ))/2; 
     v(i2+1,j2  ) = (uc(i  ,j)+uc(i  ,j+1))/2; 
     v(i2  ,j2+1) = (uc(i-1,j+1)+uc(i,j+1))/2; 
   end %j
end,   v = v(2:end-1,2:end-1);
function u_out = smooth(u,lev, steps); %==================================[3]
      global hx hy U BC COF RHS R 
  eval([' global ' COF(lev,:) ' ' RHS(lev,:) ])
  eval([' rhs = ' RHS(lev,:) '; A =' COF(lev,:) ';' ])
      [m n]=size(u); h1=hx(lev);h2=hy(lev);
      L=tril(A);UU=triu(A,1); u=u(:); rhs=rhs(:);
  for iter=1:steps, u = - L\(UU*u) + L\rhs; end, u_out=reshape(u,m,n);
function Goto_finest(Type,lev_now, gam1,gam2,nu, iout); %==================[4]
     global hx hy nnx nny levs  U BC COF RHS R check alev iter 
 alev=lev_now; for i=1:levs
      eval([' global ' U(i,:) ' ' BC(i,:) ' ' RHS(i,:) ' ' R(i,:) ])
               end
  if levs<=1&Type==0, return; end,  a_pass=0;  check(1:levs)=0; 
while (alev < levs) | a_pass<1,     a_pass=1;  
  if check(alev)<nu
   if check(alev)==0&Type==0&alev<levs 
     s=sprintf('%s=sparse(%d,%d);',U(alev,:),nnx(alev),nny(alev));eval(s)
   end
   Goto_coarest(gam1,iout); 
  else
     if alev<levs, check(alev)=0; end
  end 
     if levs<=1, return; end 
  s=sprintf('v = interpo(%s,Type,%s); ',U(alev,:),BC(alev,:)); eval(s)
  if alev<levs
    alev = alev+1;s=sprintf('%s=v+%s;',U(alev,:),U(alev,:)); eval(s) 
  end 
  s=sprintf('%s=smooth(%s,%d,gam2);',U(alev,:),U(alev,:),alev); eval(s)
  check(alev)=check(alev)+1;  
end
function k = Goto_coarest(gam1,iout); %====================================[5]
     global hx hy nnx nny levs  U BC COF RHS R alev iter 
 for i=1:levs
 eval(['global ' U(i,:) ' ' BC(i,:) ' ' R(i,:) ' ' RHS(i,:) ' ' COF(i,:)])
 end
while alev>1     
    s=sprintf('%s = smooth(%s,%d, gam1); ',U(alev,:),U(alev,:),alev);eval(s)
    s=sprintf('%s = residual(%s,%s,%s);',R(alev,:),COF(alev,:),...
                U(alev,:),RHS(alev,:)); eval(s)
    s=sprintf('%s = restrict(%s,alev); ',RHS(alev-1,:),R(alev,:)); eval(s)
    alev = alev-1;
end 
  alev=1;  s=sprintf('[ m n] = size(%s); ', RHS(alev,:));eval(s)
  s=sprintf('ue = %s \\ %s(:);',COF(alev,:),RHS(alev,:));eval(s)
  s=sprintf('%s = reshape(ue,m,n);', U(alev,:));eval(s)
function r=residual(A,u,rhs) %=============================================[6]
 [m n] = size(u);
 r=reshape(rhs,m*n,1)-A*reshape(u,m*n,1); r=reshape(r,m,n); 
function a=coef_a(x,y) %==========================================PROB
    a = 1+x^2/100;
function b=coef_b(x,y) %==========================================PROB
    b = 1-y^2/100;
function c=coef_c(x,y) %==========================================PROB
    c = -x^2/100-y^2/100;
function d=coef_d(x,y) %==========================================PROB
   global iprob
   if iprob==1
    d = -5/3+(-8+(776809/64800+(11/135+(-19/150+...
      (1/75-y/100)*y)*y)*y)*y)*y+(6-3/400*y^2+...
      (-1871/64800+(1/675+(1/120+(1/75-y/100)*y)*y)*y+...
      (-y^2/100+21/400+(3/200-x/100)*x)*x)*x)*x;
   else
    d = (6960+(-14400+(0.1436519E7/100+(144-144*y)*y)*y)*y+...
      (-14400+(0.1443479E7/100+(-144+144*x)*x)*x)*x)*...
      exp(-60*( (x-1/2)^2+(y-1/2)^2));
   end
function u=given_bc(x,y) %=========================================PROB
   global iprob
   if iprob==1
    u = (x-1/2)^3+(y-1/3)^4;
   else
    u = exp(-60*( (x-1/2)^2+(y-1/2)^2 ));
   end
function Get_A_rhs_R_BC %===============================================PROB
      global hx hy nnx nny levs  U BC COF RHS R 
for i=1:levs
 eval(['global ' U(i,:) ' ' BC(i,:) ' ' COF(i,:) ' ' RHS(i,:) ' ' R(i,:) ])
end
for i=1:levs  
 bc=[]; for j=1:nny(i), bc(j) = given_bc(0,j*hy(i)); end
  s=sprintf('%s.left = bc; ', BC(i,:)); eval(s) 
 bc=[]; for j=1:nnx(i), bc(j) = given_bc(j*hx(i),0); end
  s=sprintf('%s.bottom = bc; ', BC(i,:)); eval(s) 
 bc=[]; for j=1:nny(i), bc(j) = given_bc(1,j*hy(i)); end
  s=sprintf('%s.right = bc; ', BC(i,:)); eval(s) 
 bc=[]; for j=1:nnx(i), bc(j) = given_bc(j*hy(i),1); end
  s=sprintf('%s.top = bc; ', BC(i,:)); eval(s) 
end 
 for lev=1:levs 
      eval([' bc = ' BC(lev,:) ';']), m=nnx(lev);n=nny(lev);
      h1=hx(lev);h2=hy(lev); coeffs=sparse(m*n,m*n); rhs=sparse(m,n);   
   for jg=1:n, y = jg*hy(lev);
      for ig=1:m, x=ig*hx(lev); row = (jg-1)*m + ig;
        a=coef_a(x,y)/h1^2; c=coef_c(x,y); 
        b=coef_b(x,y)/h2^2; d=coef_d(x,y); 
      for j=max(1,jg-1):min(jg+1,n) 
        if j==jg, is=max(1,ig-1):min(ig+1,m); else, is=ig; end 
      for i=is, col = (j-1)*m + i;
        if ig==i & jg==j
        coeffs(row,col) = -2*a -2*b +c;
        elseif ig==i 
        coeffs(row,col) = a;
        elseif jg==j
        coeffs(row,col) = b;
        end %if
      end %i
      end %j
      if ig==1, rhs(ig,jg) = d - a*bc.left(jg);
      elseif ig==m, rhs(ig,jg) = d - a*bc.right(jg);
      elseif jg==1, rhs(ig,jg) = d - b*bc.bottom(ig);
      elseif jg==n, rhs(ig,jg) = d - b*bc.top(ig);
      else,         rhs(ig,jg) = d;  % The middle guys
      end 
   end, end %jg
   s=sprintf('%s=rhs; %s=coeffs;',RHS(lev,:),COF(lev,:));eval(s)
   eval([ R(lev,:) ' = rhs;']) % True residual if u=0
   eval([ U(lev,:) ' = sparse(m,n);']) % Set u=0 for now
 end %lev ---------------Use keyboard command to examine other levels!!!
