0001
0002
0003
0004
0005
0006 function [eta,Yeta] = precond_laplace_overlapGS( L, xi, xL, xR, G, Lapl )
0007
0008 r = xi.rank;
0009 n = xi.size;
0010 d = xi.order;
0011
0012 eta = xi;
0013 xi = tangent_to_TTeMPS( xi );
0014
0015 for idx=1:d
0016 eta.dU{idx} = eta.dU{idx}*0;
0017 end
0018
0019
0020
0021
0022 Y = cell(1,d);
0023
0024 right = innerprod( xR, xi, 'RL', 2 );
0025 Y{1} = tensorprod_ttemps( xi.U{1}, right, 3 );
0026
0027 for idx = 2:d-1
0028 left = innerprod( xL, xi, 'LR', idx-1 );
0029 right = innerprod( xR, xi, 'RL', idx+1 );
0030 res = tensorprod_ttemps( xi.U{idx}, left, 1 );
0031 Y{idx} = tensorprod_ttemps( res, right, 3 );
0032 end
0033
0034
0035 left = innerprod( xL, xi, 'LR', d-1 );
0036 Y{d} = tensorprod_ttemps( xi.U{d}, left, 1 );
0037
0038
0039
0040
0041
0042 X_mid = xR;
0043 eta.dU{1} = solve_inner( L{1}, X_mid, Y{1}, 1 );
0044 for idx = 2:d-1
0045 X_mid.U{idx-1} = xL.U{idx-1};
0046 X_mid.U{idx} = tensorprod_ttemps(X_mid.U{idx},G{idx-1},1);
0047
0048 Eta = tangent_to_TTeMPS( eta );
0049
0050 PEta = xi - apply(Lapl, Eta);
0051 left = innerprod( xL, PEta, 'LR', idx-1 );
0052 right = innerprod( xR, PEta, 'RL', idx+1 );
0053 res = tensorprod_ttemps( PEta.U{idx}, left, 1 );
0054 Yeta = tensorprod_ttemps( res, right, 3 );
0055
0056
0057 eta.dU{idx} = solve_inner( L{idx}, X_mid, Yeta, idx );
0058 end
0059
0060 X_mid.U{d-1} = xL.U{d-1};
0061 X_mid.U{d} = tensorprod_ttemps(X_mid.U{d},G{d-1},1);
0062
0063 Eta = tangent_to_TTeMPS( eta );
0064
0065 PEta = xi - apply(Lapl, Eta);
0066 left = innerprod( xL, PEta, 'LR', d-1 );
0067 Yeta = tensorprod_ttemps( PEta.U{d}, left, 1 );
0068
0069 eta.dU{d} = solve_inner( L{d}, X_mid, Yeta, d );
0070
0071 for idx = d-1:-1:2
0072 X_mid.U{idx+1} = xR.U{idx+1};
0073 X_mid.U{idx} = tensorprod_ttemps(X_mid.U{idx},(G{idx}'),3);
0074
0075 Eta = tangent_to_TTeMPS( eta );
0076
0077 PEta = xi - apply(Lapl, Eta);
0078 left = innerprod( xL, PEta, 'LR', idx-1 );
0079 right = innerprod( xR, PEta, 'RL', idx+1 );
0080 res = tensorprod_ttemps( PEta.U{idx}, left, 1 );
0081 Yeta = tensorprod_ttemps( res, right, 3 );
0082
0083 eta.dU{idx} = eta.dU{idx} + solve_inner( L{idx}, X_mid, Yeta, idx );
0084 end
0085
0086
0087 X_mid.U{1+1} = xR.U{1+1};
0088 X_mid.U{1} = tensorprod_ttemps(X_mid.U{1},(G{1}'),3);
0089
0090 Eta = tangent_to_TTeMPS( eta );
0091
0092 PEta = xi - apply(Lapl, Eta);
0093 right = innerprod( xR, PEta, 'RL', 2 );
0094 Yeta = tensorprod_ttemps( PEta.U{1}, right, 3 );
0095
0096 eta.dU{1} = eta.dU{1} + solve_inner( L{1}, X_mid, Yeta, 1 );
0097
0098
0099
0100
0101
0102 eta = TTeMPS_tangent_orth( xL, xR, eta );
0103
0104
0105 end
0106
0107 function X = get_mid(xL, xR, G, idx)
0108 X = xR;
0109 X.U{1:idx-1} = xL.U{1:idx-1};
0110 if idx>1
0111 X.U{idx} = tensorprod_ttemps(X.U{idx},G{idx-1},1);
0112 end
0113 end
0114
0115 function [res,BB1,BB3] = solve_inner( L0, X, Fi, idx )
0116 n = size(L0, 1);
0117 rl = X.rank(idx);
0118 rr = X.rank(idx+1);
0119
0120 B1 = zeros( rl );
0121
0122
0123 for i = 1:idx-1
0124
0125 tmp = X;
0126 tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
0127
0128 B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
0129 end
0130
0131
0132 B2 = L0;
0133
0134 B3 = zeros( rr );
0135
0136
0137 for i = idx+1:X.order
0138 tmp = X;
0139 tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
0140
0141 B3 = B3 + innerprod( X, tmp, 'RL', idx+1);
0142 end
0143
0144
0145
0146
0147
0148 [V1,E1] = eig(B1); [V3,E3] = eig(B3);
0149 V = kron(V3,V1);
0150 EE = diag(E1)*ones(1,rr) + ones(rl,1)*diag(E3)'; E = EE(:);
0151
0152 rhs = matricize( Fi, 2 ) * V;
0153 Y = zeros(size(rhs));
0154 for i=1:length(E)
0155 Y(:,i) = (B2 + E(i)*speye(n)) \ rhs(:,i);
0156 end
0157 res = tensorize( Y*V', 2, [rl, n, rr] );
0158 end
0159