function solu = MGIteration(mg,solu,Rhs,level,smooth,MGIter)
%
% Solve the linear system with the multigrid iteration:
% INPUT: mg: the multigrid data
%        solu: the initial solution
%        Rhs: the right hand side term
%        level: the current level of the sequence of space
%        smooth: the number of CG iteration steps for the pre- and post-smoothing
%        steps
%        MGIter: the number of multigrid iteration steps
% OUTPUT: solu: the resultant solution
%
% Hehu Xie 2015-10-09

for i = 1 : MGIter    
    solu = MGITStep(mg,solu,Rhs,level,smooth);
end
%--------------sub functions-----------------------------------------------
function solu = MGITStep(mg,solu,Rhs,level,smooth)

%--- First, do the smooth iteration -------------
dim = size(solu,2);
[A,Prolong, tmp, bd]  =  MGGetMatrix( mg, level );

if level == 1
    if length(bd) > 0
        Rhs(bd,:) = 0;
    end
    solu = A\Rhs;
    
else
    for i = 1:dim
        %level
        %disp('Pre smooth')
        if length(bd) > 0
            Rhs(bd, :) =0;
        end
        [solu(:,i), error, iter, flag] = cg(A, solu(:,i), Rhs(:,i), smooth, eps);
        %error
        %iter
        %flag
    end
    % do the correction 
    residual = Prolong'*(Rhs - A*solu);
    dim1 = size(Prolong,2);
    solu0 = sparse(dim1,dim);
    solu0 = MGITStep(mg,solu0,residual,level-1,smooth);
    solu = solu + Prolong*solu0;
    % do the post smoothing
    for i = 1:dim
        %disp('post smoothing')
        [solu(:,i), error, iter, flag] = cg(A, solu(:,i), Rhs(:,i), smooth,eps);
        %error
        %iter
        %flag
    end
end