

function [Z,D,it,err]=simit(A,Z0=eye(size(A)),tolr=1e-12,tola=tolr*1e-4,maxit=200)
#simultanuous iteration - finding eigenvectors and eigenvalues of symmetric A
# Z_{old}- orthogonal then next orthogonal  Z_{new} is such that A*Z_{old}=Z_{new}R
#where R upper triangular 
#stoping criterium - offdiagonal elements of Z'AZ almost zero i.e.
#if D_k- offdiagonal elements of 
# D_k=Z_k^TAZ_k then    stopping criterium : err(k)=\|D_k- diagonal(D_k)\|<tol for 
#tol=tor*err(1)+tola 
#Input 
#a - symmetric matrix
#Optional arguments
#Z0 - orthogonal matrix (default value Z0=I identity)
#tolr,tola - tolerences - default values tolr=1e-12 and tola=tolr*1e-4
#maxit- default = 200
#
#Output:
#Z - orthogonal matrix with columns as eigenvectors
#D - diagonal matrix with eigenvalues on diagonal
#err - vector with the history of convergence i.e. norms \|D_k- diagonal(D_k)\|


if(norm(A'-A,1)>1e-12)
 printf("Matrix A is NOT symmetric! Terminating\n");
 return;
endif

Z=Z0;
it=0;
err=zeros(maxit,1);
D=Z'*A*Z;
er=err(it+1)=norm(D-diag(diag(D)),'fro')
tol=tolr*er+tola;

while((er>tol)&&(it<maxit))
 Y=A*Z;
 [Z,R]=qr(Y);
 D=Z'*A*Z;
 it++;
 er=err(it+1)=norm(D-diag(diag(D)),'fro');
endwhile
 err=err(1:it+1);
endfunction
