function mg_nr(m,nu1,nu2,mu,u,ur)
%
%   Multigrid cycle
%

u = zeros(0,1);
ur = zeros(0,1);
nk = zeros(0,0);

% Start on the finest level -- m.

k = m;

while k <= m
    
    % Relax nu1 times (pre-relaxations).
    
    for i = 1:nu1
        relax(k)
    end
    
    % If k=1, start going up to next-finer grid.
    
    while k > 1
        
        % Restrict residuals and form right-hand sides on next coarser grid.
        
        restrict(k,k-1)
        
        k = k - 1;
        
        % Set nk(k) initially equal to the cycle index ivw
        % (1 for V cycle, 2 for W cycle, etc.)
        
        nk(k) = mu;
        
        % Continue down to next-coarser grid.
        
        if k == 1
            for i = 1:nu1
                relax(k)
            end
            break
        end
        
        break
    end
    
    if k > 1
        continue
    end
    
    
    while k <= m
        
        % nr2 is used simply to allow extra sweeps on coarsest level. nu2 is
        % the number of post-relaxation sweeps.
        
        nr2 = nu2;
        
        if k == 1
            nr2 = nu2;
        end
        
        % Relax nr2 times.
        
        for i = 1:nr2
            relax(k)
        end
        
        % If you're back on the finest level, then you have completed the
        % cycle.
        
        if k == m
            return
        end
        
        % nk(k) counts how many cycles you have performed on intermediate
        % level k.
        
        nk(k) = nk(k) - 1;
        
        % If nk(k)=0, then you have completed ivw cycles on this level. 
        % Otherwise, you need to go down to the next-coarser level again.
        
        if nk(k) > 0
            break
        end
        
        % Interpolate and add the correction to the solution of next-finer
        % grid.
        
        intadd(k+1,k)
        
        % Continue up
        
        k = k + 1;
    end
    
end
