function [V,info] = DMG3(U,F,par)

% function [V,info] = DMG3(U,F,par)
% dynamic multigrid iteration 
% start on coarsest grid
% 
% data
% U{1:max_level}(:,:,:): solution
% F{1:max_level}(:,:,:): right hand side
% 
% parameters with defaults
% par.max_level [5]: number of grids, h = 1/2^level
% par.tol_high [0.01]*h^2: fine grid tolerance 
% par.tol_low [0.2]: coarse grid tolerance 
% par.rate [0.8]: smoothing rate 
% par.relax [1.5]: relaxation factor for smoothing
% par.max_time [100]: maximal computing time
%
% output
% V: fine grid solution
% info.levels: levels
% info.times: computing times

% initialization
r_F = ones(par.max_level,1)*norm(F{par.max_level}(:),inf);
r_old = inf;
tol = par.tol_high*(1/4).^[1:par.max_level];
level = 1; info.levels = [1];
t_start = clock; info.times = etime(clock,t_start);

% dynamic iteration
while (level <= par.max_level) & (info.times(end) <= par.max_time);
   [U{level},r_new] = Smooth(U{level},F{level},par.relax);
   q = r_new/r_old; r_old = r_new;
   % convergence test
   if r_new < tol(level)*r_F(level);
      level = level + 1;
      if level <= par.max_level
         tol(level-1) = par.tol_low;
	 U{level} = U{level} + Extend(U{level-1});
	 r_old = inf;
      end;
   % smoothness test (smooth = slow reduction)
   elseif (q > par.rate) & (level > 1);
      F{level-1} = Restrict(F{level}-System(U{level}));
      level = level-1;
      U{level} = 0*F{level};
      r_F(level) = norm(F{level}(:),inf);
      r_old = inf;
   end;
   % statistics
   info.levels = [info.levels; level];
   info.times = [info.times; etime(clock,t_start)];
end;

V = U{par.max_level};


%%%%% local programs 

function [U,r] = Smooth(U,F,s)

% function [U,r] = Smooth(U,F,s)
% damped Jacobi step
% x_j,k,l=x_j,k,l + (s/12 h^2) (f-(6x_j,k-...)/h^2) 
% r: norm of residual

n = size(U,1); 
R = F - System(U);
r = norm(R(:),inf);
U = U + (s*(n-1)^(-2)/12)*R;


function Y = System(X)

% function Y = System(X)
% multiplication with discrete Poisson matrix 

n = size(X,1); 
I = [2:n-1]; 
Y = 0*X;
Y(I,I,I) = ((n-1)^2) * (6*X(I,I,I) - ... 
   X(I-1,I,I) - X(I+1,I,I) - ... 
   X(I,I-1,I) - X(I,I+1,I) - ...
   X(I,I,I-1) - X(I,I,I+1));


function V = Restrict(U)

% function V = Restrict(U)
% transfer to coarse grid

n = size(U,1);
V = U(1:2:n,1:2:n,1:2:n);


function U = Extend(V)

% function U = Restrict(V)
% transfer to fine grid

n = 2*size(V,1)-1; 
U = zeros(n,n);
I = [2:2:n]; J = [1:2:n];
U(J,J,J) = V;
U(I,J,J) = (U(I-1,J,J) + U(I+1,J,J))/2;
U(:,I,J) = (U(:,I-1,J) + U(:,I+1,J))/2;
U(:,:,I) = (U(:,:,I-1) + U(:,:,I+1))/2;