function [Clevels,condCA,DAmin,DAmax] = mg1dprec(fce_tau,fce_q,fce_f,a,fa,xl,el,levels,smoothingsteps)

% function [Clevels,condCA,Amin,Amax,Dmin,Dmax] = mg1dprec(fce_tau,fce_q,fce_f,a,fa,xl,el,levels,smoothingsteps)

l = 0;
[Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
Cl = inv(Al);

for l=1:levels
    [xl,el,El] = refine(xl,el);
    
    Cl = El*Cl*El';
    
    [Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
    
    Dl = diag(diag(Al)); % Jacobi smoother
    
    Il = eye(max(size(xl)));
    
    for s=1:smoothingsteps
        Cl = Cl + inv(Dl)*(Il-Al*Cl);
    end
end

Clevels = Cl;

% for test purposes only
e = abs(eig(Clevels*Al));
condCA = max(e)/min(e);

DA = inv(diag(diag(Al)))*Al;
e = abs(eig(DA));
DAmin = min(e);
DAmax = max(e);
