precond_laplace_noorth

   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 function [eta] = precond_laplace_noorth( L, xi, xL, xR, G )
0006     
0007     r = xi.rank;
0008     n = xi.size;
0009     d = xi.order;
0010 
0011     eta = xi;
0012     xi = tangent_to_TTeMPS( xi );
0013 
0014     % 1. STEP: Project right hand side
0015 
0016     Y = cell(1,d);
0017     % precompute inner products
0018     left = innerprod( xL, xi, 'LR', d-1, true );             
0019     right = innerprod( xR, xi, 'RL', 2, true );             
0020 
0021     % contract to first core
0022     Y{1} = tensorprod_ttemps( xi.U{1}, right{2}, 3 );      
0023     % contract to first core
0024     for idx = 2:d-1
0025         res = tensorprod_ttemps( xi.U{idx}, left{idx-1}, 1 );
0026         Y{idx} = tensorprod_ttemps( res, right{idx+1}, 3 );
0027     end 
0028     % contract to last core
0029     Y{d} = tensorprod_ttemps( xi.U{d}, left{d-1}, 1 );
0030 
0031     % 2. STEP: Solve ALS systems:
0032     % Instead of doing
0033     %    X_mid = orthogonalize(xR, idx);
0034     % we recursively adjust the gauge based on xL and xR
0035     X_mid = xR;
0036     eta.dU{1} = solve_inner( L{1}, X_mid, Y{1}, 1 );
0037     for idx = 2:d
0038         X_mid.U{idx-1} = xL.U{idx-1};
0039         X_mid.U{idx} = tensorprod_ttemps(X_mid.U{idx},G{idx-1},1);
0040         eta.dU{idx} = solve_inner( L{idx}, X_mid, Y{idx}, idx );  
0041     end
0042 
0043     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
0044 
0045 end
0046 
0047 function X = get_mid(xL, xR, G, idx)
0048 X = xR;
0049 X.U{1:idx-1} = xL.U{1:idx-1};
0050 if idx>1
0051     X.U{idx} = tensorprod_ttemps(X.U{idx},G{idx-1},1);
0052 end
0053 end
0054 
0055 function [res,BB1,BB3] = solve_inner( L0, X, Fi, idx )
0056     n = size(L0, 1);
0057     rl = X.rank(idx);
0058     rr = X.rank(idx+1);
0059 
0060     B1 = zeros( rl );
0061     %BB1 = {};
0062     % calculate B1 part:
0063     for i = 1:idx-1
0064         % apply L to the i'th core
0065         tmp = X;
0066         tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
0067         %BB1{i} = tmp.U{i};
0068         B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
0069     end
0070 
0071     % calculate B2 part:
0072     B2 = L0;
0073 
0074     B3 = zeros( rr );
0075     %BB3 = {};
0076     % calculate B3 part:
0077     for i = idx+1:X.order
0078         tmp = X;
0079         tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
0080         %BB3{i} = innerprod( X, tmp, 'RL', idx+1);
0081         B3 = B3 + innerprod( X, tmp, 'RL', idx+1);
0082     end
0083 
0084     % Faster below
0085     %[V,E] = eig( kron( eye(rr), B1 ) + kron( B3, eye(rl) ) );
0086     %E = diag(E);
0087             
0088     [V1,E1] = eig(B1); [V3,E3] = eig(B3);
0089     V = kron(V3,V1);    
0090     EE = diag(E1)*ones(1,rr) + ones(rl,1)*diag(E3)'; E = EE(:);
0091     
0092     rhs = matricize( Fi, 2 ) * V;
0093     Y = zeros(size(rhs));
0094     for i=1:length(E)        
0095         Y(:,i) = (B2 + E(i)*speye(n)) \ rhs(:,i);
0096     end
0097     res = tensorize( Y*V', 2, [rl, n, rr] );
0098 end
0099