force[{x_,y_}]:=-Exp[-(x-0.25)^2-(y-0.6)^2];

findRes[u_,f_]:=Module[{n,h,i,j,r},
{n,n}=Dimensions[u];
Print["in findRes,n=",n];
Print["in findRes,u=",u];
h=1/(n-1);
r=Table[f[[i,j]]-(1/h^2)*(u[[i,j-1]]+u[[i,j+1]]+u[[i+1,j]]+u[[i-1,j]]-4*u[[i,j]]),{i,2,n-1},{j,2,n-1}];
r
];

relax[u_,f_]:=Module[{i,j,A,n,h},
Print["in relax"];
{n,n}=Dimensions[u];
Print["in relax,n=",n];
h=1/(n-1);
A={{0,1,0},{1,0,1},{0,1,0}};
u=ReplacePart[u,{i_,j_}/;(i>1&&i<n&&j>1&&j<n&&EvenQ[i+j]):>(1/4)Total[A*u[[i-1;;i+1,j-1;;j+1]],2]-((h^2)/4)f[[i,j]]  ];
u=ReplacePart[u,{i_,j_}/;(i>1&&i<n&&j>1&&j<n&&OddQ[i+j]):>(1/4)Total[A*u[[i-1;;i+1,j-1;;j+1]] ],2]-((h^2)/4)f[[i,j]] ;
u
];
Attributes[relax]={HoldFirst};

f2c[f_]:=Module[{mask,i,j,c,newn},
Print["in f2c"];
{n,n}=Dimensions[c];
If[n<3,Throw["n<=3 found"]];
mask=(1/16)*{{1,2,1},{2,4,2},{1,2,1}};

For[i=3,i≤n-2,i=i+2,
For[j=3,j≤n-2,j=j+2,
f[[i,j]]=Total[f[[i-1;;i+1,j-1;;j+1]].mask,2];
]
];
newn=Ceiling[n/2];
c=Table[0,{i,newn},{j,newn}];
c[[2;;-2,2;;-2]]=f[[Range[3,n-2,2],Range[3,n-2,2]]];
c
];

Attributes[f2c]={HoldFirst};

c2f[c_]:=Module[{f,n,i,j},
Print["in c2f"];
{n,n}=Dimensions[c];
If[n<3,Throw["n<3 found"]];
n=2*(n-1)+1;
f=Table[0,{i,n},{j,n}];
i=Range[3,n-2,2];
j=Range[3,n-2,2];
f[[i-1,j-1]]=f[[i-1,j-1]]+(1/4)*c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i-1,j]]=f[[i-1,j]]+(1/2)*c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i-1,j+1]]=f[[i-1,j+1]]+(1/4)*c[[Ceiling[i/2],Ceiling[i/2]]];

f[[i,j-1]]=f[[i,j-1]]+(1/2)*c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i,j]]=f[[i,j]]+c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i,j+1]]=f[[i,j+1]]+(1/2)*c[[Ceiling[i/2],Ceiling[i/2]]];

f[[i+1,j-1]]=f[[i+1,j-1]]+(1/4)*c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i+1,j]]=f[[i+1,j]]+(1/2)*c[[Ceiling[i/2],Ceiling[i/2]]];
f[[i+1,j+1]]=f[[i+1,j+1]]+(1/4)*c[[Ceiling[i/2],Ceiling[i/2]]];

];
Attributes[c2f]={HoldFirst};


vcycle[u_,f_,mu1_,mu2_]:=Module[{i,j,r,correct,n,h },
{n,n}=Dimensions[u];
Print["in vcyle,n=",n]; 
h=1/(n-1);
If[n<3,
{
Print["in vcyle,n==3 case"]; 
u[[2,2]]=-0.25*h*f[[2,2]]
},
{
For[i=1,i≤mu1,i++,u=relax[u,f]];

r=f2c[findRes[u,f]];
{n,n}=Dimensions[r];
Print["in vcycle, after f2c, n=",n];
correct = vcycle[Table[0,{i,n},{j,n}],r,mu1,mu2];
correct = c2f[correct];
Print["in vcycle, after f2c, correct=",correct];
u=u+correct;
For[i=1,i≤mu2,i++,u=relax[u,f]];
}
];
u
];
Attributes[vcycle]={HoldFirst};

h=N[2^-5];
n=(1/h)+1;
xcoord=Table[(j-1)*h,{i,1,n},{j,1,n}];
ycoord=Table[(n-i)*h,{i,1,n},{j,1,n}];
f=Table[force[{xcoord[[i,j]],ycoord[[i,j]]}],{i,1,n},{j,1,n}];
grid=Flatten[Table[{xcoord[[i,j]],ycoord[[i,j]],f[[i,j]]},{i,1,n},{j,1,n}],1];
tol=0.01*h^2;
mu1=1;
mu2=1;
done=False;
k=0;
normF=Norm[f,2];
Print["normF=",normF];
Print["n=",n];
u=Table[0,{i,n},{i,n}];

While[Not[done],
{
k++;
Print["k=",k];
u=vcycle[u,f,mu1,mu2];
res=findRes[u,f];
relRes = Norm[res,2]/normF;
If[relRes<tol,done=True];
Print["relRes=",relRes];
If[k==10,done=True];
}
]