function solu = GeoMultiGrid(gmgsystem,solu,Rhs,level,smooth,GMGIter)
%
% Solve the linear system with the geometry multigrid method
% smooth and correction 
% 
for i = 1 : GMGIter
    solu = GMGIteration(gmgsystem,solu,Rhs,level,smooth);
end
%--------------sub functions-----------------------------------------------
function solu = GMGIteration(gmgsystem,solu,Rhs,level,smooth)

%--- First, do the smooth iteration 
dim = size(solu,2);
[A,Prolong] = GetMatrix(gmgsystem,level);
if level == 1
    solu = A\Rhs;
else
    for i = 1:dim
        %level
        %disp('Pre smooth')
        [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 = GMGIteration(gmgsystem,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