function sol=fmg(f,nu0,nu1,nu2)
% FMG recursive fmg as on page 43, builds and uses the real data structure
%
nt=length(f); l=log2(nt-1);
lt=2^(l+1)-2+l;
ifu=lt; ifl=lt-2^l;
ft=zeros(lt,1); vt=zeros(lt,1); ft(ifl:ifu)=f;
[vt,ft]=fmgx(vt,ft,nu0,nu1,nu2);
sol=vt(ifl:ifu);
%
% internal fmg code
%
function [vout,fout]=fmgx(vin,fin,nu0,nu1,nu2)
%
% cheap hack to back out the level
%
ntl=length(fin); l=max(1,floor(log2(ntl))-1);
%
% lower and upper bounds for the current mesh and the one below
%
vout=vin; fout=fin;
ifu=ntl; ifl=ntl-2^l; icu=ifl-1; icl=icu-2^(l-1);
fx=zeros(ifl,1);
nl=2^l; nt=nl+1;
lmin=3;
if l==lmin
	vout(ifl:ifu)=zeros(nt,1);
else
	fout(icl:icu)=ftoc(fin(ifl:ifu));
        [vout(1:icu),fx(1:icu)]=fmgx(vout(1:icu),fout(1:icu),nu0,nu1,nu2);
	vout(ifl:ifu)=ctof(vout(icl:icu));
end
for i=1:nu0
         [vout(1:ifu),fout(1:ifu)]=vcycle(vout(1:ifu),fin(1:ifu),nu1,nu2);
end