function[vysledek,iterace]=new(A,M,N,beta,presnost,maxit)
%na vstupu: A - daná matice
%           M,N - symetrické pozitivně definitní matice
%           beta - konstanta
%           presnost - přesnost algoritmu
%           maxit - maximální počet iterací
%na výstupu: vysledek - hledaná matice
%            iterace - konečný počet iterací
if nargin < 6
    maxit=5000;
    if nargin < 5
        presnost=1e-6;
    end;
end;
[m,n]=size(A);
r=rank(A);
[m1,m2]=size(M);
if norm(M-M')~=0 | m1~=m2 | m1~=m | sum(eig(M)>0)~=m
   error(['M musi byt symetricka pozitivne ', ...
        'definitni matice typu m x m'])
end;
[n1,n2]=size(N);
if norm(N-N')~=0 | n1~=n2 | n1~=n | sum(eig(N)>0)~=n
   error(['N musi byt symetricka pozitivne ', ...
        'definitni matice typu n x n'])
end;
Akriz=inv(N)*A'*M;
vlcisla=eig(Akriz*A);
lambda=max(vlcisla);
betamax=2/(lambda.^2);
if beta <= 0 | beta >= betamax
    beta=1/(lambda.^2);
end;
X=beta*Akriz;
k=0;
t=maxit;
I=eye(m);
while k<maxit
    psi=A*X;
    ksi=-7*I+psi*(9*I+psi*(-5*I+psi));
    omega=psi*ksi;
    Xplus=-(1/16)*X*ksi*(4*I+omega)*(8*I+omega*(4*I+omega));
    X=Xplus;
    k=k+1;
    n1=norm(A*X*A-A);
    n2=norm(X*A*X-X);
    n3=norm(M*A*X-(M*A*X)');
    n4=norm(N*X*A-(N*X*A)');
    if n1<=presnost && n2<=presnost ... 
       && n3<=presnost && n4<=presnost
        t=k;
        k=maxit;
    end;
end;
iterace=t;
vysledek=X;

%rand(100,120) - 8 iteraci (iterative3 dava 14)
%rand(1000,1200) - 11