nh=30;
T=diag(2*ones(nh,1))+diag(-ones(nh-1,1),1)+diag(-ones(nh-1,1),-1);
I=speye(nh);
A=-(kron(T,I)+kron(I,T));
n=nh^2;
B=eye(n,30);
m=1000;
tol=1e-9;
opts.tol=1e-2;
l1=eigs(-(A+A')/2,1,'sa',opts)
l2=eigs(-(A+A')/2,1,'la',opts)
[V,Z,js,s]=rksm_tan_lyap_amg_Mdir(A,l1,l2,B,m,tol,0,0,1);

fprintf('final true absolute residual norm: \n')
disp(norm(A*Z*Z'+Z*Z'*A'+B*B'))    %this matrix should never be formed for n large 

