fprintf('\n\nSPINVIT ohne shift\n');

if (~exist('tt_matrix','file')) % try to load the tensor train toolbox
  fprintf('Tensor Train Toolbox required!');
  addpath('~/matlab/TT-Toolbox/'); % change the path
  addpath('~/matlab/TT-Toolbox/core/');
  addpath('~/matlab/TT-Toolbox/exp/');
end
assert(exist('tt_matrix','file')==2); % test if tensor train toolbox is available

% some constants

%% Parameters for all examples 
dmax2 = 5;
epsinv= 1e-6;

%% Table 1
 
dd = 2;  % dimension of the Laplace-Operator, exact ev for dd<=4
dmin = 2; % 2^dmin ... 2^dmax discretization points per direction
dmax = 16; % 
eps_be = 1e-2; % truncation accuarcy for the first iteration
be = 4; % reduce the truncation accuracy by a 
        % factor of 10 for the first be steps
eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
c=0.2; % preconditioner accuracy 
s = 1+dd; % dimension of subspace
niters = 100;
display_status = 0;  % 1 some output; 2 more output
sh = 0; % 0 kein shift, 1 mit shift

%% Table 2

% dd = 2;  % dimension of the Laplace-Operator, exact ev for dd<=4
% dmin = 2; % 2^dmin ... 2^dmax discretization points per direction
% dmax = 8; % 
% eps_be = 1e-2; % truncation accuarcy for the first iteration
% be = 4; % reduce the truncation accuracy by a 
%         % factor of 10 for the first be steps
% eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
% c=0.20; % preconditioner accuracy 
% s = 1+dd; % dimension of subspace
% niters = 100;
% display_status = 0;  % 1 some output; 2 more output
% sh = 1; % 0 kein shift, 1 mit shift

%% Table 3

% dd = 3;  % dimension of the Laplace-Operator, exact ev for dd<=4
% dmin = 2; % 2^dmin ... 2^dmax discretization points per direction
% dmax = 10; % 
% eps_be = 1e-2; % truncation accuarcy for the first iteration
% be = 4; % reduce the truncation accuracy by a 
%         % factor of 10 for the first be steps
% eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
% c=0.20; % preconditioner accuracy 
% s = 1+dd; % dimension of subspace
% niters = 100;
% display_status = 0;  % 1 some output; 2 more output
% sh = 0; % 0 kein shift, 1 mit shift

%% Table 4
 
% dd = 3;  % dimension of the Laplace-Operator, exact ev for dd<=4
% dmin = 2; % 2^dmin ... 2^dmax discretization points per direction
% dmax = 5; % 
% eps_be = 1e-2; % truncation accuarcy for the first iteration
% be = 4; % reduce the truncation accuracy by a 
%         % factor of 10 for the first be steps
% eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
% c=0.20; % preconditioner accuracy 
% s = 1+dd; % dimension of subspace
% niters = 100;
% display_status = 0;  % 1 some output; 2 more output
% sh = 1; % 0 kein shift, 1 mit shift
  
%% Table 5

% dd = 4;  % dimension of the Laplace-Operator, exact ev for dd<=4
% dmin = 2; % 2^dmin ... 2^dmax discretization points per direction
% dmax = 7; % 
% eps_be = 1e-2; % truncation accuarcy for the first iteration
% be = 4; % reduce the truncation accuracy by a 
%         % factor of 10 for the first be steps
% eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
% c=0.20; % preconditioner accuracy 
% s = 1+dd; % dimension of subspace
% niters = 100;
% display_status = 0;  % 1 some output; 2 more output
% sh = 0; % 0 kein shift, 1 mit shift

%% Example 5.4 from Tobler/Kressner 2011

% dd = -1;  % dimension of the Laplace-Operator, exact ev for dd<=4
% dmin = 4; % 2^dmin ... 2^dmax discretization points per direction
% dmax = 32; % 
% eps_be = 1e-2; % truncation accuarcy for the first iteration
% be = 5; % reduce the truncation accuracy by a 
%         % factor of 10 for the first be steps
% eps_eig = 1e-5; % stop iteration if residual falls below eps_eig
% c=0.20; % preconditioner accuracy 
% s = 1+dd; % dimension of subspace
% s = 5;
% niters = 100;
% display_status = 0;  % 1 some output; 2 more output
% sh = 0; % 0 kein shift, 1 mit shift

%% Shifts
if (sh) % shift and subspace dimension suitable to dd, e.g.,
    if (dd==2)
        shift = 203.3139*ones(dmax,1);
        s = 3; % dimension of suspace for SPINVIT
    end
    if (dd==3)
        shift = 230.6195*ones(dmax,1);
        s = 6; % dimension of suspace for SPINVIT
    end
end
% variables for the measurements
outd=zeros(dmax,1);
outn=zeros(dmax,1);
outt=zeros(dmax,1);
outtinv=zeros(dmax,1);
outerr=zeros(dmax,1);
outerr2=zeros(dmax,1);
outiter=zeros(dmax,1);
outitrec=zeros(dmax,1);

d = dmin;
while (d<=dmax)
    outd(d)=d;
    n=2^d; %  One-dimensional grid size; full matrix will be n^2 x n^2 x ...
    outn(d)=n;
    fprintf('\nDimenision d = %d, 2^d = %d, n = %d\n',d,n,n^dd);

    if (dd>0)
        % compute ev by exact formula
        k=0;
        ev = zeros(s^dd,1);
        if (sh)
            if (dd==2)
                for evi=1:min(min(3,d)*s,n)
                    for evj=1:min(min(3,d)*s,n)
                        k=k+1;
                        ev(k) = 4*(sin(evi/(n+1)*pi/2)^2+sin(evj/(n+1)*pi/2)^2)*(n+1)^2;
                    end
                end
            end
            if (dd==3)
                for evi=1:min(min(3,d)*s,n)
                    for evj=1:min(min(3,d)*s,n)
                        for evk=1:min(min(3,d)*s,n)
                            k=k+1;
                            ev(k) = 4*(sin(evi/(n+1)*pi/2)^2+sin(evj/(n+1)*pi/2)^2+sin(evk/(n+1)*pi/2)^2)*(n+1)^2;
                        end
                    end
                end
            end
            [~,I] = sort(abs(ev-shift(d)));
            ev=ev(I);
        else
            if (dd==2)
                for evi=1:s
                    for evj=1:s
                        k=k+1;
                        ev(k) = 4*(sin(evi/(n+1)*pi/2)^2+sin(evj/(n+1)*pi/2)^2)*(n+1)^2;
                    end
                end
            end
            if (dd==3)
                for evi=1:s
                    for evj=1:s
                        for evk=1:s
                            k=k+1;
                            ev(k) = 4*(sin(evi/(n+1)*pi/2)^2+sin(evj/(n+1)*pi/2)^2+sin(evk/(n+1)*pi/2)^2)*(n+1)^2;
                        end
                    end
                end
            end
            if (dd==4)
                for evi=1:s
                    for evj=1:s
                        for evk=1:s
                            for evl=1:s
                                k=k+1;
                                ev(k) = 4*(sin(evi/(n+1)*pi/2)^2 + ...
                                    sin(evj/(n+1)*pi/2)^2 + ...
                                    sin(evk/(n+1)*pi/2)^2 + ...
                                    sin(evl/(n+1)*pi/2)^2)*(n+1)^2;
                            end
                        end
                    end
                end
            end
            ev = sort(ev);
        end
    end
    % Laplace by A. Kazeev  and B.N. Khoromskij
    if (dd==2)
        t=tt_qlaplace_dd([d,d]);
    end
    if (dd==3)
        t=tt_qlaplace_dd([d,d,d]);
    end
    if (dd==4)
        t=tt_qlaplace_dd([d,d,d,d]);
    end
    if (dd==-1)
      t = example43_tobler(d);
      outn(d)=1./outn(d);
    else
      t=tt_matrix(t); % Create tensor with dd-D Laplacian
% on 2^d \times 2^d grid with Dirichlet BC
      t=t*((n+1)^2); % scale tensor with 1/h^2
    end
    vm = t.m;
    vn = t.n;
    d2 = t.d;
    m = vm(1);
    n = vn(1);
    assert(vm(1)==vn(1));
    for i=2:d2
        m = m * vm(i);
        n = n * vn(i);
        assert(vm(i)==vn(i));
    end
    
    % identity tensor e
    e=cell(d2,1);
    e{1}=zeros(vm(1),vm(1),1);
    e{1}=eye(vm(1));
    for i=2:d2-1
        e{i}=zeros(vm(i),vm(i),1,1);
        e{i}(:,:,1,1)=eye(vm(i));
    end
    e{d2}=zeros(vm(d2),vm(d2),1);
    e{d2}(:,:,1)=eye(vm(d2));
    et = tt_matrix(e);
    
    vec = cell(s,1); % s TT
    mu = zeros(s,s); % sxs matrix for the gen Rayleigh Quotient
    % random vector generator
    matvec=cell(s,1); 
    scavec=cell(s,1);
    for j=1:s
        vec{j}=cell(d2,1);
        for i=1:d2
            vec{j}{i}=rand(vm(i),1);
        end
        vec{j} = tt_tensor(vec{j});
        matvec{j} = vec{j};
        scavec{j} = vec{j};
    end
    
    if (sh)
        msigma = round(t-shift(d)*et,1e-20);
        msigma2 = round(msigma*msigma,eps);
    else
        msigma2 = t;
    end
    % Inversion!
    tic
    [pre,itrec]=ttm_inv2(msigma2,c/sqrt(norm(msigma2,'fro')),0,2,epsinv,1e-100,6,4);
    outtinv(d)=toc; % measure time for inversion
    outitrec(d)=itrec;
    if (dd==-1)
      if (n<1024) ev=sort(eig(full(msigma2)));
      else ev = zeros(s,1);
      end
    end
    fprintf('tinv %10.6f\n',outtinv(d));
    if (display_status >= 1)
        input('inversion finished\n');
    end
    
    tic
    % orthonormalize
    for j=1:s
        for k=1:j-1
            beta = -dot(vec{k},vec{j});
            vec{j} = vec{j} - beta*vec{k};
        end
        vec{j} = vec{j}/norm(vec{j});
    end;
    
    % PINVIT
    % TTM-TT product
    for j=1:s
        matvec{j} = msigma2*vec{j};
    end
    % gen Rayleigh quotient
    for j=1:s
        for k=1:s
            mu(j,k) = dot(vec{j},matvec{k});
        end
    end
    if (display_status >= 2)
        disp(mu);
    end
    
    for j=1:s
        scavec{j} = (-mu(j,1)) * vec{1};
        for k=2:s
            scavec{j} = scavec{j} - (mu(k,j)*vec{k});
        end
        scavec{j} = round(scavec{j},eps);
    end
    
    % residual
    for j=1:s
        matvec{j} = matvec{j} + scavec{j}; 
    end
    alpha = 0;
    for j=1:s
        alpha = alpha + dot(matvec{j},matvec{j});
    end
    iter = 0;
    eps_v = eps_be;
    alpha2 = nan;
    % Iteration
    alpha
    while ((alpha > eps_eig)&&(iter<niters))
        iter = iter+1;
        for j=1:s
            vec{j} = round(vec{j} - round(pre*matvec{j},eps_v),eps_v); % update
        end
        
        % orthonormalize
        for j=1:s
            for k=1:j-1
                alpha = -dot(vec{k},vec{j});
                vec{j} = round( vec{j} + alpha*vec{k} ,eps_v);
            end
            vec{j} = round( vec{j}/norm(vec{j}) ,eps_v);
        end;
  
        % TTM-TT product
        for j=1:s
            matvec{j} = round( msigma2*vec{j} ,eps_v);
        end
        % gen Rayleigh quotient
        for j=1:s
            for k=1:s
                mu(j,k) = dot(vec{j},matvec{k});
            end
        end
        if (display_status >= 2)
            disp(mu);
        end
        
        % residual
        for j=1:s
            scavec{j} = round( (-mu(1,j)) * vec{1} ,eps_v);
            for k=2:s
                scavec{j} = round( scavec{j} + (-mu(k,j)) * vec{k} ,eps_v);
            end
        end
        
        for j=1:s
            matvec{j} = round( matvec{j} + scavec{j} ,eps_v);
        end
        alpha = 0;
        for j=1:s
            alpha = alpha + dot(matvec{j},matvec{j});
        end
        if (display_status>1)
            disp(alpha);
        end
        fprintf('iter: %2d, alpha=%e, mem=%e, %e, %e, eps_v=%e \n',iter,alpha, ...
            mem(vec{1}),mem(vec{2}),mem(vec{3}), eps_v);
        if (abs(alpha-alpha2)/alpha <1e-2)
            eps_v=eps_v/10;
            fprintf('eps reduced\n');
        end
	if (alpha<100*eps_v)
            eps_v=eps_v/10;
            fprintf('eps reduced\n');
        end
        if (iter<be)
            eps_v=eps_v/10;
        end
        alpha2=alpha;
    end  % end while
    outiter(d)=iter;
    clear mu;
    % gen Rayleigh quotient
    for j=1:s
        for k=1:j-1
            alpha = dot(vec{k},vec{j});
            vec{j} = round( vec{j} + alpha*vec{k}, eps_v);
        end
        vec{j} = round( vec{j}/norm(vec{j}) ,eps_v);
        matvec{j} = round( t*vec{j}, eps_v);
        mu(j) = dot(vec{j},matvec{j});
    end;
    
    for j=1:s
        for k=1:s
            mu(j,k) = dot(vec{j},matvec{k});
        end
    end
    disp(mu);
    eigmu = eig(mu);
    if (sh)
        [~,I] = sort(abs(eigmu-shift(d)));
        eigmu=eigmu(I);
    else
        eigmu=sort(eigmu);
    end
    disp([eigmu,ev(1:s)]);
    disp(eigmu-ev(1:s));
    outt(d)=toc;
    disp(outt(d));
    outerr(d) = norm(ev(1:s)-eigmu,'inf');
    outerr2(d) = norm(ev(1:s)-eigmu,2);
    disp(outerr(d));
    disp(outerr2(d));
    % fprintf('memory pre = %e\n',t_mem(pre)); % memory consumption
    % fprintf('memory mat = %e\n',tt_mem(mat));
    
    if (display_status >= 1)
        input('next d\n');
    end
    
    fprintf('%% %d & %d & %10.3f & %10.3f & %d & %.4e \\\\ \\hline \n',outn(d)^dd,outd(d),outtinv(d),outt(d),outiter(d),outerr(d));
    if (dd==-1)
      d = 2*d;
    else
      d = d+1;
    end
    if (d>dmax2)
      fprintf('%% dd %d, eps %.2e, eps_be %.2e, be %d, eps_eig %.2e, c %.2e\n', dd, eps, eps_be, be, eps_eig, c);
      fprintf('    $n$     & $d$ & $t_{\\text{inv}}$ in s & $t_{\\text{PINVIT}}$ in s & \\# it & error      \\\\ \\hline\n')
      for de=dmin:length(outd)
	fprintf(' %10d & %3d & %18.3f %2d & %24.3f & %5d & %.4e \\\\ \\hline \n',outn(de)^dd,outd(de),outtinv(de),outitrec(de),outt(de),outiter(de),outerr(de));
      end
    end      
end
% Output as LaTeX Source Code (tabular)
fprintf('%% dd %d, eps %.2e, eps_be %.2e, be %d, eps_eig %.2e, c %.2e\n', dd, eps, eps_be, be, eps_eig, c);
fprintf('    $n$     & $d$ & $t_{\\text{inv}}$ in s & $t_{\\text{PINVIT}}$ in s & \\# it & error      \\\\ \\hline\n')
for d=dmin:length(outd)
    fprintf(' %10d & %3d & %18.3f %2d & %24.3f & %5d & %.4e \\\\ \\hline \n',outn(d)^dd,outd(d),outtinv(d),outitrec(d),outt(d),outiter(d),outerr(d));
end
