function [u,k] = msg(A,b,M,epsilon)

% [u,k] = msg(A,b,M,epsilon)
%         Funkce resi soustavu rovnic A*u = b
%         metodou sdruzenych gradientu s predpodminenim.
%         Matice M je vhodny predpodminovac.
%         epsilon je rel. presnost reseni ve smyslu norm(A*u-b).
%         Algoritmus je opsan z knihy
%         Golub, van Loan - Matrix Computation, str.529.

k = 0;
u = 0;
rk = b;

while norm(rk) > eps*norm(b)
    zk = M\rk;

    if k==0
        p = zk;
    else
        beta = (rk'*zk) / (rk_1'*zk_1);
        p = zk + beta*p;
    end
    
    alpha = (rk'*zk)/(p'*A*p);
    u = u + alpha*p;
    
    rk_1 = rk;
    rk = rk - alpha*A*p;
    
    zk_1 = zk;
    
    k = k+1;
    norm(rk);
end
