function mg2d()
k1=3;
k2=3;
gamma=2;
ciklusszam=4;

% Megoldandó feladat jobb oldala
fexact=@(x,y)sin(pi*x).*sin(pi*y)*2*pi^2;

% Pontos megoldás
uexact=@(x,y)sin(pi*x).*sin(pi*y); 

p0=5; % a legfinomabb rács paramétere
n=2^p0-1;
h=1/(n+1);
A=gallery('poisson',n)/h^2;
felosztas=h:h:1-h;
[X,Y]=meshgrid(felosztas,felosztas);
uexactm=uexact(X,Y);
disp('A pontos megoldás grafikonja')
mesh(X,Y,uexactm);
pause
vm=rand(n,n);
vv=reshape(vm,n^2,1);
disp('Egy véletlen mátrixról indulunk')
mesh(X,Y,vm);
pause;
omega=2/3;

for j=1:ciklusszam
fprintf('%1.0f. CIKLUS\n',j) 
p=p0;
n=2^p-1;
h=1/(n+1);
A=gallery('poisson',n)/h^2;
phim=fexact(X,Y);
phiv=reshape(phim,n^2,1);

vv=multigrid2D(p,vv,phiv,k1,k2,gamma);

fprintf('A %1.0f. ciklus utáni grafikon\n',j)
mesh(X,Y,reshape(vv,n,n));
pause
end

function u=multigrid2D(p,u,jobboldal,k1,k2,gamma)
% Multigrid iteráció
% 2^p-1: belsõ pont van minkét irányban
% u: közelítõ megoldás - innét indul az iteráció 
% jobboldal: az egyenletrendszer jobb oldali vektora
% k1 : elõsimítások
% k2 : utósimítások
% gamma: ismétlések száma

fprintf('%1.0f. rétegen vagyunk\n',p)
n=2^p-1;
h=1/(n+1);
A=gallery('poisson',n)/h^2;
if p==1
    disp('Pontos megoldás a legalsó rétegen')
    u=A\jobboldal;
else
    
fprintf('%1.0f db elõsimítás\n',k1)
u=jacobirelax(A,jobboldal,k1,u,2/3);
disp('Maradékvektor-számítás, maradékegyenlet megoldása a MG algoritmussal')
r=jobboldal-A*u;
e=zeros((2^(p-1)-1)^2,1);
for i=1:gamma
    e=multigrid2D(p-1,e,reshape(restrikcio2D(reshape(r,2^(p)-1,2^(p)-1)),(2^(p-1)-1)^2,1),k1,k2,gamma);
end;
fprintf('%1.0f. rétegen vagyunk\n',p)
disp('Durva rács korrekció')
u=u+reshape(prolongacio2D(reshape(e,2^(p-1)-1,2^(p-1)-1)),n^2,1);
fprintf('%1.0f db utósimítás\n',k2)
u=jacobirelax(A,jobboldal,k2,u,2/3);
end

function [umfinom]=prolongacio2D(umdurva)
% Lineáris interpolációs prolongáció
meret=size(umdurva);
Z=zeros(meret(1)+2);
Z(2:end-1,2:end-1)=umdurva;
[X,Y]=meshgrid(1/2:1/2:meret(1)+1/2,meret(1)+1/2:-1/2:1/2);
umfinom=interp2(0:meret(1)+1,meret(1)+1:-1:0,Z,X,Y);

function [umdurva]=restrikcio2D(umfinom)
% Részleges restrikció
umdurva=(4*umfinom(2:2:end-1,2:2:end-1)+umfinom(1:2:end-2,2:2:end-1)+...
umfinom(3:2:end,2:2:end-1)+umfinom(2:2:end-1,3:2:end)+umfinom(2:2:end-1,1:2:end-2))/8;

function u=jacobirelax(A,f,klepes,u0,omega)
% klepes Jacobi-relaxációt csinál a u0 vektorról indulva.
% A az egyenletrendszer mátrixa, b pedig a jobb oldala.
u=u0; 
for i = 1:klepes
    u = u -omega*(A*u-f)./diag(A);
end  

