function [X,itrec]=ttm_inv2(t,eps0,reduce,warnings,eps,eps_H,mode,recursive)
%% Thomas Mach
% inverting a ttm
% copied from TT Toolbox (Version 1)
% modified to 0seST09 equation (18)
niters = 100; % 100
factor = 1e6;
itrec = 0;

if (~exist('reduce','var'))
    reduce=0;
end
if (~exist('eps','var'))
    eps=1e-15;
end
if (~exist('eps_H','var'))
    eps_H=1e-15;
end
if (~exist('eps0','var'))
    eps0=1e-6;
end
if (~exist('warnings','var'))
    warnings=0;
end
if (~exist('mode','var'))
    mode=4;
end
if (~exist('recursive','var'))
    recursive=3;
end
vm = t.m;
vn = t.n;
tt1=t.tt;
d = tt1.d;
% d = length(vn);

if (warnings>1)
  fprintf('niters %d factor %e eps0 %e reduce %d\n',niters, ...
	  factor, eps0, reduce);
  
  fprintf('warnings %d eps %e eps_H %e mode %d\n', ...
	  warnings, eps, eps_H, mode);
end

if ((eps0/sum(vn) <eps *1e5)&&warnings)
    fprintf('Warning, eps/0 seems to be too small\n');
end

et = tt_matrix(tt_eye(vn,d));
alpha = 2/norm(t);
e2=2*et;

i=1;
X=alpha*et;
% Xo1=alpha*et;
% Xo2=alpha*et;
Y=alpha*t;
H= Y;
rs=200000;
rs2=400000;
comprs=rs;
Xmin = X;
Ymin = Y;
rmin = rs;
epsmin = eps;
eps_Hmin = eps_H;
imin = i;
Xmino = Xmin;
Ymino = Ymin;
rmino = rmin;
epsmino = epsmin;
eps_Hmino = eps_Hmin;
imino = imin;
told=toc;
while(i<niters && comprs > eps0 )
  if ((mode == 1)||(mode == 3)||(mode == 4)||(mode == 6))
    t1= toc;
    if (eps_H<1e-31)
      H = round(e2 - Y,eps_H);
    else
      H = e2 - Y;
    end
    t2= toc;
    Y = round(Y*H,eps);
    t3= toc;
%    Xo2=Xo1;
%    Xo1=X;
    X = round(X*H,eps);
    t4= toc;
    rs= norm(Y-et);
    t5= toc;
    if (warnings>2)
      fprintf('t21 %10.7f, t32 %10.7f, t43 %10.7f, t54 %10.7f\n', t2-t1,t3- ...
	      t2,t4-t3,t5-t4);
    end
  end
  if (mode == 2)
%    X = round(2*X-X*t*X,eps);
%    X = round(2*X-X*t*X,eps);
%    X = round(round(e2-X*t,eps)*X,eps);
    X = round(X*round(e2-t*X,eps),eps);
    rs = norm(et-X*t);
  end
  if (mode == 5)
    t1= toc;
    H = e2 - Y;
    t2= toc;
    t3= toc;
    X = round(X*H,eps);
    t4= toc;
    t5= toc;
    if (warnings>2)
      fprintf('t21 %10.7f, t32 %10.7f, t43 %10.7f, t54 %10.7f\n', t2-t1,t3- ...
	      t2,t4-t3,t5-t4);
    end
  end
  if (warnings>1)
    tnew = toc;
    fprintf('iter: %2d, res=%e, mem(X)=%e, mem(Y)=%e, ',i,rs,mem(X),mem(Y));
    fprintf('mem(H)=%e, time=%10.7f ttime=%10.7f\n',mem(H),tnew-told,tnew);
    told = tnew;
    if (mode ~= 6)
      if (rs<eps0)
	mode=5;
      end
      if ((mod(i,5)==0)||(mode==5))
	Z = round(t*X,eps);
	comprs = norm(et-Z);
	fprintf('compare   res=%e\n',comprs);
	if (abs(comprs-rs)/rs>1.2)
	  mode=5;
	  fprintf('update Y\n');
	  Y = Z;
	  rs = comprs;
	  tnew2 = toc;
	  eps=10*eps;eps_H=10*eps_H;
	  fprintf('time for update %3.7f\n',tnew2-told);
	end
      end
    else 
      comprs=rs;
    end
  end
  i=i+1;
  
  if ((mode==1)&&(rs/rs2>10))
    mode = 2;
    X = Xmin;
    eps = epsmin;
    eps_H = eps_Hmin;
    i = imin-1;
    rs = rmin;
    fprintf('change to iteration mode 2\n');
    if (warnings>1)
      fprintf('iter: %d, res=%e, mem=%e \n',i,rs,mem(X));
    end
    i=i+1;
  else
    if (reduce)
      if (abs(rs2-rs)<factor*eps)
	eps=eps/10;
	eps_H=eps_H/10;
	if (warnings)
	  fprintf('reduce eps to %e\n',eps);
	end
      else
	if ((rs-rs2>0)&&(mode==3))
	  X = Xmino;
	  Y = Ymino;
	  %eps = epsmino;
	  %eps_H = eps_Hmino;
	  i = imino-1;
	  rs = rmino;
	  factor=10*factor;
	  eps=eps/100;
	  eps_H=eps_H/100;
	  if (warnings)
	    fprintf('reduce eps to %e and increase factor %e\n',eps,factor);
	  end
	  if (warnings>1)
	    fprintf('iter: %d, res=%e, mem=%e \n',i,rs,mem(X));
	  end
	  i=i+1;

	end
      end
    end
  end
  rs2=rs;
  
  if ((mode~=2)&&(rs<rmin))
    Xmino = Xmin;
    Ymino = Ymin;
    rmino = rmin;
    epsmino = epsmin;
    eps_Hmino = eps_Hmin;
    imino = imin;
    Xmin = X;
    Ymin = Y;
    rmin = rs;
    epsmin = eps;
    eps_Hmin = eps_H;
    imin = i;
  end
  
end
if ((recursive>=1)&&((abs(rs)>100)||((~(rs<100))&&(~(rs>=100)))))
  fprintf('Did not converge! try again.\n');
  eps=eps/10;
  [X,itrec]=ttm_inv2(t,eps0,reduce,warnings,eps,1e-100,mode,recursive-1);
  itrec=itrec+1;
  
else
  if (rs>10)
    fprintf('Did not converge! return tt_eye.\n');
    X=et;
  end
  
  if (~(rs<10)&&~(rs>=10))
    fprintf('Did not converge! return tt_eye.\n');
    X=et;
  end
end

if (warnings>2)
  Z = round(t*X,eps);
  comprs = norm(et-Z);
  fprintf('compare   res=%e %e\n',comprs, norm(et-t*X));
  fprintf('||I-XM||_F*sqrt{||M||_F}=%e\n',comprs*sqrt(norm(t)));
  if (warnings>3)
    sqrtnormt=sqrt(norm(t,2));
    normetmZ=norm(et-Z,2);
    fprintf('||I-XM||_2*sqrt{||M||_2}=%e = %e*%e\n',normetmZ*sqrtnormt,normetmZ,sqrtnormt);
  end
end

return

