function Lap = Konj_laplace(n)
% resime ve 2D ulohu Laplace na phi = 0 s okrajovými podmínkami
% phi = 0 pro y = 0;  phi = sin(x)*sinh(pi/2)  pro y = pi/2
% phi = 0 pro x = 0;  phi = sinh(y) pro x = pi/2
% Laplace phi v diferencnim schematu => 
% (phi(x+h,y) + phi(x-h,y) - 2 phi(x,y))/h^2 + (phi(x,y+h) + phi(x,y-h) - 2 phi(x,y))/h^2    
% dx = dy = pi/2/(n+1)
% n pocet vnitrnich bodu v kazdem smeru, celkem n^2 vnitrnich bodu, n^2 rovnic
% matice symetricka, ridka, tridiagonalni + 2 symetricke diagonaly
% prikad nevyuziva ridkosti matice, ale ukazuje jak rychle konverguje 
% Konjug gradient  pro n = 4 (16 rovnic) = 9 iteraci
% Konjug gradient  pro n = 5 (25 rovnic) = 13 iteraci 
n2 = n^2;
for i = 1:n2
    c(i) = 0.;
    d(i) = 0.;
    for j = 1:n2
       a(i,j) = 0.;
    end     
end
for i = 1:n2
    a(i,i) = 4;
end
n2mn = n2-n;
for i = 1:n2mn
  a(i,i+n) = -1;
end
for i = n+1:n2
  a(i,i-n) = -1;
end
for j = 0:n-1
   for i = 2:n
       ij = n*j+i;
       a(ij,ij-1) = -1;
   end
end   
for j = 0:n-1
   for i = 1:n-1
       ij = n*j+i;
       a(ij,ij+1) = -1;
   end
end   
a
dx = 0.5*pi/(n+1);
for i= 1:n
    x = i*dx;
    is = n2mn +i;
    c(is) = sin(x)*sinh(pi);
end;
for j = 1:n
    y = j*dx; 
    ir = n*j;
    c(ir) = sinh(y) + c(ir);
end 
b = c'
v1 = min(eig(a));
fprintf('min vlastni cislo %g > 0 => pozitivne definitni \n',v1)
bnorm = norm(b); 	
prec = 1.e-6;
fprintf('zadana presnost %g \n',prec)
phinew = d';
hold = d';
it = 0;
gnew = b - a*phinew;
nn = norm(gnew);
fprintf('norma (b - a phi) = %g pro pocatecni odhad phi \n',n)
while nn > prec*bnorm 
   if norm(hold) == 0
       hnew = gnew;
   else
       cit = -gnew'*(a*hold);
       jm = hold'*(a*hold);
       gam = cit/jm;
       hnew = gnew + gam*hold; 
   end
   cit = gnew'*hnew;
   jm = hnew'*(a*hnew);
   alf = cit/jm;
   phiold = phinew;
   it = it + 1;
   phinew = phiold + alf*hnew;
   hold = hnew;
   gnew = b - a*phinew;
   nn = norm(gnew);
   fprintf('Po %d iteraci je norma (b - a *phi) = %g \n',it,nn) 
end
   fprintf('Reseni nalezeno po %d iteracich \n',it) 
   phi = inv(a)*b;
   for i = 1:n2
      fprintf('i  %d  phi_iter(i) %g  phi_res %g \n',i,phinew(i),phi(i))
   end
