precond_laplace_overlapJacobi

   TTeMPS Toolbox. 
   Michael Steinlechner, 2013-2016
   Questions and contact: michael.steinlechner@epfl.ch
   BSD 2-clause license, see LICENSE.txt

0001 %   TTeMPS Toolbox.
0002 %   Michael Steinlechner, 2013-2016
0003 %   Questions and contact: michael.steinlechner@epfl.ch
0004 %   BSD 2-clause license, see LICENSE.txt
0005 
0006 function [eta, B1,B3] = precond_laplace_overlapJacobi( L, xi, xL, xR, G, B1, B3 )
0007 % L is a cell of operators
0008 
0009 r = xi.rank;
0010 n = xi.size;
0011 d = xi.order;
0012 
0013 % If B1 and B3 are not given as arguments, we need to precalculate them
0014 if nargin < 7
0015 %     % if applying L is expensive (not just tridiag), one can store all
0016 %     % applications with xL and compute the ones for xR with G.
0017 %     % You need to first store LUl
0018 %     LUl = cell(d,1);
0019 %     for idx = 1:d
0020 %         LUl{idx} = tensorprod_ttemps( xL.U{idx}, L{idx}, 2 );
0021 %     end
0022 %     % and then change to LUr in the loop for B3 below
0023 %     %         if idx+1==d
0024 %     %              LUr = tensorprod_ttemps( LUl{idx+1}, G{idx}, 1, true);
0025 %     %         else
0026 %     %             LUr = tensorprod_ttemps( tensorprod_ttemps( LUl{idx+1}, G{idx+1}', 3), G{idx}, 1, true);
0027 %     %         end
0028     
0029     B1 = cell(d,1);
0030     B1{1} = 0;
0031     for idx = 2:d
0032         LUl = tensorprod_ttemps( xL.U{idx-1}, L{idx-1}, 2 );
0033         if idx>2
0034             TT = tensorprod_ttemps( xL.U{idx-1}, B1{idx-1}, 1 );
0035         else
0036             TT = 0;
0037         end
0038         B1{idx} = unfold(xL.U{idx-1},'left')'*unfold(TT + LUl,'left');
0039     end
0040 
0041     B3 = cell(d,1);
0042     for idx = d-1:-1:1
0043         LUr = tensorprod_ttemps( xR.U{idx+1}, L{idx+1}, 2 );
0044         if idx<d-1
0045             TT = tensorprod_ttemps( xR.U{idx+1}, B3{idx+1}, 3 );
0046         else
0047             TT = 0;
0048         end          
0049         B3{idx} = unfold(xR.U{idx+1},'right')*unfold(TT + LUr,'right')';
0050     end
0051     B3{d} = 0;
0052 end
0053 
0054 eta = xi;
0055 xi = tangent_to_TTeMPS( xi );
0056 
0057 
0058 
0059 % % 1. STEP: Project right hand side
0060 % below is hard-coded version of
0061 % for ii=1:d
0062 %     eta_partial_ii = TTeMPS_partial_project_overlap( xL, xR, xi, ii);
0063 %     Y{ii} = eta_partial_ii.dU{ii};
0064 % end
0065 
0066 % TODO, it seems that the left and right cell arrays consist of a lot of
0067 % identities and zeros.
0068 Y = cell(1,d);
0069 % precompute inner products
0070 left = innerprod( xL, xi, 'LR', d-1, true );
0071 right = innerprod( xR, xi, 'RL', 2, true );
0072 
0073 % contract to first core
0074 Y{1} = tensorprod_ttemps( xi.U{1}, right{2}, 3 );
0075 % contract to first core
0076 for idx = 2:d-1
0077     res = tensorprod_ttemps( xi.U{idx}, left{idx-1}, 1 );
0078     Y{idx} = tensorprod_ttemps( res, right{idx+1}, 3 );
0079 end
0080 % contract to last core
0081 Y{d} = tensorprod_ttemps( xi.U{d}, left{d-1}, 1 );
0082 
0083 
0084 % 2. STEP: Solve ALS systems:
0085 
0086 % B1 and B3 were precalculated before
0087 for idx = 1:d
0088     rl = r(idx);
0089     rr = r(idx+1);
0090     
0091     B2 = L{idx};
0092   
0093     % Solve via the diagonalization trick
0094     [V1,E1] = eig(B1{idx}); [V3,E3] = eig(B3{idx});
0095     V = kron(V3,V1);
0096     EE = diag(E1)*ones(1,rr) + ones(rl,1)*diag(E3)'; E = EE(:);
0097     
0098     rhs = matricize( Y{idx}, 2 ) * V;
0099     Z = zeros(size(rhs));
0100     for i=1:length(E)
0101         Z(:,i) = (B2 + E(i)*speye(n(idx))) \ rhs(:,i);
0102     end
0103     eta.dU{idx} = tensorize( Z*V', 2, [rl, n(idx), rr] );
0104 end
0105 
0106 eta = TTeMPS_tangent_orth( xL, xR, eta );   % todo? Can we improve efficiency since eta is not a generic TTeMPS but shares the same x.U as xL and xR
0107 
0108 end
0109